Operators: finite difference

FiniteDifferenceOperator

Supertype of the finite difference operators, which act along the vertical (column) direction. Subtypes 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

Supertype of the advection operators, e.g. UpwindBiasedProductC2F and FCTZalesak. Given a face-valued velocity field v and a center-valued field x, for each face i an 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} (Outflow(; order = N) is the physically named constructor); 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, the operator is assumed to be a function of the velocity at the current face only. An operator that is a function of the velocity at neighboring faces defines

Operators.advection_velocity_width(::SomeAdvectionOperator) = Val(:neighboring)

and is then 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 (Outflow), 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, which reproduces the SetValue boundary stencil exactly and fuses 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 [51], eq. 4(a).

The only supported boundary condition is Extrapolate (Outflow): 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 [17].

Arguments

• v: A face-valued vector field.
• x: A center-valued field.

\[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 [17], as reported in [52] section 5.4.1.

As for all AdvectionOperators, boundary faces are computed with the interior stencil, padding ghost points with the Extrapolate (Outflow) 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 [18].

Arguments

• A: A face-valued vector field, the antidiffusive flux.
• tuple.(Φ, Φᵗᵈ): A center-valued field whose elements are 2-tuples of the field Φ and its transported-diffused value Φᵗᵈ.

\[Φ_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 [18], as reported in [52] 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 (Outflow) 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 of [16]. dt is the time step, which enters the limiter through the local upwind CFL number. Four limiter constraint options are provided:

• AlgebraicMean(): Algebraic mean, which guarantees neither positivity nor monotonicity (eq. 2, avg).
• PositiveDefinite(): Positive-definite with implicit diffusion based on local stencil extrema (eqs. 3b, 3c, 5a, 5b, posd).
• MonotoneHarmonic(): Monotonicity-preserving harmonic mean, which 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 [16], such as setting 𝜙min = 0 or 𝜙max to the saturation mixing ratio for water vapor, are not considered here in favour of the generalized local extrema in eqs. (5a, 5b).

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

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

Limit the face-valued antidiffusive flux 𝒜 with a TVD slope limiter method, using the center-valued field Φ to compute the slope ratio and the face-valued velocity 𝓊 to determine the upwind direction.

Following the notation of [52], 𝒜 = ℱʰ - ℱˡ is the antidiffusive flux, where the superscripts h and l denote the high- and low-order (monotone) fluxes. The TVD limiter adjusts the flux to

\[F_{j+1/2} = F^{l}_{j+1/2} + C_{j+1/2} (F^{h}_{j+1/2} - F^{l}_{j+1/2}),\]

where $C_{j+1/2}$ is the multiplicative limiter, a function of the ratio r of the upwind slope of Φ to the slope across the cell interface. C = 1 recovers the high-order flux and C = 0 the low-order flux. The operator returns $C_{j+1/2} 𝒜_{j+1/2}$.

The supported methods are the subtypes of AbstractTVDSlopeLimiter: RZeroLimiter() (returns the low-order flux), RHalfLimiter() (flux multiplier 1/2), RMaxLimiter() (returns the high-order flux), MinModLimiter(), KorenLimiter(), SuperbeeLimiter(), and 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 AdvectionOperators, boundary faces are computed with the interior stencil, padding ghost points with the Extrapolate (Outflow) 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.

B = BottomBiasedC2F(;boundaries)
B.(x)

Interpolate a center-valued field to a face-valued field from below.

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

Only the bottom boundary condition can be set. The supported condition is:

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

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

T = TopBiasedC2F(;boundaries)
T.(x)

Interpolate a center-valued field to a face-valued field from above.

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

Only the top boundary condition can be set. The supported condition is:

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

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

B = BottomBiasedF2C(;boundaries)
B.(x)

Interpolate a face-valued field to a center-valued field from below.

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

Only the bottom boundary condition can be set. The supported condition is:

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

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

T = TopBiasedF2C(;boundaries)
T.(x)

Interpolate a face-valued field to a center-valued field from above.

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

Only the top boundary condition can be set. The supported condition is:

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

\[T(x)[n] = x_0\]

AbstractTVDSlopeLimiter

Supertype of the TVD slope limiters used by TVDLimitedFluxC2F, which documents the general formulation. Each subtype defines the multiplicative limiter C(r) of the slope ratio r. Subtypes: RZeroLimiter, RHalfLimiter, RMaxLimiter, MinModLimiter, KorenLimiter, SuperbeeLimiter, and MonotonizedCentralLimiter.

RZeroLimiter()

AbstractTVDSlopeLimiter with C(r) = 0, which returns the low-order flux.

RHalfLimiter()

AbstractTVDSlopeLimiter with C(r) = 1/2.

RMaxLimiter()

AbstractTVDSlopeLimiter with C(r) = 1, which returns the high-order flux.

MinModLimiter()

AbstractTVDSlopeLimiter with C(r) = max(0, min(1, r)).

KorenLimiter()

AbstractTVDSlopeLimiter with C(r) = max(0, min(2r, (1 + 2r) / 3, 2)).

SuperbeeLimiter()

AbstractTVDSlopeLimiter with C(r) = max(0, min(1, r), min(2, r)).

MonotonizedCentralLimiter()

AbstractTVDSlopeLimiter with C(r) = max(0, min(2r, (1 + r) / 2, 2)).

LimiterConstraint

Supertype of the constraint options of LinVanLeerC2F, which select how the slope of the reconstructed field is limited: AlgebraicMean, PositiveDefinite, MonotoneHarmonic, MonotoneLocalExtrema.

AlgebraicMean()

LimiterConstraint for LinVanLeerC2F: the slope is the algebraic mean of the two one-sided differences, scaled by 1 - |CFL|. It guarantees neither positivity nor monotonicity (eq. 2, avg, of [16]).

PositiveDefinite()

LimiterConstraint for LinVanLeerC2F: the mean slope is bounded by twice the distance from the topmost stencil value to the local minimum and maximum, which keeps the reconstruction positive with implicit diffusion (eqs. 3b, 3c, 5a, 5b, posd, of [16]).

MonotoneHarmonic()

LimiterConstraint for LinVanLeerC2F: the slope is the harmonic mean of the two one-sided differences when they have the same sign and zero otherwise, a strong monotonicity constraint (eq. 4, mono4, of [16]).

MonotoneLocalExtrema()

LimiterConstraint for LinVanLeerC2F: the mean slope is bounded so that the reconstructed values stay within the minimum and maximum of the three-cell stencil, preserving monotonicity (eq. 5, mono5, of [16]).

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}]\]

The usual division factor $1 / \Delta z$ of a first-order finite difference operator is accounted for in the LocalVector basis. Hence, users must cast the output of 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 (no boundary condition), the value of x at the boundary face is 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, which reproduces the SetValue boundary stencil exactly and fuses 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 (no boundary condition), the value of v at the boundary face is 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[1]\]

• SetDivergence(d₀): set the divergence at the cell center closest to the boundary to be d₀. For the left boundary, this becomes:

\[D(v)[1] = d₀\]

• Extrapolate(), equivalently Outflow(): 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, which reproduces the SetValue boundary stencil exactly and fuses 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 returns 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, which reproduces the SetValue boundary stencil exactly and fuses into an enclosing broadcast with lazy boundary rows.

SetBoundaryOperator(;boundaries...)

Return the argument unchanged in the interior, and replace the value at each boundary for which a condition is given. The operator 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 do not support directly. Applying it with . calls the corresponding Dirichlet helper (gradient_c2f_dirichlet, divergence_c2f_dirichlet, curl_c2f_dirichlet or upwind_biased_product_c2f_dirichlet) with each SetValue(x₀) unwrapped to its value x₀ and every other boundary condition passed through as given. The helper's result is a lazy stencil broadcast with lazy boundary rows, so it fuses into an enclosing broadcast like a true operator application and allocates nothing.

Fields

• bcs: NamedTuple of boundary values and conditions, keyed by boundary name.

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

Return the vertical gradient of the center-valued field x at faces, with the value of x prescribed to be x₀ at each named boundary face: the Dirichlet form of GradientC2F, equivalent to GradientC2F(<boundary_name> = SetValue(x₀)).(x) and 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₀...)

Return the vertical contribution to the divergence of the center-valued vector field v at faces, with the value of v prescribed to be v₀ at each named boundary face: the Dirichlet form of DivergenceC2F, equivalent to DivergenceC2F(<boundary_name> = SetValue(v₀)).(v) and built by wrapping a plain DivergenceC2F in a SetBoundaryOperator that overrides each prescribed boundary face with the Dirichlet 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₀...)

Return the vertical-derivative contribution to the curl of the center-valued covariant vector field u at faces, with the value of u prescribed to be u₀ at each named boundary face: the Dirichlet form of CurlC2F, equivalent to CurlC2F(<boundary_name> = SetValue(u₀)).(u) and built by supplying the Dirichlet 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₀...)

Return 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 Dirichlet form of UpwindBiasedProductC2F, equivalent to UpwindBiasedProductC2F(<boundary_name> = SetValue(x₀)).(v, x) and built by wrapping a plain UpwindBiasedProductC2F in a SetBoundaryOperator that overrides each prescribed boundary face with the Dirichlet 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 the boundary conditions of both the finite-difference and the discontinuous-Galerkin operators. The two families act on different mesh directions, so every boundary condition is either a VerticalBoundaryCondition or a HorizontalBoundaryCondition, never both.

VerticalBoundaryCondition <: AbstractBoundaryCondition

Supertype for the boundary conditions of the vertical (column) FiniteDifferenceOperators, e.g. SetValue and Extrapolate. Subtypes should define:

• boundary_width
• stencil_left_boundary
• stencil_right_boundary

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 sets 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 sets 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. 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]

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.

Outflow(; order = 0)

Construct the outflow (zero-normal-gradient family) boundary condition, which is Extrapolate{order}(), so it is accepted wherever Extrapolate is. On the finite-difference advection operators it pads the ghost points the interior stencil reaches with an order-order extrapolation from the interior; near a boundary the order is reduced when fewer than order + 1 interior points are in range, as documented for Extrapolate. Outflow() is the zero-order (constant-value) closure.

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

Set ϕ_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 must be a cell-center Field or AbstractBroadcasted, and the output ϕ_top must be a horizontal Field. The default value of ϕ_bot is 0.

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

Set ᶠϕ```(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 must be a cell-center Field or AbstractBroadcasted, and the output ᶠϕ must be a cell-face Field. The default value of ϕ_bot is 0.

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

Set ᶠϕ```(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])

Apply 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])

Apply 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 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