Limiters

AbstractLimiter

Supertype for all limiters.

Subtypes:

• QuasiMonotoneLimiter: horizontal quasi-monotone flux limiter for spectral element advection.
• PositivityLimiter: Zhang-Shu mean-preserving positivity limiter for conserved DG states.
• VerticalMassBorrowingLimiter: vertical mass-borrowing limiter that removes negative tracer mass.

Subtypes implement apply_limiter!.

QuasiMonotoneLimiter

Quasi-monotone limiter for tracer densities on spectral element spaces, after the OP1 limiter of [14], eqs. (37)-(40). Quasi-monotone means monotone with respect to the spectral element nodal values.

For each element, the limiter finds the nodal field closest, in the mass-weighted l2 norm, to the input field that satisfies min/max bounds on the concentration q = ρq / ρ. As in HOMME, it clips the nodal values that violate the bounds to the nearest bound, which changes the tracer mass of the element, and then redistributes the mass change over the nodes that are not at a bound so that the l2 error is smallest. Redistribution can violate the bounds again, so the two steps are iterated until abs(Δtracer_mass) <= rtol * tracer_mass or Nq^2 iterations have been done. The optimization is local to each element; the neighbors enter only through the bounds.

Fields

• q_bounds: min and max of q in each element.
• q_bounds_nbr: min and max of q over each element and its neighbors.
• ghost_buffer: buffer for exchanging q_bounds with neighboring processes.
• rtol: relative tolerance for the tracer mass change per element [-].
• convergence_stats: LimiterConvergenceStats (or NoConvergenceStats) accumulated by apply_limiter!.

Constructor

QuasiMonotoneLimiter(
    ρq::Field;
    rtol = eps(eltype(parent(ρq))),
    convergence_stats = LimiterConvergenceStats{eltype(parent(ρq))}(),
)

Create a limiter for the tracer density field ρq, where q is the tracer concentration per unit mass; ρq can be a scalar-valued or a struct-valued Field. convergence_stats records the number of apply_limiter! calls in which some element failed to converge, the largest relative mass error, and the smallest tracer mass; print_convergence_stats(limiter) prints them.

Examples

Call compute_bounds! on the fields at the start of the step, then apply_limiter! on the updated fields:

limiter = QuasiMonotoneLimiter(ρq)
compute_bounds!(limiter, ρq, ρ)
# ... advance ρq and ρ ...
apply_limiter!(ρq, ρ, limiter)

VerticalMassBorrowingLimiter(q_min)

Vertical-only mass-borrowing limiter.

The limiter borrows tracer mass from adjacent lower layers to raise each layer's tracer value to at least its minimum. It conserves the total tracer mass in the column.

q_min is a tuple with one minimum tracer value per component of the limited field.

At level k, the limiter first borrows mass from layer k+1 (the lower level). If the mass in layer k+1 is not sufficient, it borrows from layer k+2, and so on down to the bottom layer. If the tracer mass in the bottom layer goes below the minimum, the limiter repeats the process from the bottom to the top. This makes the limiter work for any shape of mass profile.

Examples

ρ = fill(1.0, space)
q = fill((a = 0.1, b = 0.1), space)
limiter = VerticalMassBorrowingLimiter((0.0, 0.0))
Limiters.apply_limiter!(q, ρ, limiter)

Adapted from the E3SM mass borrower; see [15].

PositivityLimiter{FT, G, F} <: AbstractLimiter

Mean-preserving positivity limiter of [16] for the conserved state of a DG discretization, applied with apply_limiter!.

In each element, with Ū the WJ-weighted mean:

1. Each field with a linear floor is scaled toward its own mean, x_j ← x̄ + θ₁ (x_j − x̄), by the largest θ₁ ∈ [0, 1] that keeps its nodal minimum at the floor (per component for NamedTuple-valued fields). Fields without a floor are untouched.

2. If g is given, the whole conserved vector is then scaled toward the mean, U_j ← Ū + θ₂ (U_j − Ū), by the largest θ₂ ∈ [0, 1] with g(U_j) ≥ g_min at every node (bisection, maxiter steps per node).

Both steps are convex combinations with the element mean, so every element mean is preserved exactly.

FT is the floating-point type of the limited states, G the type of g, and F the type of floors.

Fields

• g: Concave constraint g(U, aux), or nothing to skip step 2. U is one node's state as a NamedTuple with the names of the limited states, and aux the node's value of the auxiliary field.
• floors: NamedTuple of linear floors keyed by state names, in the units of each state. A floor applies to every component of a NamedTuple-valued field.
• g_min: Lower bound of the nonlinear constraint g(U, aux) ≥ g_min, in the units of g.
• maxiter: Bisection steps per node for θ₂ [-]; θ₂ is accurate to 2^-maxiter and always on the admissible side.

Constructor

PositivityLimiter(FT; floors = (;), g = nothing, g_min = 0, maxiter = 10)

Each keyword sets the field of the same name; floors and g_min are converted to FT.

Examples

# Shallow water: depth floor only.
limiter = PositivityLimiter(Float64; floors = (; h = 0))

# Compressible Euler: density and tracer floors, and a pressure floor, with
# the geopotential Φ as the auxiliary field.
pressure(U, Φ) = 0.4 * (U.ρe - (U.ρu1^2 + U.ρu2^2 + U.ρu3^2) / 2U.ρ - U.ρ * Φ)
limiter = PositivityLimiter(
    Float64;
    floors = (; ρ = 1e-6, ρq = 0),
    g = pressure,
    g_min = 1e-3,
)

Notes

• The limiter only redistributes within an element, so it cannot fix a bad mean. Such an element is set to its mean (θ = 0) and stays below the floor. With an SSP time integrator, the mean stays admissible under the CFL condition of [16].
• The bounds hold at the nodes only; the interpolating polynomial may still dip below a floor between nodes.
• g must be concave in U (as pressure and depth are), so that the admissible θs at a node form an interval [0, θ*] and bisection finds θ*.
• States may hold ForwardDiff.Dual values. Derivatives are exact except where a floor switches on or off; the bisected θ₂ is differentiated implicitly through g(U(θ₂)) = g_min.

See also apply_limiter! and the Zhang–Shu positivity limiter (DG) how-to.

compute_bounds!(limiter::QuasiMonotoneLimiter, ρq::Field, ρ::Field)

Compute the bounds on the tracer concentration q = ρq / ρ that apply_limiter! enforces, from the tracer density ρq and the density ρ, and store them in limiter.q_bounds_nbr.

The steps are:

1. compute_element_bounds! computes the min and max of q in each element.
2. If distributed, start the ghost exchange of the element bounds.
3. compute_neighbor_bounds_local! widens the bounds with the local neighbors.
4. If distributed, complete the ghost exchange and widen the bounds with the ghost neighbors in compute_neighbor_bounds_ghost!.

apply_limiter!(ρq, ρ, limiter::QuasiMonotoneLimiter; warn = true)

Limit the tracer density ρq in place so that the concentration q = ρq / ρ at each node lies within the bounds computed by compute_bounds!, while preserving the tracer mass of each element up to the relative tolerance limiter.rtol.

Each element is processed by apply_limit_slab!, with the density ρ times the quadrature weights as the weights of the least-squares redistribution. When some element fails to converge, limiter.convergence_stats is updated. On CPU devices with warn = true, the accumulated convergence statistics are printed with @warn after each call.

Return ρq on CPU devices and nothing on CUDA devices.

apply_limiter!(states::NamedTuple, aux, limiter::PositivityLimiter)

Apply the PositivityLimiter limiter to states, element by element. Mutates the fields of states; returns nothing.

Arguments

• states: NamedTuple of the conserved Fields on one space, e.g. (; ρ, ρe, ρu1, ρu2, ρu3), with the limiter's eltype. Its names must include every key of limiter.floors, and are the names limiter.g reads. Field elements must be scalars or flat NamedTuples of scalars (several tracers in one field); pass vector quantities such as momentum as separate scalar component fields.
• aux: Field of non-conserved data passed unscaled to limiter.g, or nothing; e.g. the geopotential [m²/s²].
• limiter: The PositivityLimiter.

Examples

# h, hu1, hu2: depth and momentum `Field`s on one spectral-element space.
limiter = PositivityLimiter(Float64; floors = (; h = 0))
apply_limiter!((; h, hu1, hu2), nothing, limiter)

See Zhang–Shu positivity limiter (DG) for a runnable shallow-water example and the compressible Euler setup.

apply_limiter!(q::Fields.Field, ρ::Fields.Field, lim::VerticalMassBorrowingLimiter)

Apply the vertical mass-borrowing limiter lim to the tracer field q in place, given the density field ρ.

Each component of q is limited column by column with the corresponding entry of lim.q_min, using the cell volume from the local geometry of ρ. Returns nothing.

print_convergence_stats(limiter::QuasiMonotoneLimiter)
print_convergence_stats(io::IO, stats)

Print the convergence statistics accumulated by the Limiters.QuasiMonotoneLimiter limiter (to stdout), or the statistics object stats (to io): the number of elements in which the limiter iteration did not converge, the maximum relative error, and the minimum tracer mass. Nothing is printed if the limiter was constructed with convergence_stats = NoConvergenceStats().