Optimization-based limiters

[valeriabarra]

Implement element-local bounds-preserving limiters. Two limiters, one quasimonotone (OP1) and one sign-preserving (OP2) to prevent negative values of ρ, are defined in the reference paper Guba et al, "Optimization-based limiters for the spectral element method", JCP 2014. These limiters use a local element reconstruction, where quasimonotone and sign-preserving properties are defined with respect to the element nodal values (Gauss-Lobatto nodes).

Each limiter defined in Guba et al. solves a constraint optimization problem completely local to each element (plus neighboring elements to determine suitable min/max bound constraints see Figure below). These are least squares optimization problems with equality and inequality constraints for the reconstructions.

To find the new time step limited quantity of interest (a tracer mixing ratio), qⁿ⁺¹, we first split a forward/explicit time step into 2 steps. (This algorithm works for any explicit RK scheme. If the bounds-preserving property is preferred, then can use a Strong Stability Preserving RK-SSP scheme.)

1. The first step finds an intermediate representation of the solution, denoted by q*, defined locally in each element (may not be continuous across elements).
2. Then replace q* with a second intermediate approximation, denoted by q**. The simplest, least accurate q** can be the element average of q*. To find a better, higher-accurate approximation we want to find the optimal closest approximation of q** to q*, which requires solving a standard quadratic program (QP), where we minimize the quadratic objective function (i.e., least-square, density-weighted distance between q* and q**), subject to linear constraints. Hence, we construct the QP, by constraining q** by the the min/max of the quantity q, qⁿₘᵢₙ and qⁿₘₐₓ, as given by the previous time step n, within the given element and its neighboring elements (i.e., over all nodal nodes in the stencil used to compute ∇q).

Rough plan and current questions:

• [x] 1. Define two different callbacks for the ODEProblem. One is a stage_callback applied at each RK stage (whose rhs function is not DSS'ed ?) and one is a step_callback applied at the end of each new (n+1) step (which includes the DSS step).
• [x] 2. Call auxiliary functions for least squares optimization problems within ODE solver stage callback
  • [x] 2.a) Inside the auxiliary function, figure out the neighboring element values for min/max bounds
  • [x] 2.b) The QPs are given by eq. (37)-(40) in the ref paper (use Quadprog.jl ?)
  • [x] 2.c) To define the QP constraints, we need to preserve the mass of the quantity local to each element. We would need to define an element-wise integral function sum(ρ)_e, for element e (or perhaps, more simply, extract the Field values over elem e, ρ_e, and just use sum(ρ_e)).
  • [x] 3. Debug
  • [x] 4. Test and verify with advection problem on the 2D sphere with cosine bells, Gaussian hills, and slotted cylinders for tracer ICs with and without hyperdiffusion
  • [x] 5. Qualitative analysis: Plot tracers in lat-long grid and compare w/ paper results
  • [x] 6. Quantitative analysis: Compare tables and convergence rates w/ paper results

[simonbyrne]

Having looked at this a bit, it looks like the Trixi team have already added a stage callback mechanism to OrdinaryDiffEq: https://trixi-framework.github.io/Trixi.jl/stable/callbacks/#Stage-callbacks (thanks @ranocha!)

From https://diffeq.sciml.ai/stable/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)

The SSP coefficients of the methods can be queried as ssp_coefficient(alg). All explicit SSP methods take two optional arguments SSPXY(stage_limiter!, step_limiter!), where stage_limiter! and step_limiter! are functions taking arguments of the form limiter!(u, integrator, p, t). Here, u is the new solution value (updated inplace) after an explicit Euler stage / the whole time step , integrator the ODE integrator, and t the current time. These limiters can be used to enforce physical constraints, e.g. the positivity preserving limiters of Zhang and Shu (Zhang, Xiangxiong, and Chi-Wang Shu. "Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments." Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 2011.).

[valeriabarra]

Thanks for these pointers! I'm going to check OrdinaryDiffEq's limiter interface and how this Zhang-Shu limiter differs from the suggested one from Guba et al. (At a very first glance. Xhang-Shu seems to be formulated for DG and FV-WENO discretizations)

[valeriabarra]

But I can see how they used the stage callbacks. Thanks

[ranocha]

Yes, I think that should work. I added this limiter interface specifically for hyperbolic PDE discretizations. All explicit SSPRK methods should support them and many additional RK methods designed for hyperbolic PDEs, too. Feel free to ping me if anything is unclear.

[valeriabarra]

Okay, having looked at the OrdinaryDiffEqs docs and source code, as well as DiffEqCallbacks.jl's docs, and from the very useful example on how to define new callbacks by the Trixi team, I think I have a clearer understanding of what we want to do.

We most certainly want to set both the stage_limiter! and step_limiter!. My current understanding from the Guba paper is that stage_limiter! will not have the DSS internal call, while the step_limiter! will. We can group them (together with additional optional callbacks, like for viz and other queries) in a CallbackSet that is then passed to the solve call. Thanks again for the pointer and to the Trixi team for writing such a useful example!

[valeriabarra]

Skeleton started in this branch valeria/opt-based-limiters

[valeriabarra]

@simonbyrne for Point 2.c) above, would you agree that it would be easier to extract values of a Field ρ over an elem e, say ρ_e, and then simply find its integral using sum(ρ_e), than overriding an element-wise sum method?

[valeriabarra]

For 2.b), would this Quadprog package be ok? I'm not sure if I can use it with our Fields, both scalar-valued and vector-valued fields

[simonbyrne]

@simonbyrne for Point 2.c) above, would you agree that it would be easier to extract values of a Field ρ over an elem e, say ρ_e, and then simply find its integral using sum(ρ_e), than overriding an element-wise sum method?

I think so: hopefully sum(slab(field, e)) should work (but I'm not sure).

[simonbyrne]

For 2.b), would this Quadprog package be ok? I'm not sure if I can use it with our Fields, both scalar-valued and vector-valued fields

Possibly? The easiest option initially is probably to pull out the underlying array (using parent). Once you get it to work, we can then try to figure out a nicer interface.

[jakebolewski]

since it is local to each element, they are solved independently and the problem is fairly small do you need a full blown ipsolver? would it be possible to construct the linear system directly and do a fixed number of newton iterations?

[valeriabarra]

446 takes care of 2.a)