Operators

grad = Gradient()
grad.(f)

Compute the (strong) gradient of f on each element, returning a CovariantVector-field.

The $i$th covariant component of the gradient is the partial derivative with respect to the reference element:

\[(\nabla f)_i = \frac{\partial f}{\partial \xi^i}\]

Discretely, this can be written in matrix form as

\[D_i f\]

where $D_i$ is the derivative matrix along the $i$th dimension.

References

• [1], equation 16

div = Divergence()
div.(u)

Computes the per-element spectral (strong) divergence of a vector field $u$.

The divergence of a vector field $u$ is defined as

\[\nabla \cdot u = \sum_i \frac{1}{J} \frac{\partial (J u^i)}{\partial \xi^i}\]

where $J$ is the Jacobian determinant, $u^i$ is the $i$th contravariant component of $u$.

This is discretized by

\[\sum_i I \left\{\frac{1}{J} \frac{\partial (I\{J u^i\})}{\partial \xi^i} \right\}\]

where $I\{x\}$ is the interpolation operator that projects to the unique polynomial interpolating $x$ at the quadrature points. In matrix form, this can be written as

\[J^{-1} \sum_i D_i J u^i\]

where $D_i$ is the derivative matrix along the $i$th dimension

References

• [1], equation 15

split_div = SplitDivergence()
split_div.(ρu, ψ)

Computes the divergence of the product ρu * ψ using a split-form (entropy-stable) discretization.

This operator is designed for the advection of scalar quantities in conservation laws (e.g., thermodynamic variables or tracers). By evaluating the divergence using a specific averaging of the conservative and advective forms, this formulation cancels aliasing errors that arise from the product of two spectrally variable fields, thereby inhibiting the growth of quadratic instabilities (such as cold temperature spikes) without requiring hyperviscosity.

Arguments

• ρu: The transport vector field, typically the mass flux. It must be a vector quantity (e.g., Geometry.Contravariant12Vector).
• ψ: The specific scalar quantity to be advected (e.g., specific total energy $e_{tot}$ or specific humidity $q_{tot}$).

Mathematical Formulation

Continuous

The split form of the divergence operator is defined as the arithmetic mean of the conservative and advective forms:

\[\nabla \cdot (\rho \mathbf{u} \psi)|_\textrm{split} = \frac{1}{2} \nabla \cdot (\rho \mathbf{u} \psi) + \frac{1}{2} \left( \psi \nabla \cdot (\rho \mathbf{u}) + \rho \mathbf{u} \cdot \nabla \psi \right)\]

Discrete

The discretized split operator is equivalent to using the strong formulation of the gradient operator and the weak formulation of the divergence operator:

\[\textrm{split_div}(\rho \mathbf{u}, \psi) = \frac{1}{2} \textrm{wdiv}(\rho \mathbf{u} \psi) + \frac{1}{2} \left( \psi \textrm{wdiv}(\rho \mathbf{u}) + \rho \mathbf{u} \cdot \textrm{grad}(\psi) \right)\]

Swapping the weak and strong formulations in the last two terms also results in the same operator. The discrete form of the divergence theorem, which stems from the generalized summation-by-parts (SBP) property, guarantees that the integral of the first term vanishes,

\[\int_\Omega \textrm{wdiv}(\rho \mathbf{u} \psi) dV = 0\]

while the integrals of the other two terms cancel out,

\[\int_\Omega \psi \textrm{wdiv}(\rho \mathbf{u}) dV = -\int_\Omega \rho \mathbf{u} \cdot \textrm{grad}(\psi) dV\]

So, this discretization ensures that the split operator conserves the integral of $\rho \mathbf{u} \psi$.

Two-Point

A more compact representation of the discretized operator can be obtained with the symmetric two-point flux, whose values in one dimension are

\[(F^1)_{ij} = \frac{1}{2} (\rho_i J_i (u^1)_i + \rho_j J_j (u^1)_j) (\psi_i + \psi_j)\]

With $D$ denoting the spectral derivative matrix, the split operator in one dimension can be expressed as

\[\textrm{split_div}(\rho \mathbf{u}, \psi)_i = \frac{1}{J_i} \sum_{j \neq i} D_{ij} (F^1)_{ij}\]

In two dimensions, $F^1$ and the analogous quantity $F^2$ provide a similar expression for the split divergence, with the one-dimensional operator applied sequentially along each dimension.

Properties

1. Conservation: The split operator conserves $\rho \mathbf{u} \psi$
2. Consistency: If $\psi = 1$, the split operator degenerates to the weak formulation of $\nabla \cdot \rho \mathbf{u}$ (mass continuity)
3. Complexity: The split operator has the same $O(N^2)$ complexity per element as the strong and weak operators, but needs twice as many operations

References

• Fisher, T. C., & Carpenter, M. H. (2013). High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains. Journal of Computational Physics, 252, 518-557. https://doi.org/10.1016/j.jcp.2013.06.014
• Gassner, G. J. (2013). A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. SIAM Journal on Scientific Computing, 35, A1233-A1253. https://doi.org/10.1137/120890144

wdiv = WeakDivergence()
wdiv.(u)

Computes the "weak divergence" of a vector field u.

This is defined as the scalar field $\theta \in \mathcal{V}_0$ such that for all $\phi\in \mathcal{V}_0$

\[\int_\Omega \phi \theta \, d \Omega = - \int_\Omega (\nabla \phi) \cdot u \,d \Omega\]

where $\mathcal{V}_0$ is the space of $u$.

This arises as the contribution of the volume integral after applying integration by parts to the weak form expression of the divergence

\[\int_\Omega \phi (\nabla \cdot u) \, d \Omega = - \int_\Omega (\nabla \phi) \cdot u \,d \Omega + \oint_{\partial \Omega} \phi (u \cdot n) \,d \sigma\]

It can be written in matrix form as

\[ϕ^\top WJ θ = - \sum_i (D_i ϕ)^\top WJ u^i\]

which reduces to

\[θ = -(WJ)^{-1} \sum_i D_i^\top WJ u^i\]

where

• $J$ is the diagonal Jacobian matrix
• $W$ is the diagonal matrix of quadrature weights
• $D_i$ is the derivative matrix along the $i$th dimension

wgrad = WeakGradient()
wgrad.(f)

Compute the "weak gradient" of f on each element.

This is defined as the the vector field $\theta \in \mathcal{V}_0$ such that for all $\phi \in \mathcal{V}_0$

\[\int_\Omega \phi \cdot \theta \, d \Omega = - \int_\Omega (\nabla \cdot \phi) f \, d\Omega\]

where $\mathcal{V}_0$ is the space of $f$.

This arises from the contribution of the volume integral after by applying integration by parts to the weak form expression of the gradient

\[\int_\Omega \phi \cdot (\nabla f) \, d \Omega = - \int_\Omega f (\nabla \cdot \phi) \, d\Omega + \oint_{\partial \Omega} f (\phi \cdot n) \, d \sigma\]

In matrix form, this becomes

\[{\phi^i}^\top W J \theta_i = - ( J^{-1} D_i J \phi^i )^\top W J f\]

which reduces to

\[\theta_i = -W^{-1} D_i^\top W f\]

where $D_i$ is the derivative matrix along the $i$th dimension.

curl = Curl()
curl.(u)

Computes the per-element spectral (strong) curl of a covariant vector field $u$.

Note: The vector field $u$ needs to be excliclty converted to a CovaraintVector, as then the Curl is independent of the local metric tensor.

The curl of a vector field $u$ is a vector field with contravariant components

\[(\nabla \times u)^i = \frac{1}{J} \sum_{jk} \epsilon^{ijk} \frac{\partial u_k}{\partial \xi^j}\]

where $J$ is the Jacobian determinant, $u_k$ is the $k$th covariant component of $u$, and $\epsilon^{ijk}$ are the Levi-Civita symbols. In other words

\[\begin{bmatrix} (\nabla \times u)^1 \\ (\nabla \times u)^2 \\ (\nabla \times u)^3 \end{bmatrix} = \frac{1}{J} \begin{bmatrix} \frac{\partial u_3}{\partial \xi^2} - \frac{\partial u_2}{\partial \xi^3} \\ \frac{\partial u_1}{\partial \xi^3} - \frac{\partial u_3}{\partial \xi^1} \\ \frac{\partial u_2}{\partial \xi^1} - \frac{\partial u_1}{\partial \xi^2} \end{bmatrix}\]

In matrix form, this becomes

\[\epsilon^{ijk} J^{-1} D_j u_k\]

Note that unused dimensions will be dropped: e.g. the 2D curl of a Covariant12Vector-field will return a Contravariant3Vector.

References

• [1], equation 17

wcurl = WeakCurl()
wcurl.(u)

Computes the "weak curl" on each element of a covariant vector field u.

Note: The vector field $u$ needs to be excliclty converted to a CovaraintVector, as then the WeakCurl is independent of the local metric tensor.

This is defined as the vector field $\theta \in \mathcal{V}_0$ such that for all $\phi \in \mathcal{V}_0$

\[\int_\Omega \phi \cdot \theta \, d \Omega = \int_\Omega (\nabla \times \phi) \cdot u \,d \Omega\]

where $\mathcal{V}_0$ is the space of $f$.

This arises from the contribution of the volume integral after by applying integration by parts to the weak form expression of the curl

\[\int_\Omega \phi \cdot (\nabla \times u) \,d\Omega = \int_\Omega (\nabla \times \phi) \cdot u \,d \Omega - \oint_{\partial \Omega} (\phi \times u) \cdot n \,d\sigma\]

In matrix form, this becomes

\[{\phi_i}^\top W J \theta^i = (J^{-1} \epsilon^{kji} D_j \phi_i)^\top W J u_k\]

which, by using the anti-symmetry of the Levi-Civita symbol, reduces to

\[\theta^i = - \epsilon^{ijk} (WJ)^{-1} D_j^\top W u_k\]

FormType

Supertype of the singleton types StrongForm and WeakForm, which distinguish the variational form of a spectral element operator. The two variants share the same interior computation; they differ only in the applied derivative matrix (D vs. its integration-by-parts counterpart -Dᵀ) and in the weights that multiply the argument and divide the result. Weak variants are defined as aliases; e.g., WeakDivergence() is a Divergence{WeakForm}.

StrongForm()

The FormType of an operator that discretizes a derivative directly at the quadrature points (e.g. Divergence, Gradient, and Curl).

WeakForm()

The FormType of an operator that discretizes the volume-integral contribution of the corresponding weak-form expression, obtained after integration by parts (e.g. WeakDivergence, WeakGradient, and WeakCurl).

i = Interpolate(space)
i.(f)

Interpolates f to the space. If space has equal or higher polynomial degree as the space of f, this is exact, otherwise it will be lossy.

In matrix form, it is the linear operator

\[I = \bigotimes_i I_i\]

where $I_i$ is the barycentric interpolation matrix in the $i$th dimension.

See also Restrict.

r = Restrict(space)
r.(f)

Computes the projection of a field f on $\mathcal{V}_0$ to a lower degree polynomial space space ($\mathcal{V}_0^*$). space must be on the same topology as the space of f, but have a lower polynomial degree.

It is defined as the field $\theta \in \mathcal{V}_0^*$ such that for all $\phi \in \mathcal{V}_0^*$

\[\int_\Omega \phi \theta \,d\Omega = \int_\Omega \phi f \,d\Omega\]

In matrix form, this is

\[\phi^\top W^* J^* \theta = (I \phi)^\top WJ f\]

where $W^*$ and $J^*$ are the quadrature weights and Jacobian determinant of $\mathcal{V}_0^*$, and $I$ is the interpolation operator (see Interpolate) from $\mathcal{V}_0^*$ to $\mathcal{V}_0$. This reduces to

\[\theta = (W^* J^*)^{-1} I^\top WJ f\]

FiniteDifferenceOperator

An abstract type for finite difference operators. Instances of this should define:

• return_eltype
• return_space
• stencil_interior_width
• stencil_interior

See also AbstractBoundaryCondition for how to define the boundaries.

I = InterpolateC2F(;boundaries..)
I.(x)

Interpolate a center-valued field x to faces, using the stencil

\[I(x)[i] = \frac{1}{2} (x[i+\tfrac{1}{2}] + x[i-\tfrac{1}{2}])\]

Supported boundary conditions are:

• SetValue(x₀): set the value at the boundary face to be x₀. On the left boundary the stencil is

\[I(x)[\tfrac{1}{2}] = x₀\]

• Extrapolate: use the closest interior point as the boundary value. At the left boundary the stencil is

\[I(x)[\tfrac{1}{2}] = x[1]\]

InterpolateF2C()

Interpolate from face to center mesh. No boundary conditions are required (or supported).

WI = WeightedInterpolateC2F(; boundaries)
WI.(w, x)

Interpolate a center-valued field x to faces, weighted by a center-valued field w, using the stencil

\[WI(w, x)[i] = \frac{ w[i+\tfrac{1}{2}] x[i+\tfrac{1}{2}] + w[i-\tfrac{1}{2}] x[i-\tfrac{1}{2}]) }{ w[i+\tfrac{1}{2}] + w[i-\tfrac{1}{2}] }\]

Supported boundary conditions are:

• SetValue(val): set the value at the boundary face to be val.
• Extrapolate: use the closest interior point as the boundary value.

These have the same stencil as in InterpolateC2F.

WI = WeightedInterpolateF2C(; boundaries)
WI.(w, x)

Interpolate a face-valued field x to centers, weighted by a face-valued field w, using the stencil

\[WI(w, x)[i] = \frac{ w[i+\tfrac{1}{2}] x[i+\tfrac{1}{2}] + w[i-\tfrac{1}{2}] x[i-\tfrac{1}{2}]) }{ w[i+\tfrac{1}{2}] + w[i-\tfrac{1}{2}] }\]

No boundary conditions are required (or supported)

AdvectionOperator

An abstract type for advection operators. As of now, advection operators that do the following are supported:

Given a face-valued velocity field v and a center-valued field x, for each face i the advection operator computes a function of the form f(v[i-1], v[i], v[i+1], x[i-3/2], x[i-1/2], x[i+1/2], x[i+3/2]) or f(v[i], x[i-3/2], x[i-1/2], x[i+1/2], x[i+3/2]) and returns a contravariant3 component. On non-periodic domains, all faces are treated like interior faces, padding out-of-range stencil points with ghost values (on periodic domains, indices wrap around instead):

• The out-of-range values of the advected field are padded with the Extrapolate boundary condition for that boundary: every ghost point the stencil reaches takes the value extrapolated from the in-range interior points of the stencil (the extrapolation order is reduced at the boundary face itself, where fewer interior points are in range). The only supported boundary conditions are Extrapolate{N}; when an advection operator is constructed with no boundary conditions, Extrapolate{0} is added to its bcs, and a boundary whose name has no entry in bcs also falls back to Extrapolate{0}.
• The velocity field's out-of-range face indices are clamped to the domain.

An advection operator whose interior stencil is linear in the advected argument (see Operators.has_linear_stencil) is rewritten as an operator-matrix multiply when it is broadcasted (see MatrixFields.operator_matrix), with the ghost-point extrapolations folded into its matrix's boundary rows; every other advection operator is evaluated pointwise, through the callable interface described below.

By default, it is assumed that the operator is only a function of the velocity at the current face. If the operator is a function of the velocity at neighboring faces, then the operator should define

Operators.advection_velocity_width(::SomeAdvectionOperator) = Val(:neighboring) and the operator will be evaluated with the velocity at neighboring faces as well. The default is Val(:current).

The advected field is the broadcast argument following the velocity. An operator that is a function of the stencils of multiple center-valued quantities (e.g. FCTZalesak) should take a single center-valued field whose elements are tuples of those quantities (e.g. op.(v, tuple.(x, y))); each of the 4 stencil values passed to the operator is then such a tuple.

Subtypes of this abstract type that are evaluated pointwise should be callable, with a method of the form: (::SomeAdvectionOperator)(v, x⁻⁻, x⁻, x⁺, x⁺⁺, extra_params...) or (::SomeAdvectionOperator)(v⁻, v, v⁺, x⁻⁻, x⁻, x⁺, x⁺⁺, extra_params...) if the operator is a function of the velocity at neighboring faces. All velocity arguments are supplied as the contravariant3 component of the face-valued velocity field, and extra_params are any broadcast arguments beyond the velocity and advected field (e.g. dt), evaluated at the current face and passed through as is. In particular, a vector-valued extra parameter is not converted: an operator that needs one in contravariant form (e.g. a velocity used only to determine the upwind direction, as in TVDLimitedFluxC2F) should require its callers to supply it as contravariant data. Subtypes that are instead rewritten as operator-matrix multiplies define their matrix rows in MatrixFields/operator_matrices.jl.

U = UpwindBiasedProductC2F(;boundaries)
U.(v, x)

Compute the product of the face-valued vector field v and a center-valued field x at cell faces by upwinding x according to the direction of v.

More precisely, it is computed based on the sign of the 3rd contravariant component, and it returns a Contravariant3Vector:

\[U(\boldsymbol{v},x)[i] = \begin{cases} v^3[i] x[i-\tfrac{1}{2}]\boldsymbol{e}_3 \textrm{, if } v^3[i] > 0 \\ v^3[i] x[i+\tfrac{1}{2}]\boldsymbol{e}_3 \textrm{, if } v^3[i] < 0 \end{cases}\]

where $\boldsymbol{e}_3$ is the 3rd covariant basis vector.

The only supported boundary condition is Extrapolate, which is also added to bcs (as Extrapolate{0}) by default when no boundary conditions are given: boundary faces are computed with the interior stencil, padding the ghost point it reaches with the boundary condition's extrapolation. The stencil only reaches a ghost point at the boundary face itself, where a single interior point is in range, so every extrapolation order reduces to the value of the closest interior point: since the padded upwind and downwind values then coincide, the boundary faces reduce to $v^3[i] x_b \boldsymbol{e}_3$, where $x_b$ is the value at the center closest to the boundary.

To prescribe the value of x used on the outside of a boundary instead, pass a SetValue: the constructor then returns a DirichletOperator that applies upwind_biased_product_c2f_dirichlet, the exact replacement for the removed SetValue stencil, fused into an enclosing broadcast with lazy boundary rows.

U = Upwind3rdOrderBiasedProductC2F(;boundaries)
U.(v, x)

Compute the product of a face-valued vector field v and a center-valued field x at cell faces by upwinding x, to third-order of accuracy, according to v

\[U(v,x)[i] = \begin{cases} v[i] \left(-2 x[i-\tfrac{3}{2}] + 10 x[i-\tfrac{1}{2}] + 4 x[i+\tfrac{1}{2}] \right) / 12 \textrm{, if } v[i] > 0 \\ v[i] \left(4 x[i-\tfrac{1}{2}] + 10 x[i+\tfrac{1}{2}] -2 x[i+\tfrac{3}{2}] \right) / 12 \textrm{, if } v[i] < 0 \end{cases}\]

This stencil is based on [2], eq. 4(a).

The only supported boundary condition is Extrapolate: boundary faces are computed with the interior stencil, padding each ghost point it reaches with the condition's extrapolation from the in-range interior points (the extrapolation order is reduced at the boundary face itself, where only 2 interior points are in range). When no boundary conditions are given, Extrapolate{0} is added to bcs by default. The extrapolations are taken along the third coordinate line; on a terrain-following grid that is not the wall-normal direction (see the note on AdvectionOperator). The flux through the boundary itself is not set by this padding: it is imposed by the enclosing operator, e.g. a DivergenceF2C operator with a SetValue boundary.

U = FCTBorisBook()
U.(v, x)

Correct the flux using the flux-corrected transport formulation by Boris and Book [3].

Input arguments:

• a face-valued vector field v
• a center-valued field x

\[Ac(v,x)[i] = s[i] \max \left\{0, \min \left[ |v[i] |, s[i] \left( x[i+\tfrac{3}{2}] - x[i+\tfrac{1}{2}] \right) , s[i] \left( x[i-\tfrac{1}{2}] - x[i-\tfrac{3}{2}] \right) \right] \right\},\]

where $s[i] = +1$ if $v[i] \geq 0$ and $s[i] = -1$ if $v [i] \leq 0$, and $Ac$ represents the resulting corrected antidiffusive flux. This formulation is based on [3], as reported in [4] section 5.4.1.

As for all AdvectionOperators, boundary faces are computed with the interior stencil, padding ghost points with the Extrapolate boundary condition's extrapolation (Extrapolate{0} is added to bcs by default when no boundary conditions are given). With the default, the padded values make the one-sided difference of x on the boundary side vanish at the two faces nearest each boundary, and that difference bounds the corrected antidiffusive flux, so the flux is zero there.

U = FCTZalesak()
U.(A, tuple.(Φ, Φᵗᵈ))

Correct the flux using the flux-corrected transport formulation by Zalesak [5].

Input arguments:

• a face-valued vector field A
• a center-valued field whose elements are 2-tuples of Φ and Φᵗᵈ

\[Φ_j^{n+1} = Φ_j^{td} - (C_{j+\frac{1}{2}}A_{j+\frac{1}{2}} - C_{j-\frac{1}{2}}A_{j-\frac{1}{2}})\]

This stencil is based on [5], as reported in [4] section 5.4.2, where $C$ denotes the corrected antidiffusive flux.

As for all AdvectionOperators, boundary faces are computed with the interior stencil, padding ghost points with the Extrapolate boundary condition's extrapolation (Extrapolate{0} is added to bcs by default when no boundary conditions are given); the extrapolation of a tuple-valued field applies to each of Φ and Φᵗᵈ. No value is imposed at the faces nearest each boundary: the corrected antidiffusive flux there is whatever the padded stencil gives.

LVL = LinVanLeerC2F(; constraint)
LVL.(v, x, dt)

Compute the product of the face-valued vector field v and a center-valued field x at cell faces using a slope-limited reconstruction of x, following the van Leer class of limiters as noted in [6]. Four limiter constraint options are provided:

• AlgebraicMean: Algebraic mean, this guarantees neither positivity nor monotonicity (eq 2, avg)
• PositiveDefinite: Positive-definite with implicit diffusion based on local stencil extrema (eq 3b, 3c, 5a, 5b, posd)
• MonotoneHarmonic: Monotonicity preserving harmonic mean, this implies a strong monotonicity constraint (eq 4, mono4)
• MonotoneLocalExtrema: Monotonicity preserving, with extrema bounded by the edge cells in the stencil (eq 5, mono5)

The diffusion implied by these methods is proportional to the local upwind CFL number. The mismatch Δ𝜙 = 0 returns the first-order upwind method. Special cases (discussed in Lin et al (1994)) include setting the 𝜙min = 0 or 𝜙max = saturation mixing ratio for water vapor are not considered here in favour of the generalized local extrema in equation (5a, 5b).

As for all AdvectionOperators, boundary faces are computed with the interior stencil, padding ghost points with the Extrapolate boundary condition's extrapolation (Extrapolate{0} is added to bcs by default when no boundary conditions are given).

TVDLimitedFluxC2F{BCS, M} <: AdvectionOperator

U = TVDLimitedFluxC2F(; method)
U.(𝒜, Φ, 𝓊)

𝒜, following the notation of Durran (Numerical Methods for Fluid Dynamics, 2ⁿᵈ ed.) is the antidiffusive flux given by

𝒜 = ℱʰ - ℱˡ where h and l superscripts represent the high and lower order (monotone) fluxes respectively. The effect of the TVD limiters is then to adjust the flux

C_{j+1/2} is the multiplicative limiter which is a function of ```

the ratio of the slope of the solution across a cell interface.

 - `C=1` recovers the high order flux.
 - `C=0` recovers the low order flux.

Supported limiter types are

- RZeroLimiter (returns low order flux)
- RHalfLimiter (flux multiplier == 1/2)
- RMaxLimiter (returns high order flux)
- MinModLimiter
- KorenLimiter
- SuperbeeLimiter
- MonotonizedCentralLimiter

The face-valued velocity `𝓊` is only used to determine the upwind direction,
and must be supplied as contravariant data: either a `Contravariant3Vector`
field, or a scalar field holding the contravariant3 component (e.g.
`Geometry.contravariant3.(u, Fields.local_geometry_field(face_space))` for a
velocity field `u` in another basis).

As for all [`AdvectionOperator`](@ref)s, boundary faces are computed with the
interior stencil, padding ghost points with the [`Extrapolate`](@ref) boundary
condition's extrapolation (`Extrapolate{0}` is added to `bcs` by default when
no boundary conditions are given). No value is imposed at the faces nearest
each boundary: the limited flux there is whatever the padded stencil gives.

L = LeftBiasedC2F(;boundaries)
L.(x)

Interpolate a center-value field to a face-valued field from the left.

\[L(x)[i] = x[i-\tfrac{1}{2}]\]

Only the left boundary condition should be set. Currently supported is:

• SetValue(x₀): set the value to be x₀ on the boundary.

\[L(x)[\tfrac{1}{2}] = x_0\]

R = RightBiasedC2F(;boundaries)
R.(x)

Interpolate a center-valued field to a face-valued field from the right.

\[R(x)[i] = x[i+\tfrac{1}{2}]\]

Only the right boundary condition should be set. Currently supported is:

• SetValue(x₀): set the value to be x₀ on the boundary.

\[R(x)[n+\tfrac{1}{2}] = x_0\]

L = LeftBiasedF2C(;boundaries)
L.(x)

Interpolate a face-value field to a center-valued field from the left.

\[L(x)[i+\tfrac{1}{2}] = x[i]\]

Only the left boundary condition should be set. Currently supported is:

• SetValue(x₀): set the value to be x₀ on the boundary.

\[L(x)[1] = x_0\]

R = RightBiasedF2C(;boundaries)
R.(x)

Interpolate a face-valued field to a center-valued field from the right.

\[R(x)[i] = x[i+\tfrac{1}{2}]\]

Only the right boundary condition should be set. Currently supported is:

• SetValue(x₀): set the value to be x₀ on the boundary.

\[R(x)[n+\tfrac{1}{2}] = x_0\]

BottomBiasedC2F

Alias for LeftBiasedC2F: in the vertical direction, the left boundary is the bottom.

TopBiasedC2F

Alias for RightBiasedC2F: in the vertical direction, the right boundary is the top.

BottomBiasedF2C

Alias for LeftBiasedF2C: in the vertical direction, the left boundary is the bottom.

TopBiasedF2C

Alias for RightBiasedF2C: in the vertical direction, the right boundary is the top.

AbstractTVDSlopeLimiter

An asbtract TVD-slope limiter type. Use subtypes(AbstractTVDSlopeLimiter) to see the supported subtypes. See

TVDLimitedFluxC2F for the general formulation.

G = GradientF2C(;boundaryname=boundarycondition...)
G.(x)

Compute the gradient of a face-valued field x, returning a center-valued Covariant3 vector field, using the stencil:

\[G(x)[i]^3 = x[i+\tfrac{1}{2}] - x[i-\tfrac{1}{2}]\]

We note that the usual division factor $1 / \Delta z$ that appears in a first-order finite difference operator is accounted for in the LocalVector basis. Hence, users need to cast the output of the GradientF2C to a UVector, VVector or WVector, according to the type of domain on which the operator is defined.

The following boundary conditions are supported:

• by default, the value of x at the boundary face will be used.
• SetValue(x₀): calculate the gradient assuming the value at the boundary is x₀. For the left boundary, this becomes:

\[G(x)[1]³ = x[1+\tfrac{1}{2}] - x₀\]

• SetGradient(v₀): set the value of the gradient at the center closest to the boundary to be v₀. For the left boundary, this becomes:

\[G(x)[1] = v₀\]

As with GradientC2F, v₀ is projected onto the covariant 3 axis.

G = GradientC2F(;boundaryname=boundarycondition...)
G.(x)

Compute the gradient of a center-valued field x, returning a face-valued Covariant3 vector field, using the stencil:

\[G(x)[i]^3 = x[i+\tfrac{1}{2}] - x[i-\tfrac{1}{2}]\]

The following boundary conditions are supported:

• SetGradient(v₀): set the value of the gradient at the boundary to be v₀. For the left boundary, this becomes:

  \[G(x)[\tfrac{1}{2}] = v₀\]

To prescribe the boundary value of x instead, pass a SetValue: the constructor then returns a DirichletOperator that applies gradient_c2f_dirichlet, the exact replacement for the removed SetValue stencil, fused into an enclosing broadcast with lazy boundary rows.

D = DivergenceF2C(;boundaryname=boundarycondition...)
D.(v)

Compute the vertical contribution to the divergence of a face-valued field vector v, returning a center-valued scalar field, using the stencil

\[D(v)[i] = (Jv³[i+\tfrac{1}{2}] - Jv³[i-\tfrac{1}{2}]) / J[i]\]

where Jv³ is the Jacobian multiplied by the third contravariant component of v.

The following boundary conditions are supported:

• by default, the value of v at the boundary face will be used.
• SetValue(v₀): calculate the divergence assuming the value at the boundary is v₀. For the left boundary, this becomes:

\[D(v)[1] = (Jv³[1+\tfrac{1}{2}] - Jv³₀) / J[i]\]

• SetDivergence(v₀): set the divergence at the cell center closest to the boundary

• Extrapolate(): set the value at the center closest to the boundary to be the same as the neighbouring interior value. For the left boundary, this becomes:

\[D(v)[1] = D(v)[2]\]

D = DivergenceC2F(;boundaryname=boundarycondition...)
D.(v)

Compute the vertical contribution to the divergence of a center-valued field vector v, returning a face-valued scalar field, using the stencil

\[D(v)[i] = (Jv³[i+\tfrac{1}{2}] - Jv³[i-\tfrac{1}{2}]) / J[i]\]

where Jv³ is the Jacobian multiplied by the third contravariant component of v.

The following boundary conditions are supported:

• SetDivergence(x): set the value of the divergence at the boundary to be x.

  \[D(v)[\tfrac{1}{2}] = x\]

To prescribe the boundary value of v instead, pass a SetValue: the constructor then returns a DirichletOperator that applies divergence_c2f_dirichlet, the exact replacement for the removed SetValue stencil, fused into an enclosing broadcast with lazy boundary rows.

C = CurlC2F(;boundaryname=boundarycondition...)
C.(v)

Compute the vertical-derivative contribution to the curl of a center-valued covariant vector field v. It acts on the horizontal covariant components of v (that is it only depends on $v₁$ and $v₂$), and will return a face-valued horizontal contravariant vector field (that is $C(v)³ = 0$).

Specifically it approximates:

\[\begin{align*} C(v)^1 &= -\frac{1}{J} \frac{\partial v_2}{\partial \xi^3} \\ C(v)^2 &= \frac{1}{J} \frac{\partial v_1}{\partial \xi^3} \\ \end{align*}\]

using the stencils

\[\begin{align*} C(v)[i]^1 &= - \frac{1}{J[i]} (v₂[i+\tfrac{1}{2}] - v₂[i-\tfrac{1}{2}]) \\ C(v)[i]^2 &= \frac{1}{J[i]} (v₁[i+\tfrac{1}{2}] - v₁[i-\tfrac{1}{2}]) \end{align*}\]

where $v₁$ and $v₂$ are the 1st and 2nd covariant components of $v$, and $J$ is the Jacobian determinant.

The following boundary conditions are supported:

• SetCurl(v⁰): enforce the curl operator output at the boundary to be the contravariant vector v⁰.

To prescribe the boundary value of v instead, pass a SetValue: the constructor then returns a DirichletOperator that applies curl_c2f_dirichlet, the exact replacement for the removed SetValue stencil, fused into an enclosing broadcast with lazy boundary rows.

SetBoundaryOperator(;boundaries...)

This operator is the identity in the interior, and replaces the value at each boundary for which a condition is given. It preserves the space of its argument, so it modifies the boundary faces of a face field or the boundary center cells of a center field. A side with no condition is left untouched.

The following boundary conditions are supported:

• SetValue(val): set the value to be val on the boundary.
• SetGradient(val): set the value to be val on the boundary, projected onto the Covariant3 axis.
• SetCurl(val): set the value to be val on the boundary, projected onto the Contravariant12 axis (the axis of CurlC2F's output).
• SetDivergence(val): set the value to be val on the boundary.

The projecting conditions exist so that this operator can reapply the boundary conditions of the operator it was derived from when a broadcast is rewritten as an operator matrix multiply; see MatrixFields.modifies_output.

DirichletOperator{Op}(bcs)

The operator returned by the constructor of Op (one of GradientC2F, DivergenceC2F, CurlC2F or UpwindBiasedProductC2F) when one of the requested boundary conditions is a SetValue, which those operators no longer support directly. Applying it with . calls the corresponding Dirichlet replacement helper (gradient_c2f_dirichlet, divergence_c2f_dirichlet, curl_c2f_dirichlet or upwind_biased_product_c2f_dirichlet) with each SetValue(x₀) unwrapped to its x₀value and every other boundary condition passed through as given. The replacement is built as a lazy stencil broadcast with lazy boundary rows, so it fuses into an enclosing broadcast like a true operator application and allocates nothing.

gradient_c2f_dirichlet(x; <boundary_name> = x₀...)

The vertical gradient of the center-valued field x interpolated to faces, with the value of x prescribed to be x₀ at each named boundary face: the exact replacement for the removed GradientC2F(<boundary_name> = SetValue(x₀)).(x), built (for bottom and top boundaries) as

GradientC2F(
    bottom = SetGradient(Geometry.Covariant3Vector.(2 .* (Fields.level(x, 1) .- x₀))),
    top = SetGradient(Geometry.Covariant3Vector.(2 .* (x₀ .- Fields.level(x, nlevels)))),
).(
    x,
)

x must have a scalar eltype. Each boundary value may be a number, a Field (on the corresponding boundary level of x's space, or on a whole space, of which only the level adjacent to the boundary is read), or an unmaterialized lazy broadcast of such fields; a boundary value that is already an AbstractBoundaryCondition (e.g. a SetGradient) is instead applied as given, so a Dirichlet value on one boundary can be combined with an explicit condition on the other; and a boundary without a prescribed value is computed as by GradientC2F without a boundary condition there. The result is materialized on the face space.

divergence_c2f_dirichlet(v; <boundary_name> = v₀...)

The vertical contribution to the divergence of the center-valued vector field v interpolated to faces, with the value of v prescribed to be v₀ at each named boundary face: the exact replacement for the removed DivergenceC2F(<boundary_name> = SetValue(v₀)).(v), built by wrapping a plain DivergenceC2F in a SetBoundaryOperator that overrides each prescribed boundary face with the removed stencil's value,

\[D(v)[\tfrac{1}{2}] = (Jv³[1] - Jv³₀) \frac{2}{J[\tfrac{1}{2}]}\]

(and its mirror image at the top), where Jv³₀ is computed from v₀ and the boundary face's local geometry.

Each boundary value may be an axis tensor such as Geometry.WVector(0.0) (a number is not meaningful here), a Field of such values (on the corresponding boundary level, or on a whole space, of which only the level adjacent to the boundary is read), or an unmaterialized lazy broadcast of such fields. A boundary value that is already an AbstractBoundaryCondition (one accepted by SetBoundaryOperator, e.g. a SetValue or SetDivergence of the operator's output) is instead imposed as given on the wrapping SetBoundaryOperator, and a boundary without a prescribed value is computed as by DivergenceC2F without a boundary condition there. The result is materialized on the face space.

curl_c2f_dirichlet(u; <boundary_name> = u₀...)

The vertical-derivative contribution to the curl of the center-valued covariant vector field u interpolated to faces, with the value of u prescribed to be u₀ at each named boundary face: the exact replacement for the removed CurlC2F(<boundary_name> = SetValue(u₀)).(u), built by supplying the removed stencil's boundary rows,

\[C(u)[\tfrac{1}{2}]^1 = -(u_2[1] - u_{2,0}) \frac{2}{J[\tfrac{1}{2}]}, \quad C(u)[\tfrac{1}{2}]^2 = (u_1[1] - u_{1,0}) \frac{2}{J[\tfrac{1}{2}]}\]

(and their mirror images at the top), as SetCurl boundary conditions.

Each boundary value must have the covariant 1 and 2 components of eltype(u) (e.g. a Geometry.Covariant12Vector), and may be an axis tensor, a Field of such values (on the corresponding boundary level, or on a whole space, of which only the level adjacent to the boundary is read), or an unmaterialized lazy broadcast of such fields. A boundary value that is already an AbstractBoundaryCondition (e.g. a SetCurl) is instead applied as given, and a boundary without a prescribed value is computed as by CurlC2F without a boundary condition there. The result is materialized on the face space.

upwind_biased_product_c2f_dirichlet(v, x; <boundary_name> = x₀...)

The first-order upwind product of the face-valued vector field v and the center-valued field x, with the value of x on the outside of each named boundary prescribed to be x₀: the exact replacement for the removed UpwindBiasedProductC2F(<boundary_name> = SetValue(x₀)).(v, x), built by wrapping a plain UpwindBiasedProductC2F in a SetBoundaryOperator that overrides each prescribed boundary face with the removed stencil's value, the upwind product of v³ there with x₀ on the boundary side and the closest center value of x on the interior side.

Each boundary value may be a number, a Field (on the corresponding boundary level of x's space, or on a whole space, of which only the level adjacent to the boundary is read), or an unmaterialized lazy broadcast of such fields; a boundary value that is already an AbstractBoundaryCondition (one accepted by SetBoundaryOperator, e.g. a SetValue of the flux) is instead imposed as given on the wrapping SetBoundaryOperator; and a boundary without a prescribed value is computed as by UpwindBiasedProductC2F without a boundary condition there. The result is materialized on the face space.

AbstractBoundaryCondition

Supertype for boundary conditions of both the finite-difference and the discontinuous-Galerkin operators.

Boundary conditions for FiniteDifferenceOperators should define:

• boundary_width
• stencil_left_boundary
• stencil_right_boundary

Boundary conditions for the DG numerical-flux operators (see add_numerical_flux_boundary!) should define ghost_state.

SetCurl(val)

Set the curl at the boundary to be val.

SetValue(val)

Set the value at the boundary to be val. In the case of gradient operators, this will set the input value from which the gradient is computed.

SetGradient(val)

Set the gradient at the boundary to be val. In the case of gradient operators this will set the output value of the gradient.

SetDivergence(val)

Set the divergence at the boundary to be val.

Extrapolate{N}()
Extrapolate(N = 0)

Evaluate the same stencil as the interior, but pad each ghost point the stencil reaches with a value extrapolated (with an order-N polynomial) from the N + 1 closest interior points. Currently, only 0 <= N <= 2 is supported.

If a stencil at a face i is a function of the values at x[i-3/2], x[i-1/2], x[i+1/2], x[i+3/2], then at the face i = 3/2 the single ghost point x[0] is padded with the weighted sum of the interior points x[1], x[2], x[3], with the following weights:

| N | x[1] | x[2] | x[3] | |:- |:–– |:–– |:–– | | 0 | 1 | 0 | 0 | | 1 | 2 | -1 | 0 | | 2 | 3 | -3 | 1 |

Only the interior points that the stencil can reach are available for the extrapolation, and if a ghost point requires more interior points than are available, N is reduced until the ghost point can be extrapolated with the available interior points. For example, if N = 2 and the stencil above is evaluated at the boundary face i = 1/2, only the 2 interior points x[1], x[2] are available, so both ghost points are padded with the N = 1 extrapolation:

x[-1] = x[0] = 2 * x[1] - x[2]

Note that every ghost point of a stencil is padded with the same extrapolated value: the extrapolation continues the field along the third coordinate line with a single boundary value, rather than evaluating the extrapolating polynomial at each ghost point's own position.

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

Add the numerical flux at the internal faces of the spectral space mesh.

The numerical flux is determined by evaluating

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 should return the net flux from the "minus" side to the "plus" side.

For consistency, it should satisfy the property that

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 space mesh:

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_internal!. 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::AbstractNumericalFlux, bc::AbstractBoundaryCondition, dydt, args...)

Add numerical flux at boundaries using a typed boundary condition. Constructs the ghost state via ghost_state(bc, normal, argvals⁻) and applies the numerical flux.

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

Add symmetric face lifting terms at internal faces — 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_internal!. 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...)

WJ-normalized DG face-lifting correction field of element type T: applies add_lifting_flux_internal! with face function fn to a zero residual on the space of args[1] and divides 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 (Souza et al. 2023, JAMES, 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_internal! 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_ldg_laplacian_flux_internal!(dydt, q, G, κ, τ)

Add consistent LDG/SIPG face coupling $−\{\!\{κ G\}\!\}·n̂ + τ[[q]]$ to a WJ-weighted Laplacian residual. The method with a leading ghost_exchange consumes a shared start_dg_ghost_exchange handle started on (q, G, κ).

ldg_laplacian_tendency(q, ρ_weight, κ, τ)

WJ-normalized interior-penalty Laplacian tendency approximating $κ ∇⋅(ρ_{weight} ∇q)$ (or $κ ∇²q$ when ρ_weight === nothing): weak-form volume term plus the consistent numerical flux $−\{\!\{κ G\}\!\}·n̂ + τ [\![q]\!]$ with $G = ρ_{weight} ∇q$ (or $G = ∇q$). See LDGLaplacianFlux and ldg_penalty_parameter.

ldg_penalty_parameter(κ, space)

Interior-penalty scaling $τ = κ (2N_q − 1)^2 / h$ using the horizontal spectral-element length scale (works for extruded hybrid spaces).

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_internal!(ex, numflux, dydt, y)
Operators.add_lifting_flux_internal!(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.

AbstractNumericalFlux

Abstract type for numerical flux functions used in DG methods.

CentralNumericalFlux(fluxfn)

Evaluates the central numerical flux using fluxfn.

RusanovNumericalFlux(fluxfn, wavespeedfn)

Evaluates the Rusanov numerical flux using fluxfn with wavespeed wavespeedfn

LDGLaplacianFlux(τ)

Consistent interior-penalty flux for the LDG/SIPG Laplacian. Called through add_numerical_flux_internal! on a WJ-weighted residual of $−∇·F$ with $F = −κ G$ and $G = ∇q$ (or $ρ_{weight} ∇q$). Arguments are (q, G, κ) on each side, where G must share the face-normal basis (typically Geometry.UVVector); returns $−\{\!\{κ G\}\!\}·n̂ + τ[[q]]$ with $[[q]] = q⁻ − q⁺$.

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

Symmetric central lifting completing 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_internal! / lifting_correction.

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

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⁺, λ⁺))

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

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

Construct the exterior-side argument tuple for the given 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.

PeriodicBC <: AbstractBoundaryCondition

Periodic boundary condition (handled by topology, no ghost state needed).

ReflectingWallBC <: AbstractBoundaryCondition

Reflecting wall boundary condition (no-normal-flow). Reflects normal momentum component; preserves density and potential temperature.

column_integral_definite!(ϕ_top, ᶜ∂ϕ∂z, [ϕ_bot])

Sets ϕ_top```{}= \frac{1}{ΔA(z_{bot})}\int_{z_{bot}}^{z_{top}}\, ```ᶜ∂ϕ∂z```(z)\,ΔA(z)\,dz +{}```ϕ_bot, where $z_{bot}$ and $z_{top}$ are the values of z at the bottom and top of the domain, and where ΔA is the area differential J/Δz, with J denoting the metric Jacobian. The input ᶜ∂ϕ∂z should be a cell-center Field or AbstractBroadcasted, and the output ϕ_top should be a horizontal Field. The default value of ϕ_bot is 0.

column_integral_indefinite!(ᶠϕ, ᶜ∂ϕ∂z, [ϕ_bot])

Sets ᶠϕ```(z) = \frac{1}{ΔA(z_{bot})}\int_{z_{bot}}^z\,```ᶜ∂ϕ∂z```(z')\, ΔA(z')\,dz' +{}```ϕ_bot, where $z_{bot}$ is the value of z at the bottom of the domain, and where ΔA is the area differential J/Δz, with J denoting the metric Jacobian. The input ᶜ∂ϕ∂z should be a cell-center Field or AbstractBroadcasted, and the output ᶠϕ should be a cell-face Field. The default value of ϕ_bot is 0.

column_integral_indefinite!(∂ϕ∂z, ᶠϕ, [ϕ_bot], [rtol])

Sets ᶠϕ```(z) = \frac{1}{ΔA(z_{bot})}\int_{z_{bot}}^z\, ```∂ϕ∂z```(```ᶠϕ```(z'), z')\,ΔA(z')\,dz' +{}```ϕ_bot, where ∂ϕ∂z can be any scalar-valued two-argument function. When a shallow atmosphere approximation is used, ΔA = ΔA_{bot} at all values of z, and the output ᶠϕ satisfies ᶜgradᵥ.(ᶠϕ) ≈ ∂ϕ∂z.(ᶜint.(ᶠϕ), ᶜz) with a relative tolerance of rtol, where ᶜgradᵥ = GradientF2C() and ᶜint = InterpolateF2C(). When a deep atmosphere is used, the vertical gradient is replaced with an area-weighted gradient. The default value of ϕ_bot is 0, and the default value of rtol is 0.001.

column_reduce!(f, output, input; [init], [transform], [reverse])

Applies reduce to input along the vertical direction, storing the result in output. The input can be either a Field or an AbstractBroadcasted that performs pointwise or columnwise operations on Fields. Each reduced value is computed by iteratively applying f to the values in input, starting from the bottom of each column and moving upward, and the result of the final iteration is passed to the transform function before being stored in output. If init is specified, it is used as the initial value of the iteration; otherwise, the value at the starting boundary of each column in input is used as the initial value. By default, reduction starts at the bottom boundary and proceeds upward. When reverse = true, it starts at the top boundary and proceeds downward.

With first_level and last_level denoting the indices of the boundary levels of input, the default reduction in each column can be summarized as follows:

• If init is unspecified,

  reduced_value = input[first_level]
  for level in (first_level + 1):last_level
      reduced_value = f(reduced_value, input[level])
  end
  output[] = transform(reduced_value)

• If init is specified,

  reduced_value = init
  for level in first_level:last_level
      reduced_value = f(reduced_value, input[level])
  end
  output[] = transform(reduced_value)

column_accumulate!(f, output, input; [init], [transform], [reverse])

Applies accumulate to input along the vertical direction, storing the result in output. The input can be either a Field or an AbstractBroadcasted that performs pointwise or columnwise operations on Fields. By default, each accumulated value is computed by iteratively applying f to the values in input, starting from the bottom of each column and moving upward, and the result of each iteration is passed to the transform function before being stored in output. The init value is is optional for center-to-center, face-to-face, and face-to-center accumulation, but it is required for center-to-face accumulation. When reverse = true, accumulation starts at the top boundary and proceeds downward, with the corresponding staggered boundary offsets reversed.

With first_level and last_level denoting the indices of the boundary levels of input, the default accumulation in each column can be summarized as follows:

• For center-to-center and face-to-face accumulation with init unspecified,

  accumulated_value = input[first_level]
  output[first_level] = transform(accumulated_value)
  for level in (first_level + 1):last_level
      accumulated_value = f(accumulated_value, input[level])
      output[level] = transform(accumulated_value)
  end

• For center-to-center and face-to-face accumulation with init specified,

  accumulated_value = init
  for level in first_level:last_level
      accumulated_value = f(accumulated_value, input[level])
      output[level] = transform(accumulated_value)
  end

• For face-to-center accumulation with init unspecified,

  accumulated_value = input[first_level]
  for level in (first_level + 1):last_level
      accumulated_value = f(accumulated_value, input[level])
      output[level - half] = transform(accumulated_value)
  end

• For face-to-center accumulation with init specified,

  accumulated_value = f(init, input[first_level])
  for level in (first_level + 1):last_level
      accumulated_value = f(accumulated_value, input[level])
      output[level - half] = transform(accumulated_value)
  end

• For center-to-face accumulation,

  accumulated_value = init
  output[first_level - half] = transform(accumulated_value)
  for level in first_level:last_level
      accumulated_value = f(accumulated_value, input[level])
      output[level + half] = transform(accumulated_value)
  end

return_eltype(::Op, fields...)

Element type of the Field generated by an AbstractOperator when applied to the given fields.

return_space(::Op, spaces...)

Axes of the Field generated by an AbstractOperator when applied to fields on the given spaces.

stencil_interior_width(::Op, args...)

Defines the width of the interior stencil for the operator Op with the given arguments. Returns a tuple of 2-tuples: each 2-tuple should be the lower and upper bounds of the index offsets of the stencil for each argument in the stencil.

Example

stencil(::Op, arg1, arg2) = ((-half, 1+half), (0,0))

implies that at index i, the stencil accesses arg1 at i-half, i+half and i+1+half, and arg2 at index i.

stencil_interior(::Op, space, idx, args...)

Defines the stencil of the operator Op in the interior of the domain at idx; args are the input arguments.

boundary_width(::Op, ::BC, args...)

Defines the width of a boundary condition BC on an operator Op. This is the number of locations that are used in a modified stencil. Either this function, or left_interior_idx and right_interior_idx should be defined for a specific Op/BC combination.

stencil_left_boundary(op, bc, idx, hidx, args...)

The result of stencil operator op at horizontal index hidx and some vertical index idx near the left boundary, with boundary condition bc. For operators that cannot be evaluated without a boundary condition, a NullBoundaryCondition generates NaN values here.

Operators that are rewritten into an operator matrix multiply do not reach this method: their boundary rows come from MatrixFields instead, where a NullBoundaryCondition row is filled with NaNs, so the boundary output is NaN there as well.

stencil_right_boundary(op, bc, idx, hidx, args...)

The result of stencil operator op at horizontal index hidx and some vertical index idx near the right boundary, with boundary condition bc. For operators that cannot be evaluated without a boundary condition, a NullBoundaryCondition generates NaN values here.

Operators that are rewritten into an operator matrix multiply do not reach this method: their boundary rows come from MatrixFields instead, where a NullBoundaryCondition row is filled with NaNs, so the boundary output is NaN there as well.

left_interior_idx(space::AbstractSpace, op::FiniteDifferenceOperator, bc::AbstractBoundaryCondition, args..)

The index of the left-most interior point of the operator op with boundary bc when used with arguments args.... By default, this is

left_idx(space) + boundary_width(op, bc)

but can be overwritten for specific stencil types (e.g. if the stencil is assymetric).

right_interior_idx(space::AbstractSpace, op::FiniteDifferenceOperator, bc::AbstractBoundaryCondition, args..)

The index of the right-most interior point of the operator op with boundary bc when used with arguments args.... By default, this is

right_idx(space) - boundary_width(op, bc)

but can be overwritten for specific stencil types (e.g. if the stencil is assymetric).