# Operators

*Operators* can compute spatial derivative operations.

- for performance reasons, we need to be able to "fuse" multiple operators and function applications
- Julia provides a tool for this:
**broadcasting**, with a very flexible API

Can think of operators are "pseudo-functions": can't be called directly, but act similar to functions in the context of broadcasting. They are matrix-free, in the sense that we define the *action* of the operator directly on a field, without explicitly assembling the matrix representing the discretized operator.

## Spectral element operators

### Differential Operators

`ClimaCore.Operators.Gradient`

— Type```
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

`ClimaCore.Operators.Divergence`

— Type```
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

`ClimaCore.Operators.WeakDivergence`

— Type```
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

`ClimaCore.Operators.WeakGradient`

— Type```
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.

`ClimaCore.Operators.Curl`

— Type```
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

`ClimaCore.Operators.WeakCurl`

— Type```
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\]

### Interpolation Operators

`ClimaCore.Operators.Interpolate`

— Type```
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`

.

`ClimaCore.Operators.Restrict`

— Type```
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\]

## Finite difference operators

Finite difference operators are similar with some subtle differences:

- they can change staggering (center to face, or vice versa)
- they can span multiple elements
- no DSS is required
- boundary handling may be required

We use the following convention:

- centers are indexed by integers
`1, 2, ..., n`

- faces are indexed by half integers
`half, 1+half, ..., n+half`

`ClimaCore.Operators.FiniteDifferenceOperator`

— Type`FiniteDifferenceOperator`

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

See also `AbstractBoundaryCondition`

for how to define the boundaries.

### Interpolation operators

`ClimaCore.Operators.InterpolateC2F`

— Type```
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]\]

`ClimaCore.Operators.InterpolateF2C`

— Type`InterpolateF2C()`

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

`ClimaCore.Operators.WeightedInterpolateC2F`

— Type```
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`

.

`ClimaCore.Operators.WeightedInterpolateF2C`

— Type```
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)

`ClimaCore.Operators.UpwindBiasedProductC2F`

— Type```
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}]\]

`ClimaCore.Operators.Upwind3rdOrderBiasedProductC2F`

— Type```
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.

These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a `DivergenceF2C`

operator, with a `SetValue`

boundary.

`ClimaCore.Operators.FCTBorisBook`

— Type```
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.

Similar to the `Upwind3rdOrderBiasedProductC2F`

operator, these boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a `DivergenceF2C`

operator, with a `SetValue`

boundary.

`ClimaCore.Operators.FCTZalesak`

— Type```
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.

Similar to the `Upwind3rdOrderBiasedProductC2F`

operator, these boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a `DivergenceF2C`

operator, with a `SetValue`

boundary.

`ClimaCore.Operators.LeftBiasedC2F`

— Type```
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\]

`ClimaCore.Operators.RightBiasedC2F`

— Type```
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\]

`ClimaCore.Operators.LeftBiasedF2C`

— Type```
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\]

`ClimaCore.Operators.RightBiasedF2C`

— Type```
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\]

### Derivative operators

`ClimaCore.Operators.GradientF2C`

— Type```
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]³\]

`ClimaCore.Operators.GradientC2F`

— Type```
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₀\]

`ClimaCore.Operators.AdvectionF2F`

— Type```
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.

`ClimaCore.Operators.AdvectionC2C`

— Type```
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}] \}\]

`ClimaCore.Operators.DivergenceF2C`

— Type```
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]³\]

`ClimaCore.Operators.DivergenceC2F`

— Type```
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\]

`ClimaCore.Operators.CurlC2F`

— Type```
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⁰`

.

### Other

`ClimaCore.Operators.SetBoundaryOperator`

— Type`SetBoundaryOperator(;boundaries...)`

This operator only modifies the values at the boundary:

`SetValue(val)`

: set the value to be`val`

on the boundary.

`ClimaCore.Operators.FirstOrderOneSided`

— Type`FirstOrderOneSided()`

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

`ClimaCore.Operators.ThirdOrderOneSided`

— Type`ThirdOrderOneSided()`

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

## Finite difference boundary conditions

`ClimaCore.Operators.AbstractBoundaryCondition`

— Type`AbstractBoundaryCondition`

An abstract type for boundary conditions for `FiniteDifferenceOperator`

s.

Subtypes should define:

`ClimaCore.Operators.SetCurl`

— Type`SetCurl(val)`

Set the curl at the boundary to be `val`

.

`ClimaCore.Operators.SetValue`

— Type`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.

`ClimaCore.Operators.SetGradient`

— Type`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.

`ClimaCore.Operators.SetDivergence`

— Type`SetDivergence(val)`

Set the divergence at the boundary to be `val`

.

`ClimaCore.Operators.Extrapolate`

— Type`Extrapolate()`

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

## Integrals

`ClimaCore.Operators.column_integral_definite!`

— Function`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.

`ClimaCore.Operators.column_integral_indefinite!`

— Function`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.

`ClimaCore.Operators.column_reduce!`

— Function`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 `Field`

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

`ClimaCore.Operators.column_accumulate!`

— Function`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 `Field`

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

## Internal APIs

`ClimaCore.Operators.return_eltype`

— Function`return_eltype(::Op, fields...)`

Defines the element type of the result of operator `Op`

`ClimaCore.Operators.return_space`

— Function`return_space(::Op, spaces...)`

Defines the space upon which the operator `Op`

returns given arguments on input `spaces`

.

`ClimaCore.Operators.stencil_interior_width`

— Function`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`

.

`ClimaCore.Operators.stencil_interior`

— Function`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.

`ClimaCore.Operators.boundary_width`

— Function`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.

`ClimaCore.Operators.stencil_left_boundary`

— Function`stencil_left_boundary(::Op, ::BC, loc, idx, args...)`

Defines the stencil of operator `Op`

at `idx`

near the left boundary, with boundary condition `BC`

.

`ClimaCore.Operators.stencil_right_boundary`

— Function`stencil_right_boundary(::Op, ::BC, loc, idx, args...)`

Defines the stencil of operator `Op`

at `idx`

near the right boundary, with boundary condition `BC`

.

`ClimaCore.Operators.left_interior_idx`

— Function`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).

`ClimaCore.Operators.right_interior_idx`

— Function`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).