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

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

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

• SetGradient(v): set the value at the boundary such that the gradient is v. At the left boundary the stencil is

\[I(x)[\tfrac{1}{2}] = x[1] - \frac{1}{2} v³\]

• 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.
• SetGradient: set the value at the boundary such that the gradient is 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)

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.

Supported boundary conditions are:

• SetValue(x₀): set the value of x to be x₀ in a hypothetical ghost cell on the other side of the boundary. On the left boundary the stencil is

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

• Extrapolate(): set the value of x to be the same as the closest interior point. On the left boundary, the stencil is

  \[U(\boldsymbol{v},x)[\tfrac{1}{2}] = U(\boldsymbol{v},x)[1 + \tfrac{1}{2}]\]

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).

Supported boundary conditions are:

• FirstOrderOneSided(x₀): uses the first-order downwind scheme to compute x on the left boundary, and the first-order upwind scheme to compute x on the right boundary.
• ThirdOrderOneSided(x₀): uses the third-order downwind reconstruction to compute x on the left boundary,

and the third-order upwind reconstruction to compute x on the right boundary.

U = FCTBorisBook(;boundaries)
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.

Supported boundary conditions are:

• FirstOrderOneSided(x₀): uses the first-order downwind reconstruction to compute x on the left boundary, and the first-order upwind reconstruction to compute x on the right boundary.

U = FCTZalesak(;boundaries)
U.(A, Φ, Φᵗᵈ)

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

Input arguments:

• a face-valued vector field A
• a center-valued field Φ
• a center-valued field Φᵗᵈ

\[Φ_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.

Supported boundary conditions are:

• FirstOrderOneSided(x₀): uses the first-order downwind reconstruction to compute x on the left boundary, and the first-order upwind reconstruction to compute x on the right boundary.

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

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

• 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:

\[G(x)[1]³ = G(x)[2]³\]

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:

• SetValue(x₀): calculate the gradient assuming the value at the boundary is x₀. For the left boundary, this becomes:

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

• 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₀\]

A = AdvectionF2F(;boundaries)
A.(v, θ)

Vertical advection operator at cell faces, for a face-valued velocity field v and face-valued variables θ, approximating $v^3 \partial_3 \theta$.

It uses the following stencil

\[A(v,θ)[i] = \frac{1}{2} (θ[i+1] - θ[i-1]) v³[i]\]

No boundary conditions are currently supported.

A = AdvectionC2C(;boundaries)
A.(v, θ)

Vertical advection operator at cell centers, for cell face velocity field v cell center variables θ, approximating $v^3 \partial_3 \theta$.

It uses the following stencil

\[A(v,θ)[i] = \frac{1}{2} \{ (θ[i+1] - θ[i]) v³[i+\tfrac{1}{2}] + (θ[i] - θ[i-1])v³[i-\tfrac{1}{2}]\}\]

Supported boundary conditions:

• SetValue(θ₀): set the value of θ at the boundary face to be θ₀. At the lower boundary, this is:

\[A(v,θ)[1] = \frac{1}{2} \{ (θ[2] - θ[1]) v³[1 + \tfrac{1}{2}] + (θ[1] - θ₀)v³[\tfrac{1}{2}]\}\]

• Extrapolate: use the closest interior point as the boundary value. At the lower boundary, this is:

\[A(v,θ)[1] = (θ[2] - θ[1]) v³[1 + \tfrac{1}{2}] \}\]

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

• 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:

• SetValue(v₀): calculate the divergence assuming the value at the boundary is v₀. For the left boundary, this becomes:

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

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

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

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:

• SetValue(v₀): calculate the curl assuming the value of $v$ at the boundary is v₀. For the left boundary, this becomes:

  \[C(v)[\tfrac{1}{2}]^1 = -\frac{2}{J[i]} (v_2[1] - (v₀)_2) C(v)[\tfrac{1}{2}]^2 = \frac{2}{J[i]} (v_1[1] - (v₀)_1)\]

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

SetBoundaryOperator(;boundaries...)

This operator only modifies the values at the boundary:

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

FirstOrderOneSided()

Use a first-order up/down-wind scheme to compute the value at the boundary.

ThirdOrderOneSided()

Use a third-order up/down-wind scheme to compute the value at the boundary.

AbstractBoundaryCondition

An abstract type for boundary conditions for FiniteDifferenceOperators.

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 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()

Set the value at the boundary to be the same as the closest interior point.

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

Sets ϕ_top${}= \int_{z_{bot}}^{z_{top}}\,$ᶜ∂ϕ∂z$(z)\,dz +{}$ϕ_bot, where $z_{bot}$ and $z_{top}$ are the values of z at the bottom and top of the domain, respectively. 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) = \int_{z_{bot}}^z\,$ᶜ∂ϕ∂z$(z')\,dz' +{}$ϕ_bot, where $z_{bot}$ is the value of z at the bottom of the domain. 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) = \int_{z_{bot}}^z\,$∂ϕ∂z$($ᶠϕ$(z'), z')\,dz' +{}$ϕ_bot, where ∂ϕ∂z can be any scalar-valued two-argument function. The output ᶠϕ satisfies ᶜgradᵥ.(ᶠϕ) ≈ ∂ϕ∂z.(ᶜint.(ᶠϕ), ᶜz), where ᶜgradᵥ = GradientF2C(), ᶜint = InterpolateF2C(), and ᶜz = Fields.coordinate_field(ᶜint.(ᶠϕ)).z, and where the approximation is accurate to a relative tolerance of rtol. The default value of ϕ_bot is 0, and the default value of rtol is 0.001.

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

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 bottom of each column in input is used as the initial value.

With first_level and last_level denoting the indices of the boundary levels of input, the 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])

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

With first_level and last_level denoting the indices of the boundary levels of input, the 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...)

Defines the element type of the result of operator Op

return_space(::Op, spaces...)

Defines the space upon which the operator Op returns given arguments on input 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, loc, 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, loc, idx, args...)

Defines the stencil of operator Op at idx near the left boundary, with boundary condition BC.

stencil_right_boundary(::Op, ::BC, loc, idx, args...)

Defines the stencil of operator Op at idx near the right boundary, with boundary condition BC.

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).