Operators: discontinuous Galerkin

add_numerical_flux_interior!(fn, dydt, args...)
add_numerical_flux_interior!(ghost_exchange, fn, dydt, args...)

Add the numerical flux at the interior faces of the spectral element mesh to the mass-weighted residual dydt (WJ * ∂Y/∂t), and return dydt.

At each face node, the flux is fn(normal, argvals⁻, argvals⁺), where

• normal is the unit normal vector, pointing from the "minus" side to the "plus" side,
• argvals⁻ is the tuple of values of args on the "minus" side of the face,
• argvals⁺ is the tuple of values of args on the "plus" side of the face,

and fn returns the net flux from the "minus" side to the "plus" side, scaled by the surface Jacobian weight and subtracted from the minus side and added to the plus side. For consistency, fn must satisfy

fn(normal, argvals⁻, argvals⁺) == -fn(-normal, argvals⁺, argvals⁻)

The method with a leading ghost_exchange consumes a shared halo exchange from start_dg_ghost_exchange on distributed spaces.

See also:

• CentralNumericalFlux
• RusanovNumericalFlux

add_numerical_flux_boundary!(fn, dydt, args...)

Add the numerical flux at the domain-boundary faces of the spectral element mesh, and return dydt:

dydt -= sWJ * fn(normal, argvals⁻)

per boundary face node, where normal is the outward unit normal and argvals⁻ is the tuple of values of args at that node. dydt must be in mass-weighted residual form (WJ * ∂Y/∂t), matching add_numerical_flux_interior!. Implemented for pure 2D spectral element spaces and extruded spaces with 2D horizontal spectral elements ($sWJ$ then carries the vertical measure). No-op on domains without boundary faces (e.g. the sphere).

add_numerical_flux_boundary!(numflux, bc::HorizontalBoundaryCondition, dydt, args...)

Add the numerical flux at the domain-boundary faces using the boundary condition bc: at each boundary node the exterior state is ghost_state(bc, normal, argvals⁻), and the flux is numflux(normal, argvals⁻, argvals⁺). Returns dydt. A bc outside the HorizontalBoundaryCondition family throws.

add_lifting_flux_interior!(fn, dydt, args...)
add_lifting_flux_interior!(ghost_exchange, fn, dydt, args...)

Add symmetric face lifting terms at interior faces to dydt, and return dydt. This is the DG correction for non-conservative (gradient / curl) terms, where both sides of a face receive their own correction rather than equal-and-opposite fluxes:

dydt⁻ += sWJ * fn(n̂⁻, argvals⁻, argvals⁺)
dydt⁺ += sWJ * fn(n̂⁺, argvals⁺, argvals⁻)

with n̂⁻ = -n̂⁺ the outward unit normals. For example, the strong-form DG gradient of a scalar q is completed by fn(n̂, (q⁻,), (q⁺,)) = ((q⁺ − q⁻)/2) * n̂ (the lifting of (q* − q⁻) n̂ with a central interface value q*).

dydt must be in mass-weighted residual form (WJ * ∂Y/∂t), matching add_numerical_flux_interior!. Implemented for pure 2D spectral element spaces and for extruded spaces with 1D (plane) or 2D (e.g. cubed-sphere) horizontal spectral elements. The method with a leading ghost_exchange consumes a shared halo exchange from start_dg_ghost_exchange on distributed spaces.

lifting_correction(fn, ::Type{T}, args...)

Return the WJ-normalized DG face-lifting correction field of element type T: apply add_lifting_flux_interior! with face function fn to a zero residual on the space of args[1] and divide by WJ. The result is the correction to the corresponding element-local strong-form operator.

add_flux_differencing_divergence!(fn2pt, dydt, y)

Add the horizontal flux-differencing (split-form / FDDG) volume divergence to the mass-weighted residual dydt, and return dydt ([7], Eqs. 25–30): the collocation derivative acts on symmetric two-point fluxes along each reference direction.

fn2pt(nvec_a, nvec_b, y_a, y_b) returns the two-point flux contracted with the (non-unit) nodal metric vectors in the local orthonormal horizontal frame. It must be jointly linear in (nvec_a, nvec_b), symmetric under (nvec_a, y_a) ↔ (nvec_b, y_b), and consistent (fn2pt(n, n, y, y) == F(y)⋅n). Kinetic-energy / entropy properties are fixed by this choice (e.g. Kennedy–Gruber → KEP).

Stored in weak-equivalent form (strong flux-differencing plus one-sided boundary lifts), so it replaces dydt = hwdiv(F) * (-WJ) and composes with add_numerical_flux_interior! to give the FDDG SAT $F^* - F(y^-)⋅n̂$. SBP telescoping gives local conservation; global conservation follows from antisymmetric interface fluxes.

Supports pure 2D spectral elements and extruded spaces with 2D horizontal elements.

add_sipg_laplacian_flux_interior!(dydt, q, G, κ, τ)

Add SIPGLaplacianFlux at interior faces to a mass-weighted (WJ) Laplacian residual. τ may be a scalar or a Field. The method with a leading ghost_exchange reuses a start_dg_ghost_exchange handle started on (q, G, κ, τ).

sipg_laplacian_tendency(q, ρ_weight, κ, τ)

Horizontal Laplacian κ ∇⋅(ρ_weight ∇q) on a DG space, or κ ∇²q when ρ_weight === nothing. Neighbouring elements are coupled by SIPGLaplacianFlux, with τ from sipg_penalty_parameter.

Allocates the result and two scratch fields; sipg_laplacian_tendency! takes all three as arguments.

sipg_laplacian_tendency!(out, G_uv, R_uv, q, ρ_weight, κ, τ)

In-place sipg_laplacian_tendency: writes into out, using G_uv and R_uv as scratch. None of the three may alias q or ρ_weight. Returns out.

sipg_penalty_parameter(κ, space; weight = nothing)

Penalty field τ = κ w (2Nq − 1)^2 / h for the DG Laplacian, where h is the node spacing and w is weight (1 when nothing).

τ sets how hard neighbouring elements are pushed to agree at their shared face. Too small and the Laplacian stops damping — it can amplify the grid-scale noise it exists to remove. It is a Field rather than one number because both inputs vary over the mesh: the node spacing by 1.4x on a cubed sphere, and weight by however much the caller's does. A face uses the larger of its two sides.

Allocates; sipg_penalty_parameter! writes into a field you own.

sipg_penalty_parameter!(τ, κ; weight = nothing)

In-place sipg_penalty_parameter: fills τ, whose space supplies the geometry, and returns it.

start_dg_ghost_exchange(args...)
start_dg_ghost_exchange(space, args...)

Start the ghost-face halo exchange of the DG face operators once, to be shared by several operator calls in the same tendency evaluation:

ex = Operators.start_dg_ghost_exchange(y)
# ... element-local volume terms (the exchange overlaps them) ...
Operators.add_numerical_flux_interior!(ex, numflux, dydt, y)
Operators.add_lifting_flux_interior!(ex, lift, dydt2, y)

Without the leading handle each operator performs its own exchange, so an RHS that applies several face operators to the same state sends every halo message once per operator; a shared exchange sends each once. The space is taken from the first Field in args, or passed as a leading argument.

Contract: every operator receiving the handle must be called with the same args (the same fields, in the same order). Mismatched argument types or counts throw; two same-typed fields swapped in place cannot be detected. The handle covers one exchange round — start a fresh one whenever an argument's values change (e.g. each stage), and pass a started handle to at least one operator call before starting the next round: the underlying buffers are shared per argument type and position on a space (also across distinct fields of the same type), and starting a round while another is in flight throws. Returns a no-op handle on single-process contexts, so calling code needs no distributed-vs-single branch.

DGGhostExchange

Handle for one round of the DG ghost-face halo exchange, shared across face operators — see start_dg_ghost_exchange. The first operator that consumes the handle completes the exchange (ClimaComms.finish); later consumers read the same recv strips.

tendency_completion(dydt; numflux, boundary_numflux = nothing)

Build the interface-coupling configuration for the tendency field dydt (or any field with the tendency's space and value type), selected by dispatch on Grids.discretization: a DSSCompletion on Grids.CG() spaces, a NumericalFluxCompletion on Grids.DG() ones. Call once at model setup; the discretization choice then lives in the space, and the tendency code is shared:

completion = Operators.tendency_completion(dydt; numflux)
# ... each RHS evaluation:
@. dydt = -wdiv(physical_flux(y, params))       # element-local weak form
Operators.complete_tendency!(completion, dydt, y, params)

numflux(normal, argvals⁻, argvals⁺) is the DG interface flux (see add_numerical_flux_interior! for the contract); it must be built from the same physical flux as the weak volume term. It is required on DG spaces and unused on CG spaces, so models can pass it unconditionally. boundary_numflux(normal, argvals⁻) is the one-sided flux at domain-boundary faces (see add_numerical_flux_boundary!); without it, boundary faces contribute nothing to the DG tendency (a zero-flux closure). CG boundary conditions are imposed by the operators themselves, not here.

dydt may be a Field on either discretization, or a FieldVector on a continuous one, where the completion is a single batched Spaces.weighted_dss! over all components. A FieldVector is not supported on a discontinuous space: the interface flux is evaluated on the whole state at a face node, which requires the state to be one Field with a composite (e.g. NamedTuple) eltype rather than a collection of Fields.

complete_tendency!(completion, dydt, args...)

Complete the element-local weak-form tendency dydt across element interfaces, using the tendency_completion built for its space. The contract on dydt is the CG weak-form convention (e.g. -wdivₕ(F) for a flux-form equation, with no WJ weighting); args are the arguments of the completion's flux functions, the model state first.

On a DSSCompletion this is Spaces.weighted_dss!(dydt, buffer) (args unused). On a NumericalFluxCompletion it is the mass-weighted DG surface term:

dydt *= WJ
add_numerical_flux_interior!(numflux, dydt, args...)
add_numerical_flux_boundary!(boundary_numflux, dydt, args...)  # if given
dydt /= WJ

Mutates dydt in place and returns it.

AbstractTendencyCompletion

Supertype of DSSCompletion and NumericalFluxCompletion, the two ways complete_tendency! couples element-local weak-form tendencies across element interfaces. Construct with tendency_completion, which selects the subtype from the discretization of the space (dispatch on Grids.discretization).

DSSCompletion(buffer)

AbstractTendencyCompletion for continuous (CG) spaces: complete_tendency! applies Spaces.weighted_dss! with the stored buffer, projecting the element-local tendency onto the continuous space. Constructed by tendency_completion.

NumericalFluxCompletion(numflux, boundary_numflux)

AbstractTendencyCompletion for discontinuous (DG) spaces: complete_tendency! weights the tendency by WJ, accumulates the interface flux numflux with add_numerical_flux_interior! (and, when boundary_numflux is not nothing, the one-sided boundary flux with add_numerical_flux_boundary!), and unweights. Constructed by tendency_completion.

cartesian_tensor_divergence(T, completion, faceargs...)

Weak-form horizontal divergence ∇ₕ·T of a rank-2 flux tensor field T (e.g. the momentum flux ρu⊗u). On a curved space ∇·T carries a Christoffel connection term on T's momentum axis that the weak Divergence omits; rotating that axis into the global Cartesian basis, where the Christoffel symbols vanish, makes the omission exact. So the momentum (second) axis is rotated to Cartesian (Geometry.CartesianTensor), Divergence is applied, and the result rotated back to the local frame (Geometry.LocalVector); on a plane (CartesianGlobalGeometry) both rotations are the identity and are skipped.

Like the spectral Divergence it is built from, this differentiates along the two horizontal directions only, so on an extruded space the vertical flux divergence is a separate term. That term drops the same connection terms, so rotate its flux tensor to Cartesian as well before differencing it, as in the example below.

The result is a 3D Geometry.UVWVector, whose w component on the sphere is the curvature term — -|u|²/R for u⊗u under solid-body rotation, the same order as the tangential components. A model carrying horizontal-only momentum has to drop it, e.g. with Geometry.project(Geometry.UVAxis(), ...).

T's momentum (second) axis must be a local orthonormal axis; a UVAxis is read as a UVWAxis with w == 0. On a discontinuous space the transport (first) axis must be orthonormal too, because numflux contracts it against the local orthonormal face normal.

completion (from tendency_completion) couples element interfaces — DSS on a CG space, the numerical flux on a DG one, the same CG↔DG switch as complete_tendency!. On a DG space faceargs... are the trailing arguments of the completion's numflux (the Cartesian tensor is its first face argument, and the flux must return the momentum vector it produces from that tensor, whose UVW components are global Cartesian); on a CG space they are unused, as is the completion's DSS buffer — the buffer this needs is one for its own UVWVector result, which it owns — so a model can pass the completion it built for its own tendency.

Allocates the result, and the Cartesian-tensor scratch wherever the rotation is not the identity; cartesian_tensor_divergence! takes both buffers and is allocation-free.

Examples

The whole ∇·T on an extruded space, the vertical flux divergence carrying the same rotation. ᶜT is the flux tensor on cell centers and ᶠT the same flux on faces:

geom = Spaces.global_geometry(axes(ᶜT))
ᶜcoords = Fields.coordinate_field(axes(ᶜT))
ᶠcoords = Fields.coordinate_field(axes(ᶠT))

ᶜdivₕ = Operators.cartesian_tensor_divergence(ᶜT, completion)

ᶠTc = @. Geometry.CartesianTensor(ᶠT, geom, ᶠcoords)
ᶜdivᵥ = Operators.DivergenceF2C().(ᶠTc)

ᶜdivT = @. ᶜdivₕ + Geometry.LocalVector(ᶜdivᵥ, geom, ᶜcoords)

Differencing the local-frame ᶠT leaves the connection terms of the vertical derivative in place, an error that on a topographic sphere runs larger than the divergence itself (1.8x its peak at Ne = 4, GLL{4}, 10 levels). The two terms are rotated back one at a time here, which agrees with rotating their sum to roundoff, the rotation being linear.

References

• [53]: Cartesian momentum components remove the connection terms.

cartesian_tensor_divergence!(out, Tc, T, completion, faceargs...)

In-place cartesian_tensor_divergence: writes the divergence into the local-frame vector field out, using the rank-2 field Tc as scratch for the Cartesian-rotated flux. Neither out nor Tc may alias T. Tc is untouched on a plane, where the rotation is the identity and T is differentiated directly. Returns out.

Rotating the momentum axis widens it to the full 3D UVWAxis, so similar(T) is a wide enough scratch only when T's momentum axis already is, and a narrower Tc raises an error. Given a UVAxis momentum, allocate Tc from the rotation, as the allocating form does, or promote the momentum vector before forming the flux.

Allocation-free given the two buffers.

AbstractNumericalFlux

Supertype of the numerical flux functors of the DG face operators. A subtype is callable as numflux(normal, argvals⁻, argvals⁺); see add_numerical_flux_interior! for the contract. Subtypes: CentralNumericalFlux, RusanovNumericalFlux, and SIPGLaplacianFlux.

CentralNumericalFlux(fluxfn)

Interface numerical flux ((F⁻ + F⁺) / 2)' * normal, the average of the physical flux fluxfn(args...) on the two sides of a face. It adds no dissipation, so an advection-dominated DG run with it is unstable; use RusanovNumericalFlux for a dissipative alternative built from the same fluxfn. Callable as numflux(normal, argvals⁻, argvals⁺), the contract of add_numerical_flux_interior!.

Examples

numflux = Operators.CentralNumericalFlux(physical_flux)
completion = Operators.tendency_completion(dydt; numflux)

RusanovNumericalFlux(fluxfn, wavespeedfn)

Interface numerical flux ((F⁻ + F⁺) / 2)' * normal + (λ / 2) * (y⁻ - y⁺): the central flux of fluxfn(args...) plus a jump penalty with λ = max(wavespeedfn(args⁻...), wavespeedfn(args⁺...)), an upper bound on the normal signal speed. The penalty is the dissipation that keeps a DG transport scheme stable; wavespeedfn may overestimate the speed at the cost of extra dissipation. Callable as numflux(normal, argvals⁻, argvals⁺), the contract of add_numerical_flux_interior!; the first argument value on each side is the state y.

Examples

sw_flux(y, p) = (; ρ = y.ρu, ρu = (y.ρu ⊗ (y.ρu / y.ρ)) + (p.g * y.ρ^2 / 2) * I)
sw_wavespeed(y, p) = sqrt(p.g * y.ρ) + norm(y.ρu / y.ρ)
numflux = Operators.RusanovNumericalFlux(sw_flux, sw_wavespeed)
completion = Operators.tendency_completion(dydt; numflux)

SIPGLaplacianFlux()

Interface flux for the DG Laplacian, used with add_numerical_flux_interior!. Each side supplies (q, G, κ, τ), where G is the gradient in the same basis as the face normal (usually Geometry.UVVector) and τ is a scalar or a Field. Returns the gradient averaged across the face plus a penalty on the jump in q: -{{κG}}·n̂ + max(τ⁻, τ⁺) * (q⁻ - q⁺).

Two of the three face terms; the third rides along in the volume divergence of sipg_laplacian_tendency!, which is what makes the operator symmetric.

central_gradient_lift(normal, (q⁻,), (q⁺,))

Return the symmetric central lifting that completes the strong-form DG gradient of a scalar: each side adds $(q^* - q_{side}) n̂_{side}$ with central $q^*$, i.e. $((q⁺ - q⁻)/2)\,n̂$ on the minus side. Use with add_lifting_flux_interior! / lifting_correction.

central_curl3_lift(normal, (u⁻, v⁻), (u⁺, v⁺))

Return the central lifting for the radial component of the horizontal curl: $r̂ ⋅ (n̂ × (u^* - u_{side}))$ from the tangential jumps of the orthonormal velocity components (u, v).

jump_penalty_lift(normal, (q⁻, λ⁻), (q⁺, λ⁺))

Return the λ-scaled interface penalty max(λ⁻, λ⁺) / 2 * (q⁺ - q⁻): each side relaxes toward its neighbor at rate $\max(λ⁻, λ⁺)/2$.

HorizontalBoundaryCondition <: AbstractBoundaryCondition

Supertype for the boundary conditions of the horizontal DG numerical-flux operators (see add_numerical_flux_boundary!), e.g. ReflectingWallBC. Subtypes should define ghost_state.

ghost_state(bc::HorizontalBoundaryCondition, normal, argvals⁻)

Construct the exterior-side argument tuple for the boundary condition bc.

Returns

A tuple with the same length as argvals⁻, replacing only the prognostic state argvals⁻[1] with the ghost state; remaining arguments (e.g. equation parameters, coordinates) are forwarded unchanged. The fallback method throws for boundary conditions without a ghost_state method.

PeriodicBC <: HorizontalBoundaryCondition

Periodic boundary condition for the horizontal DG numerical-flux operators, handled by the topology (no ghost state needed).

ReflectingWallBC <: HorizontalBoundaryCondition

Reflecting-wall (no-normal-flow) boundary condition for the horizontal DG numerical-flux operators: its ghost_state reflects the normal momentum component and preserves the other prognostic fields.