`ClimaCore.Hypsography.LinearAdaption`

— Type`LinearAdaption(surface::Field)`

Locate the levels by linear interpolation between the surface field and the top of the domain, using the method of GalChen1975.

`ClimaCore.Hypsography.SLEVEAdaption`

— Type`SLEVEAdaption(surface::Field, ηₕ::FT, s::FT)`

Locate vertical levels using an exponential function between the surface field and the top of the domain, using the method of Schar2002. This method is modified such no warping is applied above some user defined parameter 0 ≤ ηₕ < 1.0, where the lower and upper bounds represent the domain bottom and top respectively. `s`

governs the decay rate. If the decay-scale is poorly specified (i.e., `s * zₜ`

is lower than the maximum surface elevation), a warning is thrown and `s`

is adjusted such that it `szₜ > maximum(z_surface)`

.

`ClimaCore.Hypsography.diffuse_surface_elevation!`

— Method`diffuse_surface_elevation!(f::Field; κ::T, iter::Int, dt::T)`

Option for 2nd order diffusive smoothing of generated terrain. Mutate (smooth) a given elevation profile `f`

before assigning the surface elevation to the `HypsographyAdaption`

type. A spectral second-order diffusion operator is applied with forward-Euler updates to generate profiles for each new iteration. Steps to generate smoothed terrain ( represented as a ClimaCore Field) are as follows:

- Compute discrete elevation profile f
- Compute diffuse
*surface*elevation!(f, κ, iter). f is mutated. - Define
`Hypsography.LinearAdaption(f)`

- Define
`ExtrudedFiniteDifferenceSpace`

with new surface elevation.

Default diffusion parameters are appropriate for spherical arrangements. For `zmax-zsfc`

== 𝒪(10^4), κ == 𝒪(10^8), dt == 𝒪(10⁻¹).

`ClimaCore.Hypsography.physical_z_to_ref_z`

— Method`physical_z_to_ref_z(adaption::HypsographyAdaption, z_ref::ZPoint, z_surface::ZPoint, z_top::ZPoint) :: ZPoint`

Convert physical `z`

s to reference `z`

s as prescribed by the given adaption.

This function has to be the inverse of `ref_z_to_physical_z`

.

`ClimaCore.Hypsography.ref_z_to_physical_z`

— Method`ref_z_to_physical_z(adaption::HypsographyAdaption, z_ref::ZPoint, z_surface::ZPoint, z_top::ZPoint) :: ZPoint`

Convert reference `z`

s to physical `z`

s as prescribed by the given adaption.

This function has to be the inverse of `physical_z_to_ref_z`

.

`ClimaCore.Limiters.AbstractLimiter`

— Type`ClimaCore.Limiters.QuasiMonotoneLimiter`

— Type`QuasiMonotoneLimiter`

This limiter is inspired by the one presented in Guba et al GubaOpt2014. In the reference paper, it is denoted by OP1, and is outlined in eqs. (37)-(40). Quasimonotone here is meant to be monotone with respect to the spectral element nodal values. This limiter involves solving a constrained optimization problem (a weighted least square problem up to a fixed tolerance) that is completely local to each element.

As in HOMME, the implementation idea here is the following: we need to find a grid field which is closest to the initial field (in terms of weighted sum), but satisfies the min/max constraints. So, first we find values that do not satisfy constraints and bring these values to a closest constraint. This way we introduce some change in the tracer mass, which we then redistribute so that the l2 error is smallest. This redistribution might violate constraints; thus, we do a few iterations (until `abs(Δtracer_mass) <= rtol * tracer_mass`

).

`ρq`

: tracer density Field, where`q`

denotes tracer concentration per unit mass. This can be a scalar field, or a struct-valued field.`ρ`

: fluid density Field (scalar).

**Constructor**

`limiter = QuasiMonotoneLimiter(ρq::Field; rtol = eps(eltype(parent(ρq))))`

Creates a limiter instance for the field `ρq`

with relative tolerance `rtol`

.

**Usage**

Call `compute_bounds!`

on the input fields:

`compute_bounds!(limiter, ρq, ρ)`

Then call `apply_limiter!`

on the output fields:

`apply_limiter!(ρq, ρ, limiter)`

`ClimaCore.Limiters.apply_limit_slab!`

— Method`apply_limit_slab!(slab_ρq, slab_ρ, slab_WJ, slab_q_bounds, rtol)`

Apply the computed bounds of the tracer concentration (`slab_q_bounds`

) in the limiter to `slab_ρq`

, given the total mass `slab_ρ`

, metric terms `slab_WJ`

, and relative tolerance `rtol`

. Return whether the tolerance condition could be satisfied.

`ClimaCore.Limiters.apply_limiter!`

— Method`apply_limiter!(ρq, ρ, limiter::QuasiMonotoneLimiter)`

Apply the limiter on the tracer density `ρq`

, using the computed desired bounds on the concentration `q`

and density `ρ`

as an optimal weight. This iterates over each element, calling `apply_limit_slab!`

. If the limiter fails to converge for any element, a warning is issued.

`ClimaCore.Limiters.compute_bounds!`

— Method`compute_bounds!(limiter::QuasiMonotoneLimiter, ρq::Field, ρ::Field)`

Compute the desired bounds for the tracer concentration per unit mass `q`

, based on the tracer density, `ρq`

, and density, `ρ`

, fields.

This is computed by

`compute_element_bounds!`

- starts the ghost exchange (if distributed)
`compute_neighbor_bounds_local!`

- completes the ghost exchange (if distributed)
`compute_neighbor_bounds_ghost!`

(if distributed)

`ClimaCore.Limiters.compute_element_bounds!`

— Method`compute_element_bounds!(limiter::QuasiMonotoneLimiter, ρq, ρ)`

Given two fields `ρq`

and `ρ`

, computes the min and max of `q`

in each element, storing it in `limiter.q_bounds`

.

Part of `compute_bounds!`

.

`ClimaCore.Limiters.compute_neighbor_bounds_ghost!`

— Method`compute_neighbor_bounds_ghost!(limiter::QuasiMonotoneLimiter, topology)`

Update the field `limiter.q_bounds_nbr`

based on `limiter.q_bounds`

in the ghost neighbors. This should be called after the ghost exchange has completed.

Part of `compute_bounds!`

.

`ClimaCore.Limiters.compute_neighbor_bounds_local!`

— Method`compute_neighbor_bounds_local!(limiter::QuasiMonotoneLimiter, topology)`

Update the field `limiter.q_bounds_nbr`

based on `limiter.q_bounds`

in the local neighbors.

Part of `compute_bounds!`

.

`ClimaCore.Grids.AbstractGrid`

— Type`Grids.AbstractGrid`

Grids should define the following

`topology`

: the topology of the grid`mesh`

: the mesh of the grid`domain`

: the domain of the grid`ClimaComms.context`

`ClimaComms.device`

`local_geometry_data`

: the`DataLayout`

object containing the local geometry data information

`ClimaCore.Grids.CellCenter`

— Type`CellCenter()`

Cell center location

`ClimaCore.Grids.CellFace`

— Type`CellFace()`

Cell face location

`ClimaCore.Grids.ColumnGrid`

— Type```
ColumnGrid(
full_grid :: ExtrudedFiniteDifferenceGrid,
colidx :: ColumnIndex,
)
```

A view into a column of a `ExtrudedFiniteDifferenceGrid`

. This can be used as an

`ClimaCore.Grids.ColumnIndex`

— Type`ColumnIndex(ij,h)`

An index into a column of a field. This can be used as an argument to `getindex`

of a `Field`

, to return a field on that column.

**Example**

```
colidx = ColumnIndex((1,1),1)
field[colidx]
```

`ClimaCore.Grids.ExtrudedFiniteDifferenceGrid`

— Type```
ExtrudedFiniteDifferenceGrid(
horizontal_space::AbstractSpace,
vertical_space::FiniteDifferenceSpace,
hypsography::HypsographyAdaption = Flat(),
)
```

Construct an `ExtrudedFiniteDifferenceGrid`

from the horizontal and vertical spaces.

`ClimaCore.Grids.FiniteDifferenceGrid`

— Type```
FiniteDifferenceGrid(topology::Topologies.IntervalTopology)
FiniteDifferenceGrid(mesh::Meshes.IntervalMesh)
```

Construct a `FiniteDifferenceGrid`

from an `IntervalTopology`

(or an `IntervalMesh`

).

This is an object which contains all the necessary geometric information.

To avoid unnecessary duplication, we memoize the construction of the grid.

`ClimaCore.Grids.Flat`

— Type`Flat()`

No surface hypsography.

`ClimaCore.Grids.SpectralElementGrid1D`

— Type`SpectralElementGrid1D(mesh::Meshes.IntervalMesh, quadrature_style::Quadratures.QuadratureStyle)`

A one-dimensional space: within each element the space is represented as a polynomial.

`ClimaCore.Grids.SpectralElementGrid2D`

— Type`SpectralElementSpace2D <: AbstractSpace`

A two-dimensional space: within each element the space is represented as a polynomial.

`ClimaCore.Grids.SpectralElementGrid2D`

— Method`SpectralElementSpace2D(topology, quadrature_style; enable_bubble)`

Construct a `SpectralElementSpace2D`

instance given a `topology`

and `quadrature`

. The flag `enable_bubble`

enables the `bubble correction`

for more accurate element areas.

**Input arguments:**

- topology: Topology2D
- quadrature_style: QuadratureStyle
- enable_bubble: Bool

The idea behind the so-called `bubble_correction`

is that the numerical area of the domain (e.g., the sphere) is given by the sum of nodal integration weights times their corresponding Jacobians. However, this discrete sum is not exactly equal to the exact geometric area (4pi*radius^2 for the sphere). To make these equal, the "epsilon bubble" approach modifies the inner weights in each element so that geometric and numerical areas of each element match.

Let $\Delta A^e := A^e_{exact} - A^e_{approx}$, then, in the case of linear elements, we correct $W_{i,j} J^e_{i,j}$ by:

\[\widehat{W_{i,j} J^e}_{i,j} = W_{i,j} J^e_{i,j} + \Delta A^e * W_{i,j} / Nq^2 .\]

and the case of non linear elements, by

\[\widehat{W_{i,j} J^e}_{i,j} = W_{i,j} J^e_{i,j} \left( 1 + \tilde{A}^e \right) ,\]

where $\tilde{A}^e$ is the approximated area given by the sum of the interior nodal integration weights.

Note: This is accurate only for cubed-spheres of the `Meshes.EquiangularCubedSphere`

and `Meshes.EquidistantCubedSphere`

type, not for `Meshes.ConformalCubedSphere`

.

`ClimaCore.Grids.local_geometry_data`

— Function```
Grids.local_geometry_data(
grid :: AbstractGrid,
staggering :: Union{Staggering, Nothing},
)
```

Get the `DataLayout`

object containing the local geometry data information of the `grid`

with staggering `staggering`

.

If the grid is not staggered, `staggering`

should be `nothing`

.

`ClimaCore.Grids.local_geometry_type`

— Function`Grids.local_geometry_type(::Type)`

Get the `LocalGeometry`

type.

`ClimaCore.Grids.topology`

— Function`Grids.topology(grid::AbstractGrid)`

Get the topology of a grid.

`ClimaCore.RecursiveApply`

— Module`RecursiveApply`

This module contains operators to recurse over nested `Tuple`

s or `NamedTuple`

s.

To extend to another type `T`

, define `RecursiveApply.rmap(fn, args::T...)`

`ClimaCore.RecursiveApply.radd`

— Method```
radd(X, Y)
X ⊞ Y
```

Recursively add elements of `X`

and `Y`

.

`ClimaCore.RecursiveApply.rconvert`

— Method`rconvert(T, X)`

Identical to `convert(T, X)`

, but with improved type stability for nested types.

`ClimaCore.RecursiveApply.rdiv`

— Method`rdiv(X, Y)`

Recursively divide each element of `X`

by `Y`

`ClimaCore.RecursiveApply.rmap`

— Method`rmap(fn, X...)`

Recursively apply `fn`

to each element of `X`

`ClimaCore.RecursiveApply.rmaptype`

— Method```
rmaptype(fn, T)
rmaptype(fn, T1, T2)
```

Recursively apply `fn`

to each type parameter of the type `T`

, or to each type parameter of the types `T1`

and `T2`

, where `fn`

returns a type.

`ClimaCore.RecursiveApply.rmatmul1`

— Method`rmatmul1(W, S, i, j)`

Recursive matrix product along the 1st dimension of `S`

. Equivalent to:

`mapreduce(⊠, ⊞, W[i,:], S[:,j])`

`ClimaCore.RecursiveApply.rmatmul2`

— Method`rmatmul2(W, S, i, j)`

Recursive matrix product along the 2nd dimension `S`

. Equivalent to:

`mapreduce(⊠, ⊞, W[j,:], S[i, :])`

`ClimaCore.RecursiveApply.rmul`

— Method```
rmul(X, Y)
X ⊠ Y
```

Recursively scale each element of `X`

by `Y`

.

`ClimaCore.RecursiveApply.rmuladd`

— Method`rmuladd(w, X, Y)`

Recursively add elements of `w * X + Y`

.

`ClimaCore.RecursiveApply.rpromote_type`

— Method`rpromote_type(Ts...)`

Recursively apply `promote_type`

to the input types.

`ClimaCore.RecursiveApply.rsub`

— Method```
rsub(X, Y)
X ⊟ Y
```

Recursively subtract elements of `Y`

from `X`

.

`ClimaCore.RecursiveApply.rzero`

— Method`rzero(T)`

Recursively compute the zero value of type `T`

.

`ClimaCore.Utilities.Cache`

— Module`Utilities.Cache`

ClimaCore maintains an internal cache of topology and grid objects: this ensures that if the constructor with the same arguments is invoked again (e.g. by reading from a file), the cached object will be returned (also known as *memoization*). This has two main advantages:

topology and metric information can be reused, reducing memory usage.

it is easy to check if two fields live on the same grid: we can just check if the underlying grid objects are the same (

`===`

), rather than checking all the fields are equal (via`==`

).

However this means that objects in the cache will not be removed from the garbage collector, so we provide an interface to remove these.

`ClimaCore.Utilities.Cache.cached_objects`

— Method`Utilities.Cache.cached_objects()`

List all currently cached objects.

`ClimaCore.Utilities.Cache.clean_cache!`

— Method`Utilities.Cache.clean_cache!(object)`

Remove `object`

from the cache of created objects.

In most cases, this function should not need to be called, unless you are constructing many grid objects, for example when doing a sweep over grid paramaters.

`ClimaCore.Utilities.Cache.clean_cache!`

— Method`Utilities.Cache.clean_cache!()`

Remove all objects from the cache of created objects.

In most cases, this function should not need to be called, unless you are constructing many grid objects, for example when doing a sweep over grid paramaters.

`ClimaCore.Geometry.AbstractAxis`

— Type`AbstractAxis`

An axis of a `AxisTensor`

.

`ClimaCore.Geometry.AbstractGlobalGeometry`

— Type`AbstractGlobalGeometry`

Determines the conversion from local coordinates and vector bases to a Cartesian basis.

`ClimaCore.Geometry.AbstractPoint`

— Type`AbstractPoint`

Represents a point in space.

The following types are supported:

`XPoint(x)`

`YPoint(y)`

`ZPoint(z)`

`XYPoint(x, y)`

`XZPoint(x, z)`

`XYZPoint(x, y, z)`

`LatPoint(lat)`

`LongPoint(long)`

`LatLongPoint(lat, long)`

`LatLongZPoint(lat, long, z)`

`Cartesian1Point(x1)`

`Cartesian2Point(x2)`

`Cartesian3Point(x3)`

`Cartesian12Point(x1, x2)`

`Cartesian13Point(x1, x3)`

`Cartesian123Point(x1, x2, x3)`

`ClimaCore.Geometry.AxisTensor`

— Type`AxisTensor(axes, components)`

An `AxisTensor`

is a wrapper around a `StaticArray`

, where each dimension is labelled with an `AbstractAxis`

. These axes must be consistent for operations such as addition or subtraction, or be dual to each other for operations such as multiplication.

**See also**

`components`

to obtain the underlying array.

`ClimaCore.Geometry.CartesianGlobalGeometry`

— Type`CartesianGlobalGeometry()`

Specifies that the local coordinates align with the Cartesian coordinates, e.g. `XYZPoint`

aligns with `Cartesian123Point`

, and `UVWVector`

aligns with `Cartesian123Vector`

.

`ClimaCore.Geometry.DeepSphericalGlobalGeometry`

— Type`DeepSphericalGlobalGeometry(radius)`

Similar to `SphericalGlobalGeometry`

, but for extruded spheres. In this case, it uses the "deep-atmosphere" assumption that circumference increases with `z`

.

`ClimaCore.Geometry.LocalGeometry`

— Type`LocalGeometry`

The necessary local metric information defined at each node.

`ClimaCore.Geometry.ShallowSphericalGlobalGeometry`

— Type`ShallowSphericalGlobalGeometry(radius)`

Similar to `SphericalGlobalGeometry`

, but for extruded spheres. In this case, it uses the "shallow-atmosphere" assumption that circumference is the same at all `z`

.

`ClimaCore.Geometry.SphericalGlobalGeometry`

— Type`SphericalGlobalGeometry(radius)`

Specifies that the local coordinates are specified in reference to a sphere of radius `radius`

. The `x1`

axis is aligned with the zero longitude line.

The local vector basis has `u`

in the zonal direction (with east being positive), `v`

in the meridonal (north positive), and `w`

in the radial direction (outward positive). For a point located at the pole, we take the limit along the zero longitude line:

- at the north pole, this corresponds to
`u`

being aligned with the`x2`

direction,`v`

being aligned with the negative`x1`

direction, and`w`

being aligned with the`x3`

direction. - at the south pole, this corresponds to
`u`

being aligned with the`x2`

direction,`v`

being aligned with the`x1`

direction, and`w`

being aligned with the negative`x3`

direction.

`ClimaCore.Geometry.SurfaceGeometry`

— Type`SurfaceGeometry`

The necessary local metric information defined at each node on each surface.

`ClimaCore.Geometry.bilinear_interpolate`

— Method`bilinear_interpolate(coords::NTuple{4}, ξ1, ξ2)`

Bilinear interpolate between `coords`

by parameters `ξ1, ξ2`

, each in the interval `[-1,1]`

.

`coords`

should be a 4-tuple of coordinates in counter-clockwise order.

```
4-------3
^ | |
| | |
ξ2 | |
1-------2
ξ1-->
```

`ClimaCore.Geometry.bilinear_invert`

— Method`bilinear_invert(cc::NTuple{4, V}) where {V<:SVector{2}}`

Solve find `ξ1`

and `ξ2`

such that

`bilinear_interpolate(coords, ξ1, ξ2) == 0`

See also `bilinear_interpolate`

.

`ClimaCore.Geometry.blockmat`

— Method`blockmat(m11, m22[, m12])`

Construct an `Axis2Tensor`

from sub-blocks

`ClimaCore.Geometry.components`

— Method`components(a::AxisTensor)`

Returns a `StaticArray`

containing the components of `a`

in its stored basis.

`ClimaCore.Geometry.curl_result_type`

— Method`curl_result_type(Val(I), Val(L), V)`

The return type when taking the curl along dimensions `I`

of a field of eltype `V`

, defined on dimensions `L`

Required for statically infering the result type of the curl operation for `AxisVector`

subtypes. Curl is only defined for `CovariantVector`

` field input types.

Input Vector | Operator direction | Curl output vector |
---|---|---|

Covariant12Vector | (1,2) | Contravariant3Vector |

Covariant3Vector | (1,2) | Contravariant12Vector |

Covariant123Vector | (1,2) | Contravariant123Vector |

Covariant1Vector | (1,) | Contravariant1Vector |

Covariant2Vector | (1,) | Contravariant3Vector |

Covariant3Vector | (1,) | Contravariant2Vector |

Covariant12Vector | (3,) | Contravariant12Vector |

Covariant1Vector | (3,) | Contravariant2Vector |

Covariant2Vector | (3,) | Contravariant1Vector |

Covariant3Vector | (3,) | Contravariant3Vector |

`ClimaCore.Geometry.divergence_result_type`

— Method`divergence_result_type(V)`

The return type when taking the divergence of a field of `V`

.

Required for statically infering the result type of the divergence operation for `AxisVector`

subtypes.

`ClimaCore.Geometry.dual`

— Function`dual(ax::AbstractAxis)`

The dual axis to `ax`

.

```
julia> using ClimaCore.Geometry
julia> Geometry.dual(Geometry.Covariant12Axis())
ClimaCore.Geometry.ContravariantAxis{(1, 2)}()
julia> Geometry.dual(Geometry.Cartesian123Axis())
ClimaCore.Geometry.CartesianAxis{(1, 2, 3)}()
```

`ClimaCore.Geometry.euclidean_distance`

— Method`euclidean_distance(pt1::XYPoint, pt2::XYPoint)`

Compute the 2D or 3D Euclidean distance between `pt1`

and `pt2`

.

`ClimaCore.Geometry.float_type`

— Method`float_type(T)`

Return the floating point type backing `T`

: `T`

can either be an object or a type.

`ClimaCore.Geometry.gradient_result_type`

— Method`gradient_result_type(Val(I), V)`

The return type when taking the gradient along dimension `I`

of a field `V`

.

Required for statically infering the result type of the gradient operation for `AxisVector`

subtypes.

`ClimaCore.Geometry.great_circle_distance`

— Method`great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint, global_geometry::SphericalGlobalGeometry)`

Compute the great circle (spherical geodesic) distance between `pt1`

and `pt2`

.

`ClimaCore.Geometry.great_circle_distance`

— Method`great_circle_distance(pt1::LatLongZPoint, pt2::LatLongZPoint, global_geometry::SphericalGlobalGeometry)`

Compute the great circle (spherical geodesic) distance between `pt1`

and `pt2`

.

`ClimaCore.Geometry.linear_interpolate`

— Method`interpolate(coords::NTuple{2}, ξ1)`

Interpolate between `coords`

by parameters `ξ1`

in the interval `[-1,1]`

. The type of interpolation is determined by the element type of `coords`

`ClimaCore.Geometry.outer`

— Function```
outer(x, y)
x ⊗ y
```

Compute the outer product of `x`

and `y`

. Typically `x`

will be a vector, and `y`

can be either a number, vector or tuple/named tuple.

```
# vector ⊗ scalar = vector
julia> [1.0,2.0] ⊗ 2.0
2-element Vector{Float64}:
2.0
4.0
# vector ⊗ vector = matrix
julia> [1.0,2.0] ⊗ [1.0,3.0]
2×2 Matrix{Float64}:
1.0 3.0
2.0 6.0
# vector ⊗ tuple = recursion
julia> [1.0,2.0] ⊗ (1.0, (a=2.0, b=3.0))
([1.0, 2.0], (a = [2.0, 4.0], b = [3.0, 6.0]))
```

`ClimaCore.Geometry.project`

— Function`project(axis, V[, local_geometry])`

Project the first axis component of the vector or tensor `V`

to `axis`

This is equivalent to `transform`

, but will not throw an error if the conversion is not exact.

`ClimaCore.Geometry.transform`

— Function`transform(axis, V[, local_geometry])`

Transform the first axis of the vector or tensor `V`

to `axis`

. This will throw an error if the conversion is not exact.

The conversion rules are defined as:

`v::Covariant`

=>`Local`

:`∂ξ∂x' * v`

`v::Local`

=>`Contravariant`

:`∂ξ∂x * v`

`v::Contravariant`

=>`Local`

:`∂x∂ξ * v`

`v::Local`

=>`Covariant`

:`∂x∂ξ' * v`

`v::Covariant`

=>`Contravariant`

:`∂ξ∂x * (∂ξ∂x' * v) = gⁱʲ * v`

`v::Contravariant`

=>`Covariant`

:`∂x∂ξ' * ∂x∂ξ * v = gᵢⱼ * v`

**Example**

Consider the conversion from a `Covariant12Vector`

to a `Contravariant12Axis`

. Mathematically, we can write this as

```
[ v¹ ] [g¹¹ g¹² g¹³ ] [ v₁ ]
[ v² ] = [g²¹ g²² g²³ ] * [ v₂ ]
[ v³ ] [g³¹ g³² g³³ ] [ 0 ]
```

`project`

will drop v³ term no matter what the value is, i.e. it returns

```
[ v¹ ] [g¹¹ v₁ + g¹² v₂ ]
[ v² ] = [g²¹ v₁ + g²² v₂ ]
[ 0 ] [<drops this>]
```

`transform`

will drop the v³ term, but throw an error if it is non-zero (i.e. if the conversion is not exact)

```
[ v¹ ] [g¹¹ v₁ + g¹² v₂ ]
[ v² ] = [g²¹ v₁ + g²² v₂ ]
[ 0 ] [<asserts g²³ v₁ + g²³ v₂ == 0>]
```

`ClimaCore.Geometry.@pointtype`

— Macro`@pointtype name fieldname1 ...`

Define a subtype `name`

of `AbstractPoint`

with appropriate conversion functions.

`ClimaCore.MatrixFields`

— Module`MatrixFields`

This module adds support for defining and manipulating `Field`

s that represent matrices. Specifically, it adds the `BandMatrixRow`

type, which can be used to store the entries of a band matrix. A `Field`

of `BandMatrixRow`

s on a `FiniteDifferenceSpace`

can be interpreted as a band matrix by vertically concatenating the `BandMatrixRow`

s. Similarly, a `Field`

of `BandMatrixRow`

s on an `ExtrudedFiniteDifferenceSpace`

can be interpreted as a collection of band matrices, one for each column of the `Field`

. Such `Field`

s are called `ColumnwiseBandMatrixField`

s, and this module adds the following functionality for them:

- Constructors, e.g.,
`matrix_field = @. BidiagonalMatrixRow(field1, field2)`

- Linear combinations, e.g.,
`@. 3 * matrix_field1 + matrix_field2 / 3`

- Matrix-vector multiplication, e.g.,
`@. matrix_field ⋅ field`

- Matrix-matrix multiplication, e.g.,
`@. matrix_field1 ⋅ matrix_field2`

- Compatibility with
`LinearAlgebra.I`

, e.g.,`@. matrix_field = (4I,)`

or`@. matrix_field - (4I,)`

- Integration with
`RecursiveApply`

, e.g., the entries of`matrix_field`

can be`Tuple`

s or`NamedTuple`

s instead of single values, which allows`matrix_field`

to represent multiple band matrices at the same time - Integration with
`Operators`

, e.g., the`matrix_field`

that gets applied to the argument of any`FiniteDifferenceOperator`

`op`

can be obtained using the`FiniteDifferenceOperator`

`operator_matrix(op)`

- Conversions to native array types, e.g.,
`field2arrays(matrix_field)`

can convert each column of`matrix_field`

into a`BandedMatrix`

from`BandedMatrices.jl`

- Custom printing, e.g.,
`matrix_field`

gets displayed as a`BandedMatrix`

, specifically, as the`BandedMatrix`

that corresponds to its first column

This module also adds support for defining and manipulating sparse block matrices of `Field`

s. Specifically, it adds the `FieldMatrix`

type, which is a dictionary that maps pairs of `FieldName`

s to `ColumnwiseBandMatrixField`

s or multiples of `LinearAlgebra.I`

. This comes with the following functionality:

- Addition and subtraction, e.g.,
`@. field_matrix1 + field_matrix2`

- Matrix-vector multiplication, e.g.,
`@. field_matrix * field_vector`

- Matrix-matrix multiplication, e.g.,
`@. field_matrix1 * field_matrix2`

- Integration with
`RecursiveApply`

, e.g., the entries of`field_matrix`

can be specified either as matrix`Field`

s of`Tuple`

s or`NamedTuple`

s, or as separate matrix`Field`

s of single values - The ability to solve linear equations using
`FieldMatrixSolver`

, which is a generalization of`ldiv!`

that is designed to optimize solver performance

`ClimaCore.MatrixFields.AbstractLazyOperator`

— Type`AbstractLazyOperator`

Supertype for "lazy operators", i.e., operators that can be called without any arguments by users, as long as they appear in broadcast expressions that contain at least one `Field`

. If `lazy_op`

is an `AbstractLazyOperator`

, the expression `lazy_op.()`

will internally be translated to `non_lazy_op.(fields...)`

, as long as it appears in a broadcast expression with at least one `Field`

. This translation is done by the function `replace_lazy_operator(space, lazy_op)`

, which must be implemented by every subtype of `AbstractLazyOperator`

.

`ClimaCore.MatrixFields.BandMatrixRow`

— Type`BandMatrixRow{ld}(entries...)`

Stores the nonzero entries in a row of a band matrix, starting with the lowest diagonal, which has index `ld`

. Supported operations include accessing the entry on the diagonal with index `d`

by calling `row[d]`

, taking linear combinations with other band matrix rows (and with `LinearAlgebra.I`

), and checking for equality with other band matrix rows (and with `LinearAlgebra.I`

). There are several aliases for commonly used subtypes of `BandMatrixRow`

:

`DiagonalMatrixRow(entry_1)`

`BidiagonalMatrixRow(entry_1, entry_2)`

`TridiagonalMatrixRow(entry_1, entry_2, entry_3)`

`QuaddiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4)`

`PentadiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4, entry_5)`

`ClimaCore.MatrixFields.BlockArrowheadPreconditioner`

— Type`BlockArrowheadPreconditioner(names₁...; [P_alg₁], [alg₂])`

A `PreconditionerAlgorithm`

for a 2×2 block matrix:

\[A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}\]

The `FieldName`

s in `names₁`

correspond to the subscript `₁`

, while all other `FieldName`

s correspond to the subscript `₂`

. The preconditioner `P`

is set to the following matrix:

\[P = \begin{bmatrix} P_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}, \quad \text{where } P_{11} \text{ is a diagonal matrix}\]

The internal preconditioner `P₁₁`

is generated by the `PreconditionerAlgorithm`

`P_alg₁`

, which is set to a `MainDiagonalPreconditioner`

by default. The Schur complement of `P₁₁`

in `P`

, `A₂₂ - A₂₁ * inv(P₁₁) * A₁₂`

, is inverted using the `FieldMatrixSolverAlgorithm`

`alg₂`

, which is set to a `BlockDiagonalSolve`

by default.

`ClimaCore.MatrixFields.BlockArrowheadSchurComplementPreconditioner`

— Type`BlockArrowheadSchurComplementPreconditioner(; [P_alg₁], [alg₂])`

A `PreconditionerAlgorithm`

that is equivalent to a `BlockArrowheadPreconditioner`

, but only applied to the Schur complement of `A₁₁`

in `A`

, `A₂₂ - A₂₁ * inv(A₁₁) * A₁₂`

, which is represented by a `LazySchurComplement`

. Specifically, the preconditioner this generates is the Schur complement of `P₁₁`

in `P`

, `A₂₂ - A₂₁ * inv(P₁₁) * A₁₂`

, where `P₁₁`

is generated by `P_alg₁`

. Unlike the `BlockArrowheadPreconditioner`

constructor, this constructor does not require `names₁`

because the block structure of `A`

can be inferred from the `LazySchurComplement`

.

`ClimaCore.MatrixFields.BlockArrowheadSolve`

— Type`BlockArrowheadSolve(names₁...; [alg₂])`

A `FieldMatrixSolverAlgorithm`

for a 2×2 block arrowhead matrix:

\[A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}, \quad \text{where } A_{11} \text{ is a diagonal matrix}\]

The `FieldName`

s in `names₁`

correspond to the subscript `₁`

, while all other `FieldName`

s correspond to the subscript `₂`

. This algorithm has only 1 step:

- Solve
`(A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂ = b₂ - A₂₁ * inv(A₁₁) * b₁`

for`x₂`

using the algorithm`alg₂`

, which is set to a`BlockDiagonalSolve`

by default, and set`x₁`

to`inv(A₁₁) * (b₁ - A₁₂ * x₂)`

.

Since `A₁₁`

is a diagonal matrix, `inv(A₁₁)`

is easy to compute, which means that the Schur complement of `A₁₁`

in `A`

, `A₂₂ - A₂₁ * inv(A₁₁) * A₁₂`

, as well as the vectors `b₂ - A₂₁ * inv(A₁₁) * b₁`

and `inv(A₁₁) * (b₁ - A₁₂ * x₂)`

, are also easy to compute.

This algorithm is equivalent to block Gaussian elimination with all operations inlined into a single step.

`ClimaCore.MatrixFields.BlockDiagonalPreconditioner`

— Type`BlockDiagonalPreconditioner()`

A `PreconditionerAlgorithm`

that sets `P`

to the block diagonal entries of `A`

. This is also called a "block Jacobi" preconditioner.

`ClimaCore.MatrixFields.BlockDiagonalSolve`

— Type`BlockDiagonalSolve()`

A `FieldMatrixSolverAlgorithm`

for a block diagonal matrix:

\[A = \begin{bmatrix} A_{11} & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & A_{22} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & A_{33} & \cdots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \cdots & A_{NN} \end{bmatrix}\]

This algorithm solves the `N`

block equations `Aₙₙ * xₙ = bₙ`

in sequence (though we might want to parallelize it in the future).

If `Aₙₙ`

is a diagonal matrix, the equation `Aₙₙ * xₙ = bₙ`

is solved by making a single pass over the data, setting each `xₙ[i]`

to `inv(Aₙₙ[i, i]) * bₙ[i]`

.

Otherwise, the equation `Aₙₙ * xₙ = bₙ`

is solved using Gaussian elimination (without pivoting), which makes two passes over the data. This is currently only implemented for tri-diagonal and penta-diagonal matrices `Aₙₙ`

. In Gaussian elimination, `Aₙₙ`

is effectively factorized into the product `Lₙ * Dₙ * Uₙ`

, where `Dₙ`

is a diagonal matrix, and where `Lₙ`

and `Uₙ`

are unit lower and upper triangular matrices, respectively. The first pass multiplies both sides of the equation by `inv(Lₙ * Dₙ)`

, replacing `Aₙₙ`

with `Uₙ`

and `bₙ`

with `Uₙxₙ`

, which is referred to as putting `Aₙₙ`

into "reduced row echelon form". The second pass solves `Uₙ * xₙ = Uₙxₙ`

for `xₙ`

with a unit upper triangular matrix solver, which is referred to as "back substitution". These operations can become numerically unstable when `Aₙₙ`

has entries with large disparities in magnitude, but avoiding this would require swapping the rows of `Aₙₙ`

(i.e., replacing `Dₙ`

with a partial pivoting matrix).

`ClimaCore.MatrixFields.BlockLowerTriangularSolve`

— Type`BlockLowerTriangularSolve(names₁...; [alg₁], [alg₂])`

A `FieldMatrixSolverAlgorithm`

for a 2×2 block lower triangular matrix:

\[A = \begin{bmatrix} A_{11} & \mathbf{0} \\ A_{21} & A_{22} \end{bmatrix}\]

The `FieldName`

s in `names₁`

correspond to the subscript `₁`

, while all other `FieldName`

s correspond to the subscript `₂`

. This algorithm has 2 steps:

- Solve
`A₁₁ * x₁ = b₁`

for`x₁`

using the algorithm`alg₁`

, which is set to a`BlockDiagonalSolve`

by default. - Solve
`A₂₂ * x₂ = b₂ - A₂₁ * x₁`

for`x₂`

using the algorithm`alg₂`

, which is also set to a`BlockDiagonalSolve`

by default.

`ClimaCore.MatrixFields.CustomPreconditioner`

— Type`CustomPreconditioner(M, alg)`

A `PreconditionerAlgorithm`

that sets `P`

to the `FieldMatrix`

`M`

and inverts `P`

using the `FieldMatrixSolverAlgorithm`

`alg`

.

`ClimaCore.MatrixFields.FieldMatrixSolver`

— Type`FieldMatrixSolver(alg, A, b)`

Combination of a `FieldMatrixSolverAlgorithm`

`alg`

and the cache it requires to solve the equation `A * x = b`

for `x`

. The values of `A`

and `b`

that get passed to this constructor should be `similar`

to the ones that get passed to `field_matrix_solve!`

in order to ensure that the cache gets allocated correctly.

`ClimaCore.MatrixFields.FieldMatrixSolverAlgorithm`

— Type`FieldMatrixSolverAlgorithm`

Description of how to solve an equation of the form `A * x = b`

for `x`

, where `A`

is a `FieldMatrix`

and where `x`

and `b`

are both `FieldVector`

s. Different algorithms can be nested inside each other, enabling the construction of specialized linear solvers that fully utilize the sparsity pattern of `A`

.

**Interface**

Every subtype of `FieldMatrixSolverAlgorithm`

must implement methods for the following functions:

`ClimaCore.MatrixFields.FieldMatrixWithSolver`

— Type`FieldMatrixWithSolver(A, b, [alg])`

A wrapper that combines a `FieldMatrix`

`A`

with a `FieldMatrixSolver`

that can be used to solve the equation `A * x = b`

for `x`

, where `x`

and `b`

are both `FieldVector`

s. Similar to a `LinearAlgebra.Factorization`

, this wrapper can be passed to `ldiv!`

, whereas a regular `FieldMatrix`

cannot be passed to `ldiv!`

.

By default, the `FieldMatrixSolverAlgorithm`

`alg`

is set to a `BlockDiagonalSolve`

, so a custom `alg`

must be specified when `A`

is not a block diagonal matrix.

`ClimaCore.MatrixFields.FieldName`

— Type`FieldName(name_chain...)`

A singleton type that represents a chain of `getproperty`

calls, which can be used to access a property or sub-property of an object `x`

using the function `get_field(x, name)`

. The entire object `x`

can also be accessed with the empty `FieldName()`

.

`ClimaCore.MatrixFields.FieldNameDict`

— Type```
FieldNameDict(keys, entries)
FieldNameDict{T}(key_entry_pairs...)
```

An `AbstractDict`

with keys of type `T`

that are stored as a `FieldNameSet{T}`

. There are two subtypes of `FieldNameDict`

:

`FieldMatrix`

, which maps a set of`FieldMatrixKeys`

to either`ColumnwiseBandMatrixField`

s or multiples of`LinearAlgebra.I`

; this is the only user-facing subtype of`FieldNameDict`

`FieldVectorView`

, which maps a set of`FieldVectorKeys`

to`Field`

s; this subtype is automatically generated when a`FieldVector`

is used in the same operation as a`FieldMatrix`

(e.g., when both appear in the same broadcast expression, or when both are passed to a`FieldMatrixSolver`

)

A `FieldNameDict`

can also be "lazy", which means that it can store `AbstractBroadcasted`

objects that become `Field`

s when they are materialized. Many internal operations generate lazy `FieldNameDict`

s to reduce the number of calls to `materialize!`

, since each call comes with a small performance penalty.

The entry at a specific key can be extracted by calling `dict[key]`

, and the entries that correspond to all the keys in a `FieldNameSet`

can be extracted by calling `dict[set]`

. If `dict`

is a `FieldMatrix`

, the corresponding identity matrix can be computed by calling `one(dict)`

.

When broadcasting over `FieldNameDict`

s, the following operations are supported:

- Addition and subtraction
- Multiplication, where the first argument must be a
`FieldMatrix`

- Inversion, where the argument must be a diagonal
`FieldMatrix`

, i.e., one in which every entry is either a`ColumnwiseBandMatrixField`

of`DiagonalMatrixRow`

s or a multiple of`LinearAlgebra.I`

`ClimaCore.MatrixFields.FieldNameSet`

— Type`FieldNameSet{T}(values, [name_tree])`

An `AbstractSet`

that contains values of type `T`

, which serves as an analogue of a `KeySet`

for a `FieldNameDict`

. There are two subtypes of `FieldNameSet`

:

`FieldVectorKeys`

, for which`T`

is set to`FieldName`

`FieldMatrixKeys`

, for which`T`

is set to`Tuple{FieldName, FieldName}`

; each tuple of type`T`

represents a pair of row-column indices

Since `FieldName`

s are singleton types, the result of almost any `FieldNameSet`

operation can be inferred during compilation. So, with the exception of `map`

, `foreach`

, and `set_string`

, functions of `FieldNameSet`

s do not have any performance cost at runtime (as long as their arguments are inferrable).

Unlike other `AbstractSet`

s, `FieldNameSet`

has special behavior for overlapping values. For example, the `FieldName`

s `@name(a.b)`

and `@name(a.b.c)`

overlap, so any set operation needs to first decompose `@name(a.b)`

into its child values before combining it with `@name(a.b.c)`

. In order to support this (and also to support the ability to compute set complements), `FieldNameSet`

stores a `FieldNameTree`

`name_tree`

, which it uses to infer child values. If `name_tree`

is not specified, it gets set to `nothing`

by default, which causes some `FieldNameSet`

operations to become disabled. For binary operations like `union`

or `setdiff`

, only one set needs to specify a `name_tree`

; if two sets both specify a `name_tree`

, the `name_tree`

s must be identical.

`ClimaCore.MatrixFields.FieldNameTree`

— Type`FieldNameTree(x)`

Tree of `FieldName`

s that can be used to access `x`

with `get_field(x, name)`

. Check whether a `name`

is valid by calling `is_valid_name(name, tree)`

, and extract the children of `name`

by calling `child_names(name, tree)`

.

`ClimaCore.MatrixFields.LazyFieldMatrixSolverAlgorithm`

— Type`LazyFieldMatrixSolverAlgorithm`

A `FieldMatrixSolverAlgorithm`

that does not require `A`

to be a `FieldMatrix`

, i.e., a "matrix-free" algorithm. Internally, a `FieldMatrixSolverAlgorithm`

(for example, `SchurComplementReductionSolve`

) might run a `LazyFieldMatrixSolverAlgorithm`

on a "lazy" representation of a `FieldMatrix`

(like a `LazySchurComplement`

).

The only operations used by a `LazyFieldMatrixSolverAlgorithm`

that depend on `A`

are `lazy_mul`

and, when required, `lazy_preconditioner`

. These and other lazy operations are used to minimize the number of calls to `Base.materialize!`

, since each call comes with a small performance penalty.

`ClimaCore.MatrixFields.LazySchurComplement`

— Type`LazySchurComplement(A₁₁, A₁₂, A₂₁, A₂₂, [alg₁, cache₁, A₁₂_x₂, invA₁₁_A₁₂_x₂])`

An analogue of a `FieldMatrix`

that represents the Schur complement of `A₁₁`

in `A`

, `A₂₂ - A₂₁ * inv(A₁₁) * A₁₂`

. Since `inv(A₁₁)`

will generally be a dense matrix, it would not be efficient to directly compute the Schur complement. So, this object only supports the "lazy" functions `lazy_mul`

, which allows it to be multiplied by the vector `x₂`

, and `lazy_preconditioner`

, which allows it to be approximated with a `FieldMatrix`

.

The values `alg₁`

, `cache₁`

, `A₁₂_x₂`

, and `invA₁₁_A₁₂_x₂`

need to be specified in order for `lazy_mul`

to be able to compute `inv(A₁₁) * A₁₂ * x₂`

. When a `LazySchurComplement`

is not passed to `lazy_mul`

, these values can be omitted.

`ClimaCore.MatrixFields.MainDiagonalPreconditioner`

— Type`MainDiagonalPreconditioner()`

A `PreconditionerAlgorithm`

that sets `P`

to the main diagonal of `A`

. This is also called a "Jacobi" preconditioner.

`ClimaCore.MatrixFields.MultiplyColumnwiseBandMatrixField`

— Type`MultiplyColumnwiseBandMatrixField()`

An operator that multiplies a `ColumnwiseBandMatrixField`

by another `Field`

, i.e., matrix-vector or matrix-matrix multiplication. The `⋅`

symbol is an alias for `MultiplyColumnwiseBandMatrixField()`

.

What follows is a derivation of the algorithm used by this operator with single-column `Field`

s. For `Field`

s on multiple columns, the same computation is done for each column.

In this derivation, we will use $M_1$ and $M_2$ to denote two `ColumnwiseBandMatrixField`

s, and we will use $V$ to denote a regular (vector-like) `Field`

. For both $M_1$ and $M_2$, we will use the array-like index notation $M[row, col]$ to denote $M[row][col-row]$, i.e., the entry in the `BandMatrixRow`

$M[row]$ located on the diagonal with index $col - row$. We will also use `outer_indices`

$($`space`

$)$ to denote the tuple $($`left_idx`

$($`space`

$),$`right_idx`

$($`space`

$))$.

**1. Matrix-Vector Multiplication**

From the definition of matrix-vector multiplication,

\[(M_1 ⋅ V)[i] = \sum_k M_1[i, k] * V[k].\]

To establish bounds on the values of $k$, let us define the following values:

- $li_1, ri_1 ={}$
`outer_indices`

$($`column_axes`

$(M_1))$ - $ld_1, ud_1 ={}$
`outer_diagonals`

$($`eltype`

$(M_1))$

Since $M_1[i, k]$ is only well-defined if $k$ is a valid column index and $k - i$ is a valid diagonal index, we know that

\[li_1 \leq k \leq ri_1 \quad \text{and} \quad ld_1 \leq k - i \leq ud_1.\]

Combining these into a single inequality gives us

\[\text{max}(li_1, i + ld_1) \leq k \leq \text{min}(ri_1, i + ud_1).\]

So, we can rewrite the expression for $(M_1 ⋅ V)[i]$ as

\[(M_1 ⋅ V)[i] = \sum_{k\ =\ \text{max}(li_1, i + ld_1)}^{\text{min}(ri_1, i + ud_1)} M_1[i, k] * V[k].\]

If we replace the variable $k$ with $d = k - i$ and switch from array-like indexing to `Field`

indexing, we find that

\[(M_1 ⋅ V)[i] = \sum_{d\ =\ \text{max}(li_1 - i, ld_1)}^{\text{min}(ri_1 - i, ud_1)} M_1[i][d] * V[i + d].\]

**1.1 Interior vs. Boundary Indices**

Now, suppose that the row index $i$ is such that

\[li_1 - ld_1 \leq i \leq ri_1 - ud_1.\]

If this is the case, then the bounds on $d$ can be simplified to

\[\text{max}(li_1 - i, ld_1) = ld_1 \quad \text{and} \quad \text{min}(ri_1 - i, ud_1) = ud_1.\]

The expression for $(M_1 ⋅ V)[i]$ then becomes

\[(M_1 ⋅ V)[i] = \sum_{d = ld_1}^{ud_1} M_1[i][d] * V[i + d].\]

The values of $i$ in this range are considered to be in the "interior" of the operator, while those not in this range (for which we cannot make the above simplification) are considered to be on the "boundary".

**2. Matrix-Matrix Multiplication**

From the definition of matrix-matrix multiplication,

\[(M_1 ⋅ M_2)[i, j] = \sum_k M_1[i, k] * M_2[k, j].\]

To establish bounds on the values of $j$ and $k$, let us define the following values:

- $li_1, ri_1 ={}$
`outer_indices`

$($`column_axes`

$(M_1))$ - $ld_1, ud_1 ={}$
`outer_diagonals`

$($`eltype`

$(M_1))$ - $li_2, ri_2 ={}$
`outer_indices`

$($`column_axes`

$(M_2))$ - $ld_2, ud_2 ={}$
`outer_diagonals`

$($`eltype`

$(M_2))$

In addition, let $ld_{prod}$ and $ud_{prod}$ denote the outer diagonal indices of the product matrix $M_1 ⋅ M_2$. We will derive the values of $ld_{prod}$ and $ud_{prod}$ in the last section.

Since $M_1[i, k]$ is only well-defined if $k$ is a valid column index and $k - i$ is a valid diagonal index, we know that

\[li_1 \leq k \leq ri_1 \quad \text{and} \quad ld_1 \leq k - i \leq ud_1.\]

Since $M_2[k, j]$ is only well-defined if $j$ is a valid column index and $j - k$ is a valid diagonal index, we also know that

\[li_2 \leq j \leq ri_2 \quad \text{and} \quad ld_2 \leq j - k \leq ud_2.\]

Finally, $(M_1 ⋅ M_2)[i, j]$ is only well-defined if $j - i$ is a valid diagonal index, so

\[ld_{prod} \leq j - i \leq ud_{prod}.\]

These inequalities can be combined to obtain

\[\begin{gather*} \text{max}(li_2, i + ld_{prod}) \leq j \leq \text{min}(ri_2, i + ud_{prod}) \\ \text{and} \\ \text{max}(li_1, i + ld_1, j - ud_2) \leq k \leq \text{min}(ri_1, i + ud_1, j - ld_2). \end{gather*}\]

So, we can rewrite the expression for $(M_1 ⋅ M_2)[i, j]$ as

\[\begin{gather*} (M_1 ⋅ M_2)[i, j] = \sum_{ k\ =\ \text{max}(li_1, i + ld_1, j - ud_2) }^{\text{min}(ri_1, i + ud_1, j - ld_2)} M_1[i, k] * M_2[k, j], \text{ where} \\[0.5em] \text{max}(li_2, i + ld_{prod}) \leq j \leq \text{min}(ri_2, i + ud_{prod}). \end{gather*}\]

If we replace the variable $k$ with $d = k - i$, replace the variable $j$ with $d_{prod} = j - i$, and switch from array-like indexing to `Field`

indexing, we find that

\[\begin{gather*} (M_1 ⋅ M_2)[i][d_{prod}] = \sum_{ d\ =\ \text{max}(li_1 - i, ld_1, d_{prod} - ud_2) }^{\text{min}(ri_1 - i, ud_1, d_{prod} - ld_2)} M_1[i][d] * M_2[i + d][d_{prod} - d], \text{ where} \\[0.5em] \text{max}(li_2 - i, ld_{prod}) \leq d_{prod} \leq \text{min}(ri_2 - i, ud_{prod}). \end{gather*}\]

**2.1 Interior vs. Boundary Indices**

Now, suppose that the row index $i$ is such that

\[\text{max}(li_1 - ld_1, li_2 - ld_{prod}) \leq i \leq \text{min}(ri_1 - ud_1, ri_2 - ud_{prod}).\]

If this is the case, then the bounds on $d_{prod}$ can be simplified to

\[\text{max}(li_2 - i, ld_{prod}) = ld_{prod} \quad \text{and} \quad \text{min}(ri_2 - i, ud_{prod}) = ud_{prod}.\]

Similarly, the bounds on $d$ can be simplified using the fact that

\[\text{max}(li_1 - i, ld_1) = ld_1 \quad \text{and} \quad \text{min}(ri_1 - i, ud_1) = ud_1.\]

The expression for $(M_1 ⋅ M_2)[i][d_{prod}]$ then becomes

\[\begin{gather*} (M_1 ⋅ M_2)[i][d_{prod}] = \sum_{ d\ =\ \text{max}(ld_1, d_{prod} - ud_2) }^{\text{min}(ud_1, d_{prod} - ld_2)} M_1[i][d] * M_2[i + d][d_{prod} - d], \text{ where} \\[0.5em] ld_{prod} \leq d_{prod} \leq ud_{prod}. \end{gather*}\]

The values of $i$ in this range are considered to be in the "interior" of the operator, while those not in this range (for which we cannot make these simplifications) are considered to be on the "boundary".

**2.2 $ld_{prod}$ and $ud_{prod}$**

We only need to compute $(M_1 ⋅ M_2)[i][d_{prod}]$ for values of $d_{prod}$ that correspond to a nonempty sum in the interior, i.e, those for which

\[\text{max}(ld_1, d_{prod} - ud_2) \leq \text{min}(ud_1, d_{prod} - ld_2).\]

This can be broken down into the four inequalities

\[ld_1 \leq ud_1, \qquad ld_1 \leq d_{prod} - ld_2, \qquad d_{prod} - ud_2 \leq ud_1, \quad \text{and} \quad d_{prod} - ud_2 \leq d_{prod} - ld_2.\]

By definition, $ld_1 \leq ud_1$ and $ld_2 \leq ud_2$, so the first and last inequality are always true. Rearranging the remaining two inequalities tells us that

\[ld_1 + ld_2 \leq d_{prod} \leq ud_1 + ud_2.\]

In other words, the outer diagonal indices of $M_1 ⋅ M_2$ are

\[ld_{prod} = ld_1 + ld_2 \quad \text{and} \quad ud_{prod} = ud_1 + ud_2.\]

This means that we can express the bounds on the interior values of $i$ as

\[\text{max}(li_1, li_2 - ld_2) - ld_1 \leq i \leq \text{min}(ri_1, ri_2 - ud_2) - ud_1.\]

`ClimaCore.MatrixFields.PreconditionerAlgorithm`

— Type`PreconditionerAlgorithm`

Description of how to approximate a `FieldMatrix`

or something similar like a `LazySchurComplement`

with a preconditioner `P`

for which `P * x = b`

is easy to solve for `x`

. If `P`

is a diagonal matrix, then `x`

can be computed as `@. inv(P) * b`

; otherwise, the `PreconditionerAlgorithm`

must specify a `FieldMatrixSolverAlgorithm`

that can be used to solve `P * x = b`

for `x`

.

**Interface**

Every subtype of `PreconditionerAlgorithm`

must implement methods for the following functions:

`ClimaCore.MatrixFields.SchurComplementReductionSolve`

— Type`SchurComplementReductionSolve(names₁...; [alg₁], alg₂)`

A `FieldMatrixSolverAlgorithm`

for any 2×2 block matrix:

\[A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}\]

The `FieldName`

s in `names₁`

correspond to the subscript `₁`

, while all other `FieldName`

s correspond to the subscript `₂`

. This algorithm has 3 steps:

- Solve
`A₁₁ * x₁′ = b₁`

for`x₁′`

using the algorithm`alg₁`

, which is set to a`BlockDiagonalSolve`

by default. - Solve
`(A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂ = b₂ - A₂₁ * x₁′`

for`x₂`

using the algorithm`alg₂`

. - Solve
`A₁₁ * x₁ = b₁ - A₁₂ * x₂`

for`x₁`

using the algorithm`alg₁`

.

Since `A₁₁`

is not necessarily a diagonal matrix, `inv(A₁₁)`

will generally be a dense matrix, which means that the Schur complement of `A₁₁`

in `A`

, `A₂₂ - A₂₁ * inv(A₁₁) * A₁₂`

, cannot be computed efficiently. So, `alg₂`

must be set to a `LazyFieldMatrixSolverAlgorithm`

, which can evaluate the matrix-vector product `(A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂`

without actually computing the Schur complement matrix. This involves representing the Schur complement matrix by a `LazySchurComplement`

, which uses `alg₁`

to invert `A₁₁`

when computing the matrix-vector product.

This algorithm is equivalent to block Gaussian elimination, where steps 1 and 2 put `A`

into reduced row echelon form, and step 3 performs back substitution. For more information on this algorithm, see Section 5 of Numerical solution of saddle point problems.

`ClimaCore.MatrixFields.StationaryIterativeSolve`

— Type`StationaryIterativeSolve(; [kwargs...])`

A `LazyFieldMatrixSolverAlgorithm`

that solves `A * x = b`

by setting `x`

to some initial value `x[0]`

(usually just the zero vector, $\mathbf{0}$) and then iteratively updating it to

\[x[n] = x[n - 1] + \textrm{inv}(P) * (b - A * x[n - 1]).\]

The matrix `P`

is called a "left preconditioner" for `A`

. In general, this algorithm converges more quickly when `P`

is a close approximation of `A`

, although more complicated approximations often come with a performance penalty.

**Background**

Let `x'`

denote the value of `x`

for which `A * x = b`

. Replacing `b`

with `A * x'`

in the formula for `x[n]`

tells us that

\[x[n] = x' + (I - \textrm{inv}(P) * A) * (x[n - 1] - x').\]

In other words, the error on iteration `n`

, `x[n] - x'`

, can be expressed in terms of the error on the previous iteration, `x[n - 1] - x'`

, as

\[x[n] - x' = (I - \textrm{inv}(P) * A) * (x[n - 1] - x').\]

By induction, this means that the error on iteration `n`

is

\[x[n] - x' = (I - \textrm{inv}(P) * A)^n * (x[0] - x').\]

If we pick some norm $||\cdot||$, we find that the norm of the error is bounded by

\[||x[n] - x'|| ≤ ||(I - \textrm{inv}(P) * A)^n|| * ||x[0] - x'||.\]

For any matrix $M$, the spectral radius of $M$ is defined as

\[\rho(M) = \max\{|λ| : λ \text{ is an eigenvalue of } M\}.\]

The spectral radius has the property that

\[||M^n|| \sim \rho(M)^n, \quad \text{i.e.,} \quad \lim_{n \to \infty} \frac{||M^n||}{\rho(M)^n} = 1.\]

So, as the value of `n`

increases, the norm of the error becomes bounded by

\[||x[n] - x'|| \leq \rho(I - \textrm{inv}(P) * A)^n * ||x[0] - x'||.\]

This indicates that `x[n]`

will converge to `x'`

(i.e., that the norm of the error will converge to 0) when `ρ(I - inv(P) * A) < 1`

, and that the convergence rate is roughly bounded by `ρ(I - inv(P) * A)`

for large values of `n`

. More precisely, it can be shown that `x[n]`

will converge to `x'`

if and only if `ρ(I - inv(P) * A) < 1`

. In practice, though, the convergence eventually stops due to the limits of floating point precision.

Also, if we assume that `x[n] ≈ x'`

, we can use the formula for `x[n]`

to approximate the error on the previous iteration as

\[x[n - 1] - x' ≈ x[n - 1] - x[n] = \textrm{inv}(P) * (A * x[n - 1] - b).\]

**Debugging**

This algorithm has support for 2 debugging message group names, which can be passed to the environment variable `JULIA_DEBUG`

:

`error_norm`

: prints`||x[n] - x'||₂`

on every iteration, approximating the error`x[n] - x'`

as described above`spectral_radius`

: prints`ρ(I - inv(P) * A)`

, approximating this value with the`eigsolve`

function from KrylovKit.jl

Because the `eigsolve`

function is not compatible with CUDA, debugging the spectral radius is currently not possible on GPUs.

**Keyword Arguments**

There are 4 values that can be included in `kwargs...`

:

`P_alg = nothing`

: a`PreconditionerAlgorithm`

that specifies how to compute`P`

and solve`P * x = b`

for`x`

, or`nothing`

if preconditioning is not required (in which case`P`

is effectively set to`one(A)`

)`n_iters = 1`

: the number of iterations`correlated_solves = false`

: whether to set`x[0]`

to a value of`x`

that was generated during an earlier call to`field_matrix_solve!`

, instead of setting it to $\mathbf{0}$ (it is always set to $\mathbf{0}$ on the first call to`field_matrix_solve!`

)`eigsolve_kwargs = (;)`

: keyword arguments for the`eigsolve`

function that can be used to tune its accuracy and speed (only applicable when debugging the spectral radius)`debug = nothing`

: keyword argument for printing debug information to`@debug`

. By default,`debug`

is set to true if`"error_norm"`

or`"spectral_radius"`

is in`ENV["JULIA_DEBUG"]`

, which must be set by users.

`ClimaCore.MatrixFields.WeightedPreconditioner`

— Type`WeightedPreconditioner(M, unweighted_P_alg)`

A `PreconditionerAlgorithm`

that sets `P`

to `M * P′`

, where `M`

is a diagonal `FieldMatrix`

and `P′`

is the preconditioner generated by `unweighted_P_alg`

. When the entries of `M`

are larger than 1, this is called "relaxation" or "damping"; when the entries are smaller than 1, this is called "extrapolation".

`ClimaCore.MatrixFields.ApproximateBlockArrowheadIterativeSolve`

— Method`ApproximateBlockArrowheadIterativeSolve(names₁...; [P_alg₁], [alg₁], [alg₂], [kwargs...])`

Shorthand for constructing a `SchurComplementReductionSolve`

whose `alg₂`

is set to a `StationaryIterativeSolve`

with a `BlockArrowheadSchurComplementPreconditioner`

. The keyword argument `alg₁`

is passed to the constructor for `SchurComplementReductionSolve`

, the keyword arguments `P_alg₁`

and `alg₂`

are passed to the constructor for `BlockArrowheadSchurComplementPreconditioner`

, and all other keyword arguments are passed to the constructor for `StationaryIterativeSolve`

.

This algorithm is somewhat similar to a `StationaryIterativeSolve`

with a `BlockArrowheadPreconditioner`

, but it usually converges much more quickly, i.e., the spectral radius of its iteration matrix (`I - inv(P) * A`

) tends to be significantly smaller. Roughly speaking, this is because it runs the iterative solver on an equation with fewer variables (the Schur complement equation, `(A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂ = b₂′`

), which means that, on each iteration, it accumulates less error due to coupling between variables. However, even though it converges more quickly, its iterations take longer because they involve using `alg₁`

to invert `A₁₁`

. So, when only a few iterations are needed, a `StationaryIterativeSolve`

with a `BlockArrowheadPreconditioner`

might be more performant.

This algorithm is an example of a "segregated" solve, in contrast to the alternative "coupled" solve. In the context of computational fluid dynamics, this algorithm can also be viewed as a "SIMPLE" (Semi-Implicit Method for Pressure-Linked Equations) scheme. For more information, see Sections 4, 5, and 10 of Numerical solution of saddle point problems.

`ClimaCore.MatrixFields.apply_preconditioner`

— Method`apply_preconditioner(P_alg, P_cache, P, lazy_b)`

Constructs a lazy `FieldMatrix`

(or a concrete one when possible) that represents the product `@. inv(P) * b`

. Here, `lazy_b`

is a (possibly lazy) `FieldVectorView`

that represents `b`

.

`ClimaCore.MatrixFields.band_matrix_row_type`

— Method`band_matrix_row_type(ld, ud, T)`

A shorthand for getting the subtype of `BandMatrixRow`

that has entries of type `T`

on the diagonals with indices in the range `ld:ud`

.

`ClimaCore.MatrixFields.check_field_matrix_solver`

— Function`check_field_matrix_solver(alg, cache, A, b)`

Checks that the sparsity structure of `A`

is supported by the `FieldMatrixSolverAlgorithm`

`alg`

, and that `A`

is compatible with `b`

in the equation `A * x = b`

.

`ClimaCore.MatrixFields.check_preconditioner`

— Method`check_preconditioner(P_alg, P_cache, A, b)`

Checks that `P`

is compatible with `b`

in the equation `P * x = b`

, where `P`

is the preconditioner generated by the `PreconditionerAlgorithm`

`P_alg`

for `A`

. If `P_alg`

requires a `FieldMatrixSolverAlgorithm`

`alg`

to solve the equation, this also calls `check_field_matrix_solver`

on `alg`

.

`ClimaCore.MatrixFields.column_axes`

— Function`column_axes(matrix_field, [matrix_space])`

Returns the space that corresponds to the columns of `matrix_field`

, i.e., the `axes`

of the `Field`

s by which `matrix_field`

can be multiplied. The `matrix_space`

, on the other hand, is the space that corresponds to the rows of `matrix_field`

. By default, `matrix_space`

is set to `axes(matrix_field)`

.

`ClimaCore.MatrixFields.column_field2array`

— Method`column_field2array(field)`

Converts a field defined on a `FiniteDifferenceSpace`

into either a `Vector`

or a `BandedMatrix`

, depending on whether the elements of the field are single values or `BandMatrixRow`

s. This involves copying the data stored in the field. Because `BandedMatrix`

does not currently support operations with `CuArray`

s, all GPU data is copied to the CPU.

`ClimaCore.MatrixFields.column_field2array_view`

— Method`column_field2array_view(field)`

Similar to `column_field2array(field)`

, except that this version avoids copying the data stored in the field.

`ClimaCore.MatrixFields.concrete_field_vector`

— Method`concrete_field_vector(vector)`

Converts the `FieldVectorView`

`vector`

back into a `FieldVector`

.

`ClimaCore.MatrixFields.field2arrays`

— Method`field2arrays(field)`

Converts a field defined on a `FiniteDifferenceSpace`

or on an `ExtrudedFiniteDifferenceSpace`

into a collection of arrays, each of which corresponds to a column of the field. This is done by calling `column_field2array`

on each of the field's columns.

`ClimaCore.MatrixFields.field2arrays_view`

— Method`field2arrays_view(field)`

Similar to `field2arrays(field)`

, except that this version calls `column_field2array_view`

instead of `column_field2array`

.

`ClimaCore.MatrixFields.field_matrix_solve!`

— Method`field_matrix_solve!(solver, x, A, b)`

Solves the equation `A * x = b`

for `x`

using the `FieldMatrixSolver`

`solver`

.

`ClimaCore.MatrixFields.field_matrix_solver_cache`

— Function`field_matrix_solver_cache(alg, A, b)`

Allocates the cache required by the `FieldMatrixSolverAlgorithm`

`alg`

to solve the equation `A * x = b`

.

`ClimaCore.MatrixFields.field_vector_view`

— Function`field_vector_view(x, [name_tree])`

Constructs a `FieldVectorView`

that contains all of the `Field`

s in the `FieldVector`

`x`

. The default `name_tree`

is `FieldNameTree(x)`

, but this can be modified if needed.

`ClimaCore.MatrixFields.is_lazy`

— Method`is_lazy(dict)`

Checks whether the `FieldNameDict`

`dict`

contains any un-materialized `AbstractBroadcasted`

entries.

`ClimaCore.MatrixFields.lazy_main_diagonal`

— Method`lazy_main_diagonal(matrix)`

Constructs a lazy `FieldMatrix`

that contains the main diagonal of `matrix`

.

`ClimaCore.MatrixFields.lazy_mul`

— Method`lazy_mul(A, args...)`

Constructs a lazy `FieldMatrix`

that represents the product `@. *(A, args...)`

. This involves regular broadcasting when `A`

is a `FieldMatrix`

, but it has more complex behavior for other objects like the `LazySchurComplement`

.

`ClimaCore.MatrixFields.lazy_or_concrete_preconditioner`

— Method`lazy_or_concrete_preconditioner(P_alg, P_cache, A)`

A wrapper for `lazy_preconditioner`

that turns the lazy `FieldMatrix`

`P`

into a concrete `FieldMatrix`

when the `PreconditionerAlgorithm`

`P_alg`

requires a `FieldMatrixSolverAlgorithm`

to invert it.

`ClimaCore.MatrixFields.lazy_preconditioner`

— Method`lazy_preconditioner(P_alg, A)`

Constructs a lazy `FieldMatrix`

(or a concrete one when possible) that approximates `A`

according to the `PreconditionerAlgorithm`

`P_alg`

. If `P_alg`

is `nothing`

instead of a `PreconditionerAlgorithm`

, this returns `one(A)`

.

`ClimaCore.MatrixFields.matrix_shape`

— Function`matrix_shape(matrix_field, [matrix_space])`

Returns either `Square()`

, `FaceToCenter()`

, or `CenterToFace()`

, depending on whether the diagonal indices of `matrix_field`

are `Int`

s or `PlusHalf`

s and whether `matrix_space`

is on cell centers or cell faces. By default, `matrix_space`

is set to `axes(matrix_field)`

.

`ClimaCore.MatrixFields.mul_return_type`

— Method`mul_return_type(X, Y)`

Computes the return type of `mul_with_projection(x, y, lg)`

, where `x isa X`

and `y isa Y`

. This can also be used to obtain the return type of `x * y`

, although `x * y`

will throw an error when projection is necessary.

Note that this is equivalent to calling the internal function `_return_type`

: `Base._return_type(mul_with_projection, Tuple{X, Y, LG})`

, where `lg isa LG`

.

`ClimaCore.MatrixFields.mul_with_projection`

— Method`mul_with_projection(x, y, lg)`

Similar to `x * y`

, except that this version automatically projects `y`

to avoid `DimensionMismatch`

errors for `AxisTensor`

s. For example, if `x`

is a covector along the `Covariant3Axis`

(e.g., `Covariant3Vector(1)'`

), then `y`

will be projected onto the `Contravariant3Axis`

. In general, the first axis of `y`

will be projected onto the dual of the last axis of `x`

.

`ClimaCore.MatrixFields.operator_matrix`

— Method`operator_matrix(op)`

Constructs a new operator (or operator-like object) that generates the matrix applied by `op`

to its final argument. If `op_matrix = operator_matrix(op)`

, we can use the following identities:

- When
`op`

takes one argument,`@. op(arg) == @. op_matrix() ⋅ arg`

. - When
`op`

takes multiple arguments,`@. op(args..., arg) == @. op_matrix(args...) ⋅ arg`

.

When `op`

takes more than one argument, `operator_matrix(op)`

constructs a `FiniteDifferenceOperator`

that generates the operator matrix. When `op`

only takes one argument, it instead constructs an `AbstractLazyOperator`

, which is internally converted into a `FiniteDifferenceOperator`

when used in a broadcast expression. Implementing `op_matrix`

as a lazy operator allows us to add an argument to the expression `op_matrix.()`

, and we then use this argument to infer the space and element type of the operator matrix.

As an example, the `InterpolateF2C()`

operator on a space with $n$ cell centers applies an $n \times (n + 1)$ bidiagonal matrix:

\[\textrm{interp}(arg) = \begin{bmatrix} 0.5 & 0.5 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0.5 & 0.5 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0.5 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 0.5 & 0.5 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0.5 & 0.5 \end{bmatrix} ⋅ arg\]

The `GradientF2C()`

operator applies a similar matrix, but with different entries:

\[\textrm{grad}(arg) = \begin{bmatrix} -\textbf{e}^3 & \textbf{e}^3 & 0 & \cdots & 0 & 0 & 0 \\ 0 & -\textbf{e}^3 & \textbf{e}^3 & \cdots & 0 & 0 & 0 \\ 0 & 0 & -\textbf{e}^3 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -\textbf{e}^3 & \textbf{e}^3 & 0 \\ 0 & 0 & 0 & \cdots & 0 & -\textbf{e}^3 & \textbf{e}^3 \end{bmatrix} ⋅ arg\]

The unit vector $\textbf{e}^3$, which can also be thought of as the differential along the third coordinate axis ($\textrm{d}\xi^3$), is implemented as a `Geometry.Covariant3Vector(1)`

.

Not all operators have well-defined operator matrices. For example, the operator `GradientC2F(; bottom = SetGradient(grad_b), top = SetGradient(grad_t))`

applies an affine transformation:

\[\textrm{grad}(arg) = \begin{bmatrix} grad_b \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \\ grad_t \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 0 \\ -\textbf{e}^3 & \textbf{e}^3 & 0 & \cdots & 0 & 0 \\ 0 & -\textbf{e}^3 & \textbf{e}^3 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & \textbf{e}^3 & 0 \\ 0 & 0 & 0 & \cdots & -\textbf{e}^3 & \textbf{e}^3 \\ 0 & 0 & 0 & \cdots & 0 & 0 \end{bmatrix} ⋅ arg\]

However, this simplifies to a linear transformation when $grad_b$ and $grad_t$ are both 0:

\[\textrm{grad}(arg) = \begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 0 \\ -\textbf{e}^3 & \textbf{e}^3 & 0 & \cdots & 0 & 0 \\ 0 & -\textbf{e}^3 & \textbf{e}^3 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & \textbf{e}^3 & 0 \\ 0 & 0 & 0 & \cdots & -\textbf{e}^3 & \textbf{e}^3 \\ 0 & 0 & 0 & \cdots & 0 & 0 \end{bmatrix} ⋅ arg\]

In general, when `op`

has nonzero boundary conditions that make it apply an affine transformation, `operator_matrix(op)`

will print out a warning and zero out the boundary conditions before computing the operator matrix.

In addition to affine transformations, there are also some operators that apply nonlinear transformations to their arguments; that is, transformations which cannot be accurately approximated without using more terms of the form

\[\textrm{op}(\textbf{0}) + \textrm{op}'(\textbf{0}) ⋅ arg + \textrm{op}''(\textbf{0}) ⋅ arg ⋅ arg + \ldots.\]

When `op`

is such an operator, `operator_matrix(op)`

will throw an error. In the future, we may want to modify `operator_matrix(op)`

so that it will instead return $\textrm{op}'(\textbf{0})$, where $\textbf{0} ={}$`zero.(arg)`

.

`ClimaCore.MatrixFields.outer_diagonals`

— Method`outer_diagonals(::Type{<:BandMatrixRow})`

Gets the indices of the lower and upper diagonals, `ld`

and `ud`

, of the given subtype of `BandMatrixRow`

.

`ClimaCore.MatrixFields.preconditioner_cache`

— Method`preconditioner_cache(P_alg, A, b)`

Allocates the cache required to solve the equation `P * x = b`

, where `P`

is the preconditioner generated by the `PreconditionerAlgorithm`

`P_alg`

for `A`

.

`ClimaCore.MatrixFields.replace_lazy_operator`

— Method`replace_lazy_operator(space, lazy_op)`

Generates an instance of `Base.AbstractBroadcasted`

that corresponds to the expression `lazy_op.()`

, where the broadcast in which this expression appears is being evaluated on the given `space`

. Note that the staggering (`CellCenter`

or `CellFace`

) of this `space`

depends on the specifics of the broadcast and is not predetermined.

`ClimaCore.MatrixFields.rmul_return_type`

— Method`rmul_return_type(X, Y)`

Computes the return type of `rmul_with_projection(x, y, lg)`

, where `x isa X`

and `y isa Y`

. This can also be used to obtain the return type of `rmul(x, y)`

, although `rmul(x, y)`

will throw an error when projection is necessary.

Note that this is equivalent to calling the internal function `_return_type`

: `Base._return_type(rmul_with_projection, Tuple{X, Y, LG})`

, where `lg isa LG`

.

`ClimaCore.MatrixFields.rmul_with_projection`

— Method`rmul_with_projection(x, y, lg)`

Similar to `rmul(x, y)`

, except that this version calls `mul_with_projection`

instead of `*`

.

`ClimaCore.MatrixFields.run_field_matrix_solver!`

— Function`run_field_matrix_solver!(alg, cache, x, A, b)`

Sets `x`

to the value that solves the equation `A * x = b`

using the `FieldMatrixSolverAlgorithm`

`alg`

.

`ClimaCore.MatrixFields.solver_algorithm`

— Method`solver_algorithm(P_alg)`

A `FieldMatrixSolverAlgorithm`

that can be used to solve `P * x = b`

for `x`

, where `P`

is the preconditioner generated by the `PreconditionerAlgorithm`

`P_alg`

. If `P_alg`

is `nothing`

instead of a `PreconditionerAlgorithm`

, or if `P`

is a diagonal matrix (and no solver is required to invert it), this returns `nothing`

.

`ClimaCore.MatrixFields.@name`

— Macro`@name(expr)`

Shorthand for constructing a `FieldName`

. Some examples include

`name = @name()`

, in which case`get_field(x, name)`

returns`x`

`name = @name(a)`

, in which case`get_field(x, name)`

returns`x.a`

`name = @name(a.b.c)`

, in which case`get_field(x, name)`

returns`x.a.b.c`

`name = @name(a.b.c.:(1).d)`

, in which case`get_field(x, name)`

returns`x.a.b.c.:(1).d`

This macro is preferred over the `FieldName`

constructor because it checks whether `expr`

is a syntactically valid chain of `getproperty`

calls before calling the constructor.

`ClimaCore.Utilities.UnrolledFunctions`

— Module`UnrolledFunctions`

A collection of generated functions that get unrolled during compilation, which make it possible to iterate over nonuniform collections without sacrificing type-stability.

The functions exported by this module are

`unrolled_map(f, values, [values2])`

: alternative to`map`

`unrolled_any(f, values)`

: alternative to`any`

`unrolled_all(f, values)`

: alternative to`all`

`unrolled_filter(f, values)`

: alternative to`filter`

`unrolled_foreach(f, values)`

: alternative to`foreach`

`unrolled_in(value, values)`

: alternative to`in`

`unrolled_unique(values)`

: alternative to`unique`

`unrolled_flatten(values)`

: alternative to`Iterators.flatten`

`unrolled_flatmap(f, values)`

: alternative to`Iterators.flatmap`

`unrolled_product(values1, values2)`

: alternative to`Iterators.product`

`unrolled_findonly(f, values)`

: checks that only one value satisfies`f`

, and then returns that value`unrolled_split(f, values)`

: returns a tuple that contains the result of calling`unrolled_filter`

with`f`

and the result of calling it with`!f`

`unrolled_take(values, ::Val{N})`

: alternative to`Iterators.take`

, but with an`Int`

wrapped in a`Val`

as the second argument instead of a regular`Int`

; this usually compiles more quickly than`values[1:N]`

`unrolled_drop(values, ::Val{N})`

: alternative to`Iterators.drop`

, but with an`Int`

wrapped in a`Val`

as the second argument instead of a regular`Int`

; this usually compiles more quickly than`values[(end - N + 1):end]`

`ClimaCore.Spaces`

— Module`ClimaCore.Spaces.AbstractSpace`

— Type`AbstractSpace`

Should define

`grid`

`staggering`

`space`

constructor

`ClimaCore.Spaces.FiniteDifferenceSpace`

— Type```
FiniteDifferenceSpace(
grid::Grids.FiniteDifferenceGrid,
staggering::Staggering,
)
```

`ClimaCore.Spaces.PointSpace`

— Type`PointSpace <: AbstractSpace`

A zero-dimensional space.

`ClimaCore.Spaces.SpectralElementSpace1D`

— Type`SpectralElementSpace1D(grid::SpectralElementGrid1D)`

`ClimaCore.Spaces.SpectralElementSpace2D`

— Type`SpectralElementSpace2D(grid::SpectralElementGrid1D)`

`ClimaCore.Spaces.SpectralElementSpaceSlab`

— Type`SpectralElementSpaceSlab <: AbstractSpace`

A view into a `SpectralElementSpace2D`

for a single slab.

`ClimaCore.Spaces.area`

— Method`Spaces.area(space::Spaces.AbstractSpace)`

The length/area/volume of `space`

. This is computed as the sum of the quadrature weights $W_i$ multiplied by the Jacobian determinants $J_i$:

\[\sum_i W_i J_i \approx \int_\Omega \, d \Omega\]

If `space`

is distributed, this uses a `ClimaComms.allreduce`

operation.

`ClimaCore.Spaces.local_area`

— Method`Spaces.local_area(space::Spaces.AbstractSpace)`

The length/area/volume of `space`

local to the current context. See `Spaces.area`

`ClimaCore.Spaces.node_horizontal_length_scale`

— Method`Spaces.node_horizontal_length_scale(space::AbstractSpectralElementSpace)`

The approximate length scale of the distance between nodes. This is defined as the length scale of the mesh (see `Meshes.element_horizontal_length_scale`

), divided by the number of unique quadrature points along each dimension.

`ClimaCore.Spaces.unique_nodes`

— Method`unique_nodes(space::SpectralElementSpace2D)`

An iterator over the unique nodes of `space`

. Each node is represented by the first `((i,j), e)`

triple.

This function is experimental, and may change in future.

`ClimaCore.Spaces.weighted_dss!`

— Method```
function weighted_dss!(
data::Union{
DataLayouts.IFH,
DataLayouts.VIFH,
DataLayouts.IJFH,
DataLayouts.VIJFH,
},
space::Union{
AbstractSpectralElementSpace,
ExtrudedFiniteDifferenceSpace,
},
dss_buffer::Union{DSSBuffer, Nothing},
)
```

Computes weighted dss of `data`

.

It comprises of the following steps:

1). `Spaces.weighted_dss_start!`

`ClimaCore.Spaces.weighted_dss_ghost!`

— Method```
weighted_dss_ghost!(
data::Union{
DataLayouts.IFH,
DataLayouts.VIFH,
DataLayouts.IJFH,
DataLayouts.VIJFH,
},
space::Union{
AbstractSpectralElementSpace,
ExtrudedFiniteDifferenceSpace,
},
dss_buffer::Union{DSSBuffer, Nothing},
)
```

1). Finish communications.

2). Call `Spaces.load_from_recv_buffer!`

After the communication is complete, this adds data from the recv buffer to the corresponding location in `perimeter_data`

. For ghost vertices, this data is added only to the representative vertices. The values are then scattered to other local vertices corresponding to each unique ghost vertex in `dss_local_ghost`

.

3). Call `Spaces.dss_untransform!`

on all local elements. This transforms the DSS'd local vectors back to Covariant12 vectors, and copies the DSS'd data from the `perimeter_data`

to `data`

.

`ClimaCore.Spaces.weighted_dss_internal!`

— Method```
weighted_dss_internal!(
data::Union{
DataLayouts.IFH,
DataLayouts.VIFH,
DataLayouts.IJFH,
DataLayouts.VIJFH,
},
space::Union{
AbstractSpectralElementSpace,
ExtrudedFiniteDifferenceSpace,
},
dss_buffer::DSSBuffer,
)
```

1). Apply `Spaces.dss_transform!`

on interior elements. Local elements are split into interior and perimeter elements to facilitate overlapping of communication with computation.

2). Probe communication

3). `Spaces.dss_local!`

computes the weighted DSS on local vertices and faces.

`ClimaCore.Spaces.weighted_dss_start!`

— Method```
weighted_dss_start!(
data::Union{
DataLayouts.IFH,
DataLayouts.VIFH,
DataLayouts.IJFH,
DataLayouts.VIJFH,
},
space::Union{
AbstractSpectralElementSpace,
ExtrudedFiniteDifferenceSpace,
},
dss_buffer::Union{DSSBuffer, Nothing},
)
```

It comprises of the following steps:

1). Apply `Spaces.dss_transform!`

on perimeter elements. This weights and tranforms vector fields to physical basis if needed. Scalar fields are weighted. The transformed and/or weighted perimeter `data`

is stored in `perimeter_data`

.

2). Apply `Spaces.dss_local_ghost!`

This computes partial weighted DSS on ghost vertices, using only the information from `local`

vertices.

3). `Spaces.fill_send_buffer!`

Loads the send buffer from `perimeter_data`

. For unique ghost vertices, only data from the representative ghost vertices which store result of "ghost local" DSS are loaded.

4). Start DSS communication with neighboring processes

`ClimaCore.Spaces.Δz_data`

— Method`Δz_data(space::AbstractSpace)`

A DataLayout containing the `Δz`

on a given space `space`

.

`ClimaCore.Spaces.Δz_metric_component`

— Method`Δz_metric_component(::Type{<:Goemetry.AbstractPoint})`

The index of the z-component of an abstract point in an `AxisTensor`

.

`ClimaCore.Topologies.create_dss_buffer`

— Method```
create_dss_buffer(
data::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
hspace::AbstractSpectralElementSpace,
)
```

Creates a `DSSBuffer`

for the field data corresponding to `data`

`ClimaCore.Quadratures.ClosedUniform`

— Type`ClosedUniform{Nq}()`

Uniformly-spaced quadrature including boundary.

`ClimaCore.Quadratures.GL`

— Type`GL{Nq}()`

Gauss-Legendre quadrature using `Nq`

quadrature points.

`ClimaCore.Quadratures.GLL`

— Type`GLL{Nq}()`

Gauss-Legendre-Lobatto quadrature using `Nq`

quadrature points.

`ClimaCore.Quadratures.QuadratureStyle`

— Type`ClimaCore.Quadratures.Uniform`

— Type`Uniform{Nq}()`

Uniformly-spaced quadrature.

`ClimaCore.Quadratures.barycentric_weights`

— Method`barycentric_weights(x::SVector{Nq}) where {Nq}`

The barycentric weights associated with the array of point locations `x`

:

\[w_j = \frac{1}{\prod_{k \ne j} (x_i - x_j)}\]

See Berrut2004, equation 3.2.

`ClimaCore.Quadratures.degrees_of_freedom`

— Function`degrees_of_freedom(QuadratureStyle) -> Int`

Returns the degrees*of*freedom of the `QuadratureStyle`

concrete type

`ClimaCore.Quadratures.differentiation_matrix`

— Method`differentiation_matrix(FT, quadstyle::QuadratureStyle)`

The spectral differentiation matrix at the quadrature points of `quadstyle`

, using floating point types `FT`

.

`ClimaCore.Quadratures.differentiation_matrix`

— Method`differentiation_matrix(r::SVector{Nq, T}) where {Nq, T}`

The spectral differentiation matrix for the Lagrange polynomial of degree `Nq-1`

interpolating at points `r`

.

The matrix coefficients are computed using the Berrut2004, section 9.3:

\[D_{ij} = \begin{cases} \displaystyle \frac{w_j}{w_i (x_i - x_j)} &\text{ if } i \ne j \\ -\sum_{k \ne j} D_{kj} &\text{ if } i = j \end{cases}\]

where $w_j$ are the barycentric weights, see `barycentric_weights`

.

`ClimaCore.Quadratures.interpolation_matrix`

— Method`interpolation_matrix(x::SVector, r::SVector{Nq})`

The matrix which interpolates the Lagrange polynomial of degree `Nq-1`

through the points `r`

, to points `x`

. The matrix coefficients are computed using the Barycentric formula of Berrut2004, section 4:

\[I_{ij} = \begin{cases} 1 & \text{if } x_i = r_j, \\ 0 & \text{if } x_i = r_k \text{ for } k \ne j, \\ \frac{\displaystyle \frac{w_j}{x_i - r_j}}{\displaystyle \sum_k \frac{w_k}{x_i - r_k}} & \text{otherwise,} \end{cases}\]

where $w_j$ are the barycentric weights, see `barycentric_weights`

.

`ClimaCore.Quadratures.orthonormal_poly`

— Method`V = orthonormal_poly(points, quad)`

`V_{ij}`

contains the `j-1`

th Legendre polynomial evaluated at `points[i]`

. i.e. it is the mapping from the modal to the nodal representation.

`ClimaCore.Quadratures.polynomial_degree`

— Function`polynomial_degree(QuadratureStyle) -> Int`

Returns the polynomial degree of the `QuadratureStyle`

concrete type

`ClimaCore.Quadratures.quadrature_points`

— Function`points, weights = quadrature_points(::Type{FT}, quadrature_style)`

The points and weights of the quadrature rule in floating point type `FT`

.

`ClimaCore.Domains.AbstractDomain`

— Type`AbstractDomain`

A domain represents a region of space.

`ClimaCore.Domains.IntervalDomain`

— Method```
IntervalDomain(coord⁻, coord⁺; periodic=true)
IntervalDomain(coord⁻, coord⁺; boundary_names::Tuple{Symbol,Symbol})
```

Construct a `IntervalDomain`

, the closed interval is given by `coord⁻`

, `coord⁺`

coordinate arguments.

Either a `periodic`

or `boundary_names`

keyword argument is required.

`ClimaCore.Domains.RectangleDomain`

— Method```
RectangleDomain(x1::ClosedInterval, x2::ClosedInterval;
x1boundary::Tuple{Symbol,Symbol},
x2boundary::Tuple{Symbol,Symbol},
x1periodic = false,
x2periodic = false,
)
```

Construct a `RectangularDomain`

in the horizontal. If a given x1 or x2 boundary is not periodic, then `x1boundary`

or `x2boundary`

boundary name keyword arguments must be supplied.

`ClimaCore.Domains.SphereDomain`

— Type`SphereDomain(radius)`

A domain representing the surface of a sphere with radius `radius`

.

`ClimaCore.Domains.boundary_names`

— Function`boundary_names(obj::Union{AbstractDomain, AbstractMesh, AbstractTopology})`

A tuple or vector of unique boundary names of a spatial domain.

`ClimaCore.Fields.AbstractFieldStyle`

— Type`AbstractFieldStyle`

The supertype of all broadcasting-like operations on Fields.

`ClimaCore.Fields.Field`

— Type`Field(values, space)`

A set of `values`

defined at each point of a `space`

.

`ClimaCore.Fields.FieldStyle`

— Type`FieldStyle{DS <: DataStyle}`

Standard broadcasting on Fields. Delegates the actual work to `DS`

.

`ClimaCore.Fields.FieldVector`

— Type`FieldVector`

A `FieldVector`

is a wrapper around one or more `Field`

s that acts like vector of the underlying arrays.

It is similar in spirit to `ArrayPartition`

from RecursiveArrayTools.jl, but allows referring to fields by name.

**Constructors**

`FieldVector(;name1=field1, name2=field2, ...)`

Construct a `FieldVector`

, wrapping `field1, field2, ...`

using the names `name1, name2, ...`

.

`ClimaCore.Fields.ScalarWrapper`

— Type`Fields.ScalarWrapper(val) <: AbstractArray{T,0}`

This is a wrapper around scalar values that allows them to be mutated as part of a FieldVector. A call `getproperty`

on a `FieldVector`

with this component will return a scalar, instead of the boxed object.

`Base.fill!`

— Method`fill!(field::Field, value)`

Fill `field`

with `value`

.

`Base.fill`

— Method`fill(value, space::AbstractSpace)`

Create a new `Field`

on `space`

and fill it with `value`

.

`Base.maximum`

— Method`maximum([f=identity,]v::Field)`

Approximate maximum of `v`

or `f.(v)`

over the domain.

If `v`

is a distributed field, this uses a `ClimaComms.allreduce`

operation.

`Base.ones`

— Method`ones(space::AbstractSpace)`

Construct a field on `space`

that is one everywhere.

`Base.sum`

— Method`sum([f=identity,] v::Field)`

Approximate integration of `v`

or `f.(v)`

over the domain. In an `AbstractSpectralElementSpace`

, an integral over the entire space is computed by summation over the elements of the integrand multiplied by the Jacobian determinants and the quadrature weights at each node within an element. Hence, `sum`

is computed by summation of the field values multiplied by the Jacobian determinants and quadrature weights:

\[\sum_i f(v_i) W_i J_i \approx \int_\Omega f(v) \, d \Omega\]

where $v_i$ is the value at each node, and $f$ is the identity function if not specified.

If `v`

is a distributed field, this uses a `ClimaComms.allreduce`

operation.

`Base.zeros`

— Method`zeros(space::AbstractSpace)`

Construct a field on `space`

that is zero everywhere.

`ClimaCore.Fields.array2field`

— Method`array2field(array, space)`

Wraps `array`

in a `ClimaCore`

`Field`

that is defined over `space`

. Can be used to simplify the process of getting and setting values in an `RRTMGPModel`

; e.g.

```
array2field(center_temperature, center_space) .= center_temperature_field
face_flux_field .= array2field(model.face_flux, face_space)
```

The dimensions of `array`

are assumed to be `([number of vertical nodes], number of horizontal nodes)`

. Also, `array`

must represent a `Field`

of scalars, so that the struct type of the resulting `Field`

is the same as the element type of `array`

. If this restriction were removed, one would also need to pass the desired `Field`

struct type as an argument to `array2field`

, which would then need to permute the dimensions of `array`

to match the target `DataLayout`

.

`ClimaCore.Fields.backing_array`

— Method`backing_array(x)`

The `AbstractArray`

that is backs an object `x`

, allowing it to be treated as a component of a `FieldVector`

.

`ClimaCore.Fields.bycolumn`

— Method`Fields.bycolumn(fn, space)`

Call `fn(colidx)`

to every `ColumnIndex`

`colidx`

of `space`

. This can be used to apply multiple column-wise operations in a single pass, making use of multiple threads.

On GPUs this will simply evaluate `f`

once with `colidx=:`

(i.e. it doesn't perform evaluation by columns). This may change in future.

**Example**

```
∇ = GradientF2C()
div = DivergenceC2F()
bycolumn(axes(f)) do colidx
@. ∇f[colidx] = ∇(f[colidx])
@. df[colidx] = div(∇f[colidx])
end
```

`ClimaCore.Fields.coordinate_field`

— Method`coordinate_field(space::AbstractSpace)`

Construct a `Field`

of the coordinates of the space.

`ClimaCore.Fields.field2array`

— Method`field2array(field)`

Extracts a view of a `ClimaCore`

`Field`

's underlying array. Can be used to simplify the process of getting and setting values in an `RRTMGPModel`

; e.g.

```
center_temperature .= field2array(center_temperature_field)
field2array(face_flux_field) .= face_flux
```

The dimensions of the resulting array are `([number of vertical nodes], number of horizontal nodes)`

. Also, `field`

must be a `Field`

of scalars, so that the element type of the array is the same as the struct type of `field`

.

`ClimaCore.Fields.field_iterator`

— Method`field_iterator(::Union{Field, FieldVector})`

Returns an iterable field, that recursively calls `getproperty`

for all of the `propertynames`

, and returns a `Tuple`

of

- the individual scalar field and
- the tuple chain used to reach the scalar field

`ClimaCore.Fields.local_geometry_field`

— Method`local_geometry_field(space::AbstractSpace)`

Construct a `Field`

of the `LocalGeometry`

of the space.

`ClimaCore.Fields.local_sum`

— Method`Fields.local_sum(v::Field)`

Compute the approximate integral of `v`

over the domain local to the current context.

See `sum`

for the integral over the full domain.

`ClimaCore.Fields.ncolumns`

— Method```
ncolumns(::Field)
ncolumns(::Space)
```

Number of columns in a given space.

`ClimaCore.Fields.property_chains`

— Method`property_chains`

An array of "property chains", used to recursively `getproperty`

until a single scalar field is reached.

A property chain may be, for example `(:model, :submodel, :temperature)`

where `model.submodel.temperature`

is a scalar field.

`ClimaCore.Fields.set!`

— Function`set!(f::Function, field::Field, args = ())`

Apply function `f`

to populate values in field `field`

. `f`

must have a function signature with signature `f(::LocalGeometry[, args...])`

. Additional arguments may be passed to `f`

with `args`

.

`ClimaCore.Fields.unwrap`

— Method`Fields.unwrap(x::T)`

This is called when calling `getproperty`

on a `FieldVector`

property of element type `T`

.

`ClimaCore.Fields.wrap`

— Method`Fields.wrap(x)`

Construct a mutable wrapper around `x`

. This can be extended for new types (especially immutable ones).

`ClimaCore.Fields.Δz_field`

— Method```
Δz_field(field::Field)
Δz_field(space::AbstractSpace)
```

A `Field`

containing the `Δz`

values on the same space as the given field.

`ClimaCore.Spaces.weighted_dss!`

— Function`Spaces.weighted_dss!(f::Field[, ghost_buffer = Spaces.create_dss_buffer(field)])`

Apply weighted direct stiffness summation (DSS) to `f`

. This operates in-place (i.e. it modifies the `f`

). `ghost_buffer`

contains the necessary information for communication in a distributed setting, see `Spaces.create_dss_buffer`

.

This is a projection operation from the piecewise polynomial space $\mathcal{V}_0$ to the continuous space $\mathcal{V}_1 = \mathcal{V}_0 \cap \mathcal{C}_0$, defined as the field $\theta \in \mathcal{V}_1$ such that for all $\phi \in \mathcal{V}_1$

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

In matrix form, we define $\bar \theta$ to be the unique global node representation, and $Q$ to be the "scatter" operator which maps to the redundant node representation $\theta$

\[\theta = Q \bar \theta\]

Then the problem can be written as

\[(Q \bar\phi)^\top W J Q \bar\theta = (Q \bar\phi)^\top W J f\]

which reduces to

\[\theta = Q \bar\theta = Q (Q^\top W J Q)^{-1} Q^\top W J f\]

`ClimaCore.Spaces.weighted_dss!`

— Method`Spaces.weighted_dss!(field1 => ghost_buffer1, field2 => ghost_buffer2, ...)`

Call `Spaces.weighted_dss!`

on multiple fields at once, overlapping communication as much as possible.

`ClimaCore.Topologies.create_dss_buffer`

— Method`Spaces.create_dss_buffer(field::Field)`

Create a buffer for communicating neighbour information of `field`

.

`LinearAlgebra.norm`

— Function`norm(v::Field, p=2; normalize=true)`

The approximate $L^p$ norm of `v`

, where $L^p$ represents the space of measurable functions for which the p-th power of the absolute value is Lebesgue integrable, that is:

\[\| v \|_p = \left( \int_\Omega |v|^p d \Omega \right)^{1/p}\]

where $|v|$ is defined to be the absolute value if $v$ is a scalar-valued Field, or the 2-norm if it is a vector-valued Field or composite Field (see LinearAlgebra.norm). Similar to `sum`

and `mean`

, in an `AbstractSpectralElementSpace`

, this is computed by summation of the field values multiplied by the Jacobian determinants and quadrature weights. If `normalize=true`

(the default), then internally the discrete norm is divided by the sum of the Jacobian determinants and quadrature weights:

\[\left(\frac{\sum_i |v_i|^p W_i J_i}{\sum_i W_i J_i}\right)^{1/p} \approx \left(\frac{\int_\Omega |v|^p \, d \Omega}{\int_\Omega \, d \Omega}\right)^{1/p}\]

If `p=Inf`

, then the norm is the maximum of the absolute values

\[\max_i |v_i| \approx \sup_{\Omega} |v|\]

Consequently all norms should have the same units for all $p$ (being the same as calling `norm`

on a single value).

If `normalize=false`

, then the denominator term is omitted, and so the result will be the norm as described above multiplied by the length/area/volume of $\Omega$ to the power of $1/p$.

`Statistics.mean`

— Method`mean([f=identity, ]v::Field)`

The mean value of `field`

or `f.(field)`

over the domain, weighted by area. Similar to `sum`

, in an `AbstractSpectralElementSpace`

, this is computed by summation of the field values multiplied by the Jacobian determinants and quadrature weights:

\[\frac{\sum_i f(v_i) W_i J_i}{\sum_i W_i J_i} \approx \frac{\int_\Omega f(v) \, d \Omega}{\int_\Omega \, d \Omega}\]

where $v_i$ is the Field value at each node, and $f$ is the identity function if not specified.

If `v`

is a distributed field, this uses a `ClimaComms.allreduce`

operation.

`ClimaCore.Meshes.AbstractCubedSphere`

— Type`AbstractCubedSphere <: AbstractMesh2D`

This is an abstract type of cubed-sphere meshes on `SphereDomain`

s. A cubed-sphere mesh has 6 panels, laid out as follows:

```
: Panel 1 :
+-------------+-------------+
| +x1 | +x1 |
| | |
| Panel | Panel |
|+x3 5 -x3|-x2 6 +x2|
| -x2 | -x3 |
| | |
| -x1 | -x1 |
+-------------+-------------+-------------+
| -x2 | -x2 |
| | |
| Panel | Panel |
|+x1 3 -x1|+x3 4 -x3|
| +x3 | -x1 |
| | |
| +x2 | +x2 |
+-------------+-------------+-------------+
| +x3 | +x3 |
| | |
| Panel | Panel |
|-x2 1 +x2|+x1 2 -x1|
| +x1 | +x2 |
| | |
| -x3 | -x3 |
+-------------+-------------+
: Panel 6 :
```

This is the same panel ordering used by the S2 Geometry library (though we use 1-based instead of 0-based numering).

Elements are indexed by a `CartesianIndex{3}`

object, where the components are:

- horizontal element index (left to right) within each panel.
- vertical element index (bottom to top) within each panel.
- panel number

Subtypes should have the following fields:

`domain`

: a`SphereDomain`

`ne`

: number of elements across each panel

**External links**

`ClimaCore.Meshes.AbstractMesh`

— Type`AbstractMesh{dim}`

A `Mesh`

is an object which represents how we discretize a domain into elements.

It should be lightweight (i.e. exists on all MPI ranks), e.g for meshes stored in a file, it would contain the filename.

**Face and vertex numbering**

In 1D, faces and vertices are the same, and both are numbered `[1,2]`

.

In 2D, a face is a line segment between to vertices, and both are numbered `[1,2,3,4]`

, in a counter-clockwise direction.

```
v4 f3 v3
o-----------------o
| | face vertices
| | f1 => v1 v2
f4 | | f2 f2 => v2 v3
| | f3 => v3 v4
| | f4 => v4 v1
| |
o-----------------o
v1 f1 v2
```

**Interface**

A subtype of `AbstractMesh`

should define the following methods:

`domain(mesh)`

`elements(mesh)`

`is_boundary_face(mesh, elem, face)`

`boundary_face_name(mesh, elem, face)`

`opposing_face(mesh, elem, face)`

`coordinates(mesh, elem, vert)`

`containing_element`

(optional)

The following types/methods are provided by `AbstractMesh`

:

`ClimaCore.Meshes.ConformalCubedSphere`

— Type`ConformalCubedSphere <: AbstractCubedSphere`

A conformal mesh outlined in Rancic1996. Uses the element indexing convention of `AbstractCubedSphere`

.

**Constructors**

`ConformalCubedSphere(domain::Domains.SphereDomain, ne::Integer)`

Constuct a `ConformalCubedSphere`

on `domain`

with `ne`

elements across each panel.

`ClimaCore.Meshes.EquiangularCubedSphere`

— Type`EquiangularCubedSphere <: AbstractCubedSphere`

An equiangular gnomonic mesh proposed by Ronchi1996. Uses the element indexing convention of `AbstractCubedSphere`

.

**Constructors**

```
EquiangularCubedSphere(
domain::Domains.SphereDomain,
ne::Integer,
localelementmap=NormalizedBilinearMap()
)
```

Constuct an `EquiangularCubedSphere`

on `domain`

with `ne`

elements across each panel.

`ClimaCore.Meshes.EquidistantCubedSphere`

— Type`EquidistantCubedSphere <: AbstractCubedSphere`

An equidistant gnomonic mesh outlined in Rancic1996 and Nair2005. Uses the element indexing convention of `AbstractCubedSphere`

.

**Constructors**

`EquidistantCubedSphere(domain::Domains.SphereDomain, ne::Integer)`

Constuct an `EquidistantCubedSphere`

on `domain`

with `ne`

elements across each panel.

`ClimaCore.Meshes.ExponentialStretching`

— Type`ExponentialStretching(H::FT)`

Apply exponential stretching to the domain when constructing elements. `H`

is the scale height (a typical atmospheric scale height `H ≈ 7.5`

km).

For an interval $[z_0,z_1]$, this makes the elements uniformally spaced in $\zeta$, where

\[\zeta = \frac{1 - e^{-\eta/h}}{1-e^{-1/h}},\]

where $\eta = \frac{z - z_0}{z_1-z_0}$, and $h = \frac{H}{z_1-z_0}$ is the non-dimensional scale height. If `reverse_mode`

is `true`

, the smallest element is at the top, and the largest at the bottom (this is typical for land model configurations).

Then, the user can define a stretched mesh via

`ClimaCore.Meshes.IntervalMesh(interval_domain, ExponentialStretching(H); nelems::Int, reverse_mode = false)`

`ClimaCore.Meshes.GeneralizedExponentialStretching`

— Type`GeneralizedExponentialStretching(dz_bottom::FT, dz_top::FT)`

Apply a generalized form of exponential stretching to the domain when constructing elements. `dz_bottom`

and `dz_top`

are target element grid spacings at the bottom and at the top of the vertical column domain (m). In typical atmosphere configurations, `dz_bottom`

is the smallest grid spacing and `dz_top`

the largest one. On the other hand, for typical land configurations, `dz_bottom`

is the largest grid spacing and `dz_top`

the smallest one.

For land configurations, use `reverse_mode`

= `true`

(default value `false`

).

Then, the user can define a generalized stretched mesh via

`ClimaCore.Meshes.IntervalMesh(interval_domain, GeneralizedExponentialStretching(dz_bottom, dz_top); nelems::Int, reverse_mode = false)`

`ClimaCore.Meshes.IntervalMesh`

— Type`IntervalMesh <: AbstractMesh`

A 1D mesh on an `IntervalDomain`

.

**Constuctors**

`IntervalMesh(domain::IntervalDomain, faces::AbstractVector)`

Construct a 1D mesh with face locations at `faces`

.

`IntervalMesh(domain::IntervalDomain[, stretching=Uniform()]; nelems=)`

Constuct a 1D mesh on `domain`

with `nelems`

elements, using `stretching`

. Possible values of `stretching`

are:

`ClimaCore.Meshes.IntrinsicMap`

— Type`IntrinsicMap()`

This `LocalElementMap`

uses the intrinsic mapping of the cubed sphere to map the reference element to the physical domain.

`ClimaCore.Meshes.LocalElementMap`

— Type`LocalElementMap`

An abstract type of mappings from the reference element to a physical domain.

`ClimaCore.Meshes.NormalizedBilinearMap`

— Type`NormalizedBilinearMap()`

The `LocalElementMap`

for meshes on spherical domains of Guba2014. It uses bilinear interpolation between the Cartesian coordinates of the element vertices, then normalizes the result to lie on the sphere.

`ClimaCore.Meshes.RectilinearMesh`

— Type`RectilinearMesh <: AbstractMesh2D`

**Constructors**

`RectilinearMesh(domain::RectangleDomain, n1, n2)`

Construct a `RectilinearMesh`

of equally-spaced `n1`

by `n2`

elements on `domain`

.

`RectilinearMesh(intervalmesh1::IntervalMesh1, intervalmesh2::IntervalMesh2)`

Construct the product mesh of `intervalmesh1`

and `intervalmesh2`

.

`ClimaCore.Meshes.SharedVertices`

— Type`Meshes.SharedVertices(mesh, elem, vert)`

An iterator over (element, vertex) pairs that are shared with `(elem,vert)`

.

`ClimaCore.Meshes.Uniform`

— Type`Uniform()`

Use uniformly-sized elements.

`ClimaCore.Meshes.boundary_face_name`

— Function`Meshes.boundary_face_name(mesh::AbstractMesh, elem, face::Int)::Union{Symbol,Nothing}`

The name of the boundary facing `face`

of element `elem`

, or `nothing`

if it is not on the boundary.

`ClimaCore.Meshes.containing_element`

— Function`elem = Meshes.containing_element(mesh::AbstractMesh, coord)`

The element `elem`

in `mesh`

containing the coordinate `coord`

. If the coordinate falls on the boundary between two or more elements, an arbitrary element is chosen.

`ClimaCore.Meshes.containing_panel`

— Method`panel = cubedspherepanel(coord::Geometry.Cartesian123Point)`

Given a point `coord`

, return its panel number (an integer between 1 and 6).

`ClimaCore.Meshes.coordinates`

— Method```
Meshes.coordinates(mesh, elem, vert::Int)
Meshes.coordinates(mesh, elem, ξ::SVector)
```

Return the physical coordinates of a point in an element `elem`

of `mesh`

. The position of the point can either be a vertex number `vert`

or the coordinates `ξ`

in the reference element.

`ClimaCore.Meshes.domain`

— Function`Meshes.domain(mesh::AbstractMesh)`

The domain (a subtype of `Domains.AbstractDomain`

) on which the mesh is defined.

`ClimaCore.Meshes.element_horizontal_length_scale`

— Function`Meshes.element_horizontal_length_scale(mesh::AbstractMesh)`

The approximate length scale (in units of distance) of the elements of the mesh.

`ClimaCore.Meshes.elements`

— Function`Meshes.elements(mesh::AbstractMesh)`

An iterator over the elements of a mesh. Elements of a mesh can be of any type.

`ClimaCore.Meshes.face_connectivity_matrix`

— Function`M = Meshes.face_connectivity_matrix(mesh, elemorder = elements(mesh))`

Construct a `Bool`

-valued `SparseCSCMatrix`

containing the face connections of `mesh`

. Elements are indexed according to `elemorder`

.

Note that `M[i,i] == true`

only if two distinct faces of element `i`

are connected.

`ClimaCore.Meshes.from_panel`

— Method`coord1 = from_panel(panel::Int, coord::Geometry.Cartesian123Point)`

Given a point `coord`

and its panel number `panel`

, return its coordinates in panel 1 (`coord1`

).

`ClimaCore.Meshes.is_boundary_face`

— Function`Meshes.is_boundary_face(mesh::AbstractMesh, elem, face::Int)::Bool`

Determine whether face `face`

of element `elem`

is on the boundary of `mesh`

.

`elem`

should be an element of `elements(mesh)`

.

`ClimaCore.Meshes.linearindices`

— Method`Meshes.linearindices(elemorder)`

Given a data structure `elemorder[i] = elem`

that orders elements, construct the inverse map from `orderindex = linearindices(elemorder)`

such that `orderindex[elem] = i`

.

This will try to use the most efficient structure available.

`ClimaCore.Meshes.nelements`

— Method`nelements(mesh::AbstractMesh)`

The number of elements in the mesh.

`ClimaCore.Meshes.opposing_face`

— Function`opelem, opface, reversed = Meshes.opposing_face(mesh::AbstractMesh, elem, face::Int)`

The element and face (`opelem`

, `opface`

) that oppose face `face`

of element `elem`

.

`ClimaCore.Meshes.reference_coordinates`

— Function`ξ = Meshes.reference_coordinates(mesh::AbstractMesh, elem, coord)`

An `SVector`

of coordinates in the reference element such that

`Meshes.coordinates(mesh, elem, ξ) == coord`

This can be used for interpolation to a specific point.

`ClimaCore.Meshes.refindex`

— Method`i = refindex(ϕ, n)`

Given a reference coordinate `ϕ`

in the interval `[-1, 1]`

`, divide the interval into`

n`evenly spaced subintervals, and return the index of the subinterval (`

i`), and the position in the subinterval (`

ϕs`), normalized to normalized`

[-1, 1]`.

`ClimaCore.Meshes.to_panel`

— Method`to_panel(panel, coord1::Geometry.Cartesian123Point)`

Given a point at `coord1`

on panel 1 of a sphere, transform it to panel `panel`

.

`ClimaCore.Meshes.truncate_mesh`

— Method```
truncate_mesh(
parent_mesh::AbstractMesh,
trunc_domain::IntervalDomain{CT},
)
```

Constructs an `IntervalMesh`

, truncating the given `parent_mesh`

defined on a truncated `trunc_domain`

. The truncation preserves the number of degrees of freedom covering the space from the `trunc_domain`

's `z_bottom`

to `z_top`

, adjusting the stretching.

`ClimaCore.Meshes.vertex_connectivity_matrix`

— Function`M = Meshes.vertex_connectivity_matrix(mesh, elemorder = elements(mesh))`

Construct a `Bool`

-valued `SparseCSCMatrix`

containing the vertex connections of `mesh`

. Elements are indexed according to `elemorder`

.

Note that `M[i,i] == true`

only if two distinct vertices of element `i`

are connected.

`ClimaCore.Topologies.AbstractTopology`

— TypeAbstractTopology

Subtypes of `AbstractHorizontalTopology`

define connectiveness of a mesh in the horizontal domain.

**Interfaces**

`nelems`

`domain(topology::AbstractTopology)`

`mesh`

`nlocalelems`

`nneighbors`

`nsendelems`

`nghostelems`

`localelemindex`

`vertex_coordinates`

`opposing_face`

`face_node_index`

`interior_faces`

`ghost_faces`

`vertex_node_index`

`local_neighboring_elements`

`ghost_neighboring_elements`

`local_vertices`

`ghost_vertices`

`neighbors`

`boundary_tags`

`boundary_tag`

`boundary_faces`

`ClimaCore.Topologies.DSSBuffer`

— Type`DSSBuffer{G, D, A, B}`

**Fields**

`graph_context`

: ClimaComms graph context for communication`perimeter_data`

: Perimeter`DataLayout`

object: typically a`VIFH{TT,Nv,Np}`

, where`TT`

is the transformed type,`Nv`

is the number of vertical levels, and`Np`

is the length of the perimeter

`send_data`

: send buffer`AbstractVector{FT}`

`recv_data`

: recv buffer`AbstractVector{FT}`

`send_buf_idx`

: indexing array for loading send buffer from`perimeter_data`

`recv_buf_idx`

: indexing array for loading (and summing) data from recv buffer to`perimeter_data`

`scalarfidx`

: field id for all scalar fields stored in the`data`

array`covariant12fidx`

: field id for all covariant12vector fields stored in the`data`

array`contravariant12fidx`

: field id for all contravariant12vector fields stored in the`data`

array`internal_elems`

: internal local elements (lidx)`perimeter_elems`

: local elements (lidx) located on process boundary

`ClimaCore.Topologies.IntervalTopology`

— Type`IntervalTopology([context::SingletonCommsContext,] mesh::IntervalMesh)`

A sequential topology on an `Meshes.IntervalMesh`

.

`ClimaCore.Topologies.Perimeter2D`

— Method`Perimeter2D(Nq)`

Construct a perimeter iterator for a 2D spectral element with `Nq`

nodes per dimension (i.e. polynomial degree `Nq-1`

).

`ClimaCore.Topologies.Topology2D`

— Type`Topology2D(mesh::AbstractMesh2D, elemorder=Mesh.elements(mesh))`

This is a distributed topology for 2D meshes. `elemorder`

is a vector or other linear ordering of the `Mesh.elements(mesh)`

. `elempid`

is a sorted vector of the same length as `elemorder`

, each element of which contains the `pid`

of the owning process.

Internally, we can refer to elements in several different ways:

`elem`

: an element of the`mesh`

. Often a`CartesianIndex`

object.`gidx`

: "global index": an enumeration of all elements:`elemorder[gidx] == elem`

`orderindex[elem] == gidx`

`lidx`

: "local index": an enumeration of local elements.`local_elem_gidx[lidx] == gidx`

`sidx`

: "send index": an index into the send buffer of a local element. A single local element may have multiple`sidx`

s if it needs to be send to multiple processes.`send_elem_lidx[sidx] == lidx`

`ridx`

: "receive index": an index into the receive buffer of a ghost element.`recv_elem_gidx[ridx] == gidx`

`ClimaCore.Topologies.VertexIterator`

— Type`Topologies.VertexIterator`

An iterator over the unique (shared) vertices of the topology `topology`

. Each vertex returns a `Vertex`

object, which is itself an iterator.

`ClimaCore.Topologies.boundary_faces`

— Function`boundary_faces(topology, boundarytag)`

An iterator over the faces of `topology`

which face the boundary with tag `boundarytag`

. Each element of the iterator is an `(elem, face)`

pair.

`ClimaCore.Topologies.boundary_tag`

— Function`boundary_tag(topology, name::Symbol)`

The boundary tag of the topology for boundary name `name`

. A boundary tag is an integer that uniquely identifies a boundary.

`ClimaCore.Topologies.boundary_tags`

— Function`boundary_tags(topology)`

A `Tuple`

or `NamedTuple`

of the boundary tags of the topology. A boundary tag is an integer that uniquely identifies a boundary.

`ClimaCore.Topologies.create_dss_buffer`

— Method```
create_dss_buffer(
data::Union{DataLayouts.IJFH{S}, DataLayouts.VIJFH{S}},
topology::Topology2D,
local_geometry = nothing,
local_weights = nothing,
) where {S}
```

Creates a `DSSBuffer`

for the field data corresponding to `data`

`ClimaCore.Topologies.dss!`

— Method`dss!(data, topology)`

Computed unweighted/pure DSS of `data`

.

`ClimaCore.Topologies.dss_ghost!`

— Method```
dss_ghost!(
device::ClimaComms.AbstractCPUDevice,
perimeter_data::DataLayouts.VIFH,
perimeter::AbstractPerimeter,
topology::AbstractTopology,
)
```

Sets the value for all local vertices of each unique ghost vertex, in `perimeter_data`

, to that of the representative ghost vertex.

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.dss_local!`

— Method```
function dss_local!(
::ClimaComms.AbstractCPUDevice,
perimeter_data::DataLayouts.VIFH,
perimeter::AbstractPerimeter,
topology::AbstractTopology,
)
```

Performs DSS on local vertices and faces.

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.dss_local_ghost!`

— Method```
function dss_local_ghost!(
::ClimaComms.AbstractCPUDevice,
perimeter_data::DataLayouts.VIFH,
perimeter::AbstractPerimeter,
topology::AbstractTopology,
)
```

Computes the "local" part of ghost vertex dss. (i.e. it computes the summation of all the shared local vertices of a unique ghost vertex and stores the value in each of the local vertex locations in `perimeter_data`

)

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.dss_local_vertices!`

— Method```
dss_local_vertices!(
perimeter_data::DataLayouts.VIFH,
perimeter::Perimeter2D,
topology::Topology2D,
)
```

Apply dss to local vertices.

`ClimaCore.Topologies.dss_transform!`

— Method```
function dss_transform!(
device::ClimaComms.AbstractDevice,
dss_buffer::DSSBuffer,
data::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
local_geometry::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
weight::DataLayouts.IJFH,
perimeter::AbstractPerimeter,
localelems::Vector{Int},
)
```

Transforms vectors from Covariant axes to physical (local axis), weights the data at perimeter nodes, and stores result in the `perimeter_data`

array. This function calls the appropriate version of `dss_transform!`

based on the data layout of the input arguments.

Arguments:

`dss_buffer`

:`DSSBuffer`

generated by`create_dss_buffer`

function for field data`data`

: field data`local_geometry`

: local metric information defined at each node`weight`

: local dss weights for horizontal space`perimeter`

: perimeter iterator`localelems`

: list of local elements to perform transformation operations on

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.dss_transform!`

— Method```
function dss_transform!(
::ClimaComms.AbstractCPUDevice,
perimeter_data::DataLayouts.VIFH,
data::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
∂ξ∂x::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
∂x∂ξ::Union{DataLayouts.VIJFH, DataLayouts.IJFH},
weight::DataLayouts.IJFH,
perimeter::AbstractPerimeter,
scalarfidx::Vector{Int},
covariant12fidx::Vector{Int},
contravariant12fidx::Vector{Int},
localelems::Vector{Int},
)
```

Transforms vectors from Covariant axes to physical (local axis), weights the data at perimeter nodes, and stores result in the `perimeter_data`

array.

Arguments:

`perimeter_data`

: contains the perimeter field data, represented on the physical axis, corresponding to the full field data in`data`

`data`

: field data`∂ξ∂x`

: partial derivatives of the map from`x`

to`ξ`

:`∂ξ∂x[i,j]`

is ∂ξⁱ/∂xʲ`weight`

: local dss weights for horizontal space`perimeter`

: perimeter iterator`scalarfidx`

: field index for scalar fields in the data layout`covariant12fidx`

: field index for Covariant12 vector fields in the data layout`localelems`

: list of local elements to perform transformation operations on

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.dss_transform`

— Method`dss_transform(arg, local_geometry, weight, I...)`

Transfrom `arg[I...]`

to a basis for direct stiffness summation (DSS). Transformations only apply to vector quantities.

`local_geometry[I...]`

is the relevant`LocalGeometry`

object. If it is`nothing`

, then no transformation is performed`weight[I...]`

is the relevant DSS weights. If`weight`

is`nothing`

, then the result is simply summation.

`ClimaCore.Topologies.dss_untransform!`

— Method```
dss_untransform!(
device::ClimaComms.AbstractDevice,
dss_buffer::DSSBuffer,
data::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
local_geometry::Union{DataLayouts.IJFH, DataLayouts.VIJFH},
perimeter::AbstractPerimeter,
)
```

Transforms the DSS'd local vectors back to Covariant12 vectors, and copies the DSS'd data from the `perimeter_data`

to `data`

. This function calls the appropriate version of `dss_transform!`

function based on the data layout of the input arguments.

Arguments:

`dss_buffer`

:`DSSBuffer`

generated by`create_dss_buffer`

function for field data`data`

: field data`local_geometry`

: local metric information defined at each node`perimeter`

: perimeter iterator`localelems`

: list of local elements to perform transformation operations on

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.dss_untransform`

— Method`dss_untransform(T, targ, local_geometry, I...)`

Transform `targ[I...]`

back to a value of type `T`

after performing direct stiffness summation (DSS).

`ClimaCore.Topologies.face_node_index`

— Function`i,j = face_node_index(face, Nq, q, reversed=false)`

The node indices of the `q`

th node on face `face`

, where `Nq`

is the number of face nodes in each direction.

`ClimaCore.Topologies.fill_send_buffer!`

— Method`fill_send_buffer!(::ClimaComms.AbstractCPUDevice, dss_buffer::DSSBuffer; synchronize=true)`

Loads the send buffer from `perimeter_data`

. For unique ghost vertices, only data from the representative vertices which store result of "ghost local" DSS are loaded.

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.ghost_faces`

— Method`ghost_faces(topology::AbstractTopology)`

An iterator over the ghost faces of `topology`

. Each element of the iterator is a 5-tuple the form

`(elem1, face1, elem2, face2, reversed)`

where `elemX, faceX`

are the element and face numbers, and `reversed`

indicates whether they have opposing orientations.

`ClimaCore.Topologies.ghost_neighboring_elements`

— Function`ghost_neighboring_elements(topology::AbstractTopology, ridx::Integer)`

An iterator of the receive buffer indices (ridx) of the ghost elements which are neighbors of the local element `lidx`

in `topology`

.

`ClimaCore.Topologies.ghost_vertices`

— Function`ghost_vertices(topology)`

An iterator over the ghost vertices of `topology`

. Each vertex is an iterator over `(isghost, lidx/ridx, vert)`

pairs.

`ClimaCore.Topologies.interior_faces`

— Method`interior_faces(topology::AbstractTopology)`

An iterator over the interior faces of `topology`

. Each element of the iterator is a 5-tuple the form

`(elem1, face1, elem2, face2, reversed)`

where `elemX, faceX`

are the element and face numbers, and `reversed`

indicates whether they have opposing orientations.

`ClimaCore.Topologies.load_from_recv_buffer!`

— Method`load_from_recv_buffer!(::ClimaComms.AbstractCPUDevice, dss_buffer::DSSBuffer)`

Adds data from the recv buffer to the corresponding location in `perimeter_data`

. For ghost vertices, this data is added only to the representative vertices. The values are then scattered to other local vertices corresponding to each unique ghost vertex in `dss_local_ghost`

.

Part of `ClimaCore.Spaces.weighted_dss!`

.

`ClimaCore.Topologies.local_neighboring_elements`

— Function`local_neighboring_elements(topology::AbstractTopology, lidx::Integer)`

An iterator of the local element indices (lidx) of the local elements which are neighbors of the local element `lidx`

in `topology`

(excluding `lidx`

itself).

`ClimaCore.Topologies.local_vertices`

— Function`local_vertices(topology)`

An iterator over the interior vertices of `topology`

. Each vertex is an iterator over `(lidx, vert)`

pairs.

`ClimaCore.Topologies.localelemindex`

— Function`localelemindex(topology, elem)`

The local index for the specified element; useful for distributed topologies.

`ClimaCore.Topologies.mesh`

— Function`mesh(topology)`

Returns the mesh underlying the `topology`

`ClimaCore.Topologies.neighbors`

— Function`neighbors(topology)`

Returns an array of the PIDs of the neighbors of this process.

`ClimaCore.Topologies.nelems`

— Function`nelems(topology)`

The total number of elements in `topology`

.

`ClimaCore.Topologies.nghostelems`

— Function`nghostelems(topology)`

The number of ghost elements in `topology`

.

`ClimaCore.Topologies.nlocalelems`

— Function`nlocalelems(topology)`

The number of local elements in `topology`

.

`ClimaCore.Topologies.nneighbors`

— Function`nneighbors(topology)`

The number of neighbors of this process in `topology`

.

`ClimaCore.Topologies.nsendelems`

— Function`nsendelems(topology)`

The number of elements to send to neighbors in `topology`

.

`ClimaCore.Topologies.opposing_face`

— Function`(opelem, opface, reversed) = opposing_face(topology, elem, face)`

The opposing face of face number `face`

of element `elem`

in `topology`

.

`opelem`

is the opposing element number, 0 for a boundary, negative for a ghost element`opface`

is the opposite face number, or boundary face number if a boundary`reversed`

indicates whether the opposing face has the opposite orientation.

`ClimaCore.Topologies.spacefillingcurve`

— Method`spacefillingcurve(mesh::Meshes.AbstractCubedSphere)`

Generate element ordering, `elemorder`

, based on a space filling curve for a `CubedSphere`

mesh.

`ClimaCore.Topologies.spacefillingcurve`

— Method`spacefillingcurve(mesh::Meshes.RectilinearMesh)`

Generate element ordering, `elemorder`

, based on a space filling curve for a `Rectilinear`

mesh.

`ClimaCore.Topologies.vertex_coordinates`

— Function`(c1,c2,c3,c4) = vertex_coordinates(topology, elem)`

The coordinates of the 4 vertices of element `elem`

.

`ClimaCore.Topologies.vertex_node_index`

— Method`i,j = vertex_node_index(vertex_num, Nq)`

The node indices of `vertex_num`

, where `Nq`

is the number of face nodes in each direction.

`ClimaCore.Utilities.half`

— Constant`const half = PlusHalf(0)`

`ClimaCore.Utilities.PlusHalf`

— Type`PlusHalf(i)`

Represents `i + 1/2`

, but stored as internally as an integer value. Used for indexing into staggered finite difference meshes: the convention "half" values are indexed at cell faces, whereas centers are indexed at cell centers.

Supports `+`

, `-`

and inequalities.

See also `half`

.

`ClimaCore.InputOutput.HDF5Reader`

— Type`HDF5Reader(filename::AbstractString[, context::ClimaComms.AbstractCommsContext])`

An `AbstractReader`

for reading from HDF5 files created by `HDF5Writer`

. The reader object contains an internal cache of domains, meshes, topologies and spaces that are read so that duplicate objects are not created.

The optional `context`

can be used for reading distributed fields: in this case, the `MPICommsContext`

used passed as an argument: resulting `Field`

s will be distributed using this context. As with `HDF5Writer`

, this requires a HDF5 library with MPI support.

**Interface**

**Usage**

```
reader = InputOutput.HDF5Reader(filename)
Y = read_field(reader, "Y")
Y.c |> propertynames
Y.f |> propertynames
ρ_field = read_field(reader, "Y.c.ρ")
w_field = read_field(reader, "Y.f.w")
close(reader)
```

To explore the contents of the `reader`

, use either

`julia> reader |> propertynames`

e.g, to explore the components of the `space`

,

```
julia> reader.space_cache
Dict{Any, Any} with 3 entries:
"center_extruded_finite_difference_space" => CenterExtrudedFiniteDifferenceSpace:…
"horizontal_space" => SpectralElementSpace2D:…
"face_extruded_finite_difference_space" => FaceExtrudedFiniteDifferenceSpace:…
```

Once "unpacked" as shown above, `ClimaCorePlots`

or `ClimaCoreMakie`

can be used to visualise fields. `ClimaCoreTempestRemap`

supports interpolation onto user-specified grids if necessary.

`ClimaCore.InputOutput.HDF5Writer`

— Type```
HDF5Writer(filename::AbstractString[, context::ClimaComms.AbstractCommsContext];
overwrite::Bool = true)
```

An `AbstractWriter`

for writing to HDF5-formatted files using the ClimaCore storage conventions. An internal cache is used to avoid writing duplicate domains, meshes, topologies and spaces to the file. Use `HDF5Reader`

to load the data from the file.

The optional `context`

can be used for writing distributed fields: in this case, the `MPICommsContext`

used passed as an argument: this must match the context used for distributing the `Field`

.

The writer overwrites or appends to existing files depending on the value of the `overwrite`

keyword argument. When `overwrite`

is `false`

, the writer appends to `filename`

if the file already exists, otherwise it creates a new one.

The default Julia HDF5 binaries are not built with MPI support. To use the distributed functionality, you will need to configure HDF5.jl with an MPI-enabled HDF5 library, see the HDF5.jl documentation.

**Interface**

**Usage**

```
writer = InputOutput.HDF5Writer(filename)
InputOutput.write!(writer, Y, "Y")
close(writer)
```

`ClimaCore.InputOutput.defaultname`

— Function`defaultname(obj)`

Default name of object for InputOutput writers.

`ClimaCore.InputOutput.matrix_to_cartesianindices`

— Method`matrix_to_cartesianindices(elemorder_matrix)`

Converts the `elemorder_matrix`

to cartesian indices.

`ClimaCore.InputOutput.read_attributes`

— Method`read_attributes(reader::AbstractReader, name::AbstractString, data::Dict)`

Return the attributes associated to the object at `name`

in the given HDF5 file.

`ClimaCore.InputOutput.read_domain`

— Method`read_domain(reader::AbstractReader, name)`

Reads a domain named `name`

from `reader`

. Domain objects are cached in the reader to avoid creating duplicate objects.

`ClimaCore.InputOutput.read_field`

— Method`read_field(reader, name)`

Reads a `Field`

or `FieldVector`

named `name`

from `reader`

. Fields are *not* cached, so that reading the same field multiple times will create multiple distinct objects.

`ClimaCore.InputOutput.read_grid`

— Method`read_grid(reader::AbstractReader, name)`

Reads a space named `name`

from `reader`

, or from the reader cache if it has already been read.

`ClimaCore.InputOutput.read_mesh`

— Method`read_mesh(reader::AbstractReader, name)`

Reads a mesh named `name`

from `reader`

, or from the reader cache if it has already been read.

`ClimaCore.InputOutput.read_space`

— Method`read_space(reader::AbstractReader, name)`

Reads a space named `name`

from `reader`

, or from the reader cache if it has already been read.

`ClimaCore.InputOutput.read_topology`

— Method`read_topology(reader::AbstractReader, name)`

Reads a topology named `name`

from `reader`

, or from the reader cache if it has already been read.

`ClimaCore.InputOutput.write!`

— Function`write!(writer::AbstractWriter, obj[, preferredname])`

Write the object `obj`

using `writer`

. An optional `preferredname`

can be provided, otherwise `defaultname`

will be used to generate a name. The name of the object will be returned.

A cache of domains, meshes, topologies and spaces is kept: if one of these objects has already been written, then the file will not be modified: instead the name under which the object was first written will be returned. Note that `Field`

s and `FieldVector`

s are *not* cached, and so can be written multiple times.

`ClimaCore.InputOutput.write!`

— Method`write!(filename::AbstractString, name => value...)`

Write one or more `name => value`

pairs to the HDF5 file `filename`

.

`ClimaCore.InputOutput.write!`

— Method`write!(writer::HDF5Writer, name => value...)`

Write one or more `name => value`

pairs to `writer`

.

`ClimaCore.InputOutput.write_attributes!`

— Method`write_attributes!(writer::AbstractWriter, name::AbstractString, data::Dict)`

Write `data`

as attributes to the object at `name`

in the given HDF5 file.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, mesh, name)`

Write `CubedSphereMesh`

data to HDF5.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, topology, name)`

Write `IntervalTopology`

data to HDF5.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, domain, name)`

Writes an object of type 'SphereDomain' and name 'name' to the HDF5 file.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, topology, name)`

Write `Topology2D`

data to HDF5.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, domain, name)`

Writes an object of type 'Hypsography' and name 'name' to the HDF5 file.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, space, name)`

Write `SpectralElementSpace1D`

data to HDF5.

`ClimaCore.InputOutput.write_new!`

— Function`write_new!(writer, domain, name)`

Writes an object of type 'IntervalDomain' and name 'name' to the HDF5 file.

`ClimaCore.DeviceSideContext`

— Type`DeviceSideContext()`

Context associated to data defined on the device side of an accelerator.

The most common example is data defined on a GPU. DeviceSideContext() is used for operations within the accelerator.

`ClimaCore.DeviceSideDevice`

— Type`DeviceSideDevice()`

This device represents data defined on the device side of an accelerator.

The most common example is data defined on a GPU. DeviceSideDevice() is used for operations within the accelerator.

`ClimaCore.column`

— Function`column(data::AbstractData, i::Integer)`

A contiguous "column" view into an underlying data layout `data`

at nodal point index `i`

.

`ClimaCore.slab`

— Function`slab(data::AbstractData, h::Integer)`

A "pancake" view into an underlying data layout `data`

at location `h`

.

`ClimaCore.Operators.AbstractBoundaryCondition`

— Type`AbstractBoundaryCondition`

An abstract type for boundary conditions for `FiniteDifferenceOperator`

s.

Subtypes should define:

`ClimaCore.Operators.AbstractOperator`

— Type`AbstractOperator`

Supertype for ClimaCore operators. An operator is a pseudo-function: it can't be called directly, but can be broadcasted over `Field`

s.

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

— Type`CentralNumericalFlux(fluxfn)`

Evaluates the central numerical flux using `fluxfn`

.

`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**

- Taylor2010, equation 17

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

.

`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**

- Taylor2010, equation 15

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

— Type`Extrapolate()`

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

`ClimaCore.Operators.FCTBorisBook`

— Type```
U = FCTBorisBook(;boundaries)
U.(v, x)
```

Correct the flux using the flux-corrected transport formulation by Boris and Book BorisBook1973.

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 BorisBook1973, as reported in durran2010 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 zalesak1979fully.

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 zalesak1979fully, as reported in durran2010 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.FiniteDifferenceOperator`

— Type`FiniteDifferenceOperator`

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

See also `AbstractBoundaryCondition`

for how to define the boundaries.

`ClimaCore.Operators.FirstOrderOneSided`

— Type`FirstOrderOneSided()`

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

`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**

- Taylor2010, equation 16

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

— Type```
L = LeftBiased3rdOrderC2F(;boundaries)
L.(x)
```

Interpolate a center-value field to a face-valued field from the left, using a 3rd-order reconstruction.

\[L(x)[i] = \left(-2 x[i-\tfrac{3}{2}] + 10 x[i-\tfrac{1}{2}] + 4 x[i+\tfrac{1}{2}] \right) / 12\]

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

— Type```
L = LeftBiased3rdOrderF2C(;boundaries)
L.(x)
```

Interpolate a face-value field to a center-valued field from the left, using a 3rd-order reconstruction.

\[L(x)[i+\tfrac{1}{2}] = \left(-2 x[i-1] + 10 x[i] + 4 x[i+1] \right) / 12\]

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

— Method`LinearRemap(target::AbstractSpace, source::AbstractSpace)`

A remapping operator from the `source`

space to the `target`

space.

See Ullrich2015 eqs. 3 and 4.

`ClimaCore.Operators.NullBoundaryCondition`

— Type`NullBoundaryCondition()`

This is used as a placeholder when no other boundary condition can be applied.

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

`ClimaCore.Operators.RightBiased3rdOrderC2F`

— Type```
R = RightBiased3rdOrderC2F(;boundaries)
R.(x)
```

Interpolate a center-valued field to a face-valued field from the right, using a 3rd-order reconstruction.

\[R(x)[i] = \left(4 x[i-\tfrac{1}{2}] + 10 x[i+\tfrac{1}{2}] -2 x[i+\tfrac{3}{2}] \right) / 12\]

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

— Type```
R = RightBiased3rdOrderF2C(;boundaries)
R.(x)
```

Interpolate a face-valued field to a center-valued field from the right, using a 3rd-order reconstruction.

\[R(x)[i] = \left(4 x[i] + 10 x[i+1] -2 x[i+2] \right) / 12\]

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

`ClimaCore.Operators.RusanovNumericalFlux`

— Type`RusanovNumericalFlux(fluxfn, wavespeedfn)`

Evaluates the Rusanov numerical flux using `fluxfn`

with wavespeed `wavespeedfn`

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

— Type`SetCurl(val)`

Set the curl at the boundary to be `val`

.

`ClimaCore.Operators.SetDivergence`

— Type`SetDivergence(val)`

Set the divergence at the boundary to be `val`

.

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

— Type`SlabBlockSpectralStyle()`

Applies spectral-element operations using by making use of intermediate temporaries for each operator. This is used for CPU kernels.

`ClimaCore.Operators.SpectralBroadcasted`

— Type`SpectralBroadcasted{Style}(op, args[,axes[, work]])`

This is similar to a `Base.Broadcast.Broadcasted`

object, except it contains space for an intermediate `work`

storage.

This is returned by `Base.Broadcast.broadcasted(op::SpectralElementOperator)`

.

`ClimaCore.Operators.SpectralElementOperator`

— Type`SpectralElementOperator`

Represents an operation that is applied to each element.

Subtypes `Op`

of this should define the following:

`operator_return_eltype(::Op, ElTypes...)`

`allocate_work(::Op, args...)`

`apply_operator(::Op, work, args...)`

Additionally, the result type `OpResult <: OperatorSlabResult`

of `apply_operator`

should define `get_node(::OpResult, ij, slabidx)`

.

`ClimaCore.Operators.SpectralStyle`

— Type`SpectralStyle()`

Broadcasting requires use of spectral-element operations.

`ClimaCore.Operators.StencilBroadcasted`

— Type`StencilBroadcasted{Style}(op, args[,axes[, work]])`

This is similar to a `Base.Broadcast.Broadcasted`

object.

This is returned by `Base.Broadcast.broadcasted(op::FiniteDifferenceOperator)`

.

`ClimaCore.Operators.ThirdOrderOneSided`

— Type`ThirdOrderOneSided()`

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

`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 WickerSkamarock2002, 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.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.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\]

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

— Method`_resolve_operator(slabidx, args...)`

Calls `resolve_operator(arg, slabidx)`

for each `arg`

in `args`

`ClimaCore.Operators.add_numerical_flux_internal!`

— Method`add_numerical_flux_internal!(fn, dydt, args...)`

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

The numerical flux is determined by evaluating

`fn(normal, argvals⁻, argvals⁺)`

where:

`normal`

is the unit normal vector, pointing from the "minus" side to the "plus" side`argvals⁻`

is the tuple of values of`args`

on the "minus" side of the face`argvals⁺`

is the tuple of values of`args`

on the "plus" side of the face

and should return the net flux from the "minus" side to the "plus" side.

For consistency, it should satisfy the property that

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

See also:

`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.column_integral_definite!`

— Method`column_integral_definite!(∫field::Field, ᶜfield::Field)`

Sets `∫field`

${}= \int_0^{z_{max}}\,$`ᶜfield`

$(z)\,dz$, where $z_{max}$ is the value of `z`

at the top of the domain. The input `ᶜfield`

must lie on a cell-center space, and the output `∫field`

must lie on the corresponding horizontal space.

`ClimaCore.Operators.column_integral_indefinite!`

— Method`column_integral_indefinite!(ᶠ∫field::Field, ᶜfield::Field)`

Sets `ᶠ∫field`

$(z) = \int_0^z\,$`ᶜfield`

$(z')\,dz'$. The input `ᶜfield`

must lie on a cell-center space, and the output `ᶠ∫field`

must lie on the corresponding cell-face space.

```
column_integral_indefinite!(
f::Function,
ᶠ∫field::Fields.ColumnField,
ϕ₀ = 0,
average = (ϕ⁻, ϕ⁺) -> (ϕ⁻ + ϕ⁺) / 2,
)
```

The indefinite integral `ᶠ∫field = ϕ(z) = ∫ f(ϕ,z) dz`

given:

`f`

the integral integrand function (which may be a function)`ᶠ∫field`

the resulting (scalar) field`ϕ(z)`

`ϕ₀`

(optional) the boundary condition`average`

(optional) a function to compute the cell center average between two cell faces (`ϕ⁻, ϕ⁺`

).

`ClimaCore.Operators.column_mapreduce!`

— Method`column_mapreduce!(fn, op, reduced_field::Field, fields::Field...)`

Applies mapreduce along the vertical direction. The input `fields`

must all lie on the same space, and the output `reduced_field`

must lie on the corresponding horizontal space. The function `fn`

is mapped over every input, and the function `op`

is used to reduce the outputs of `fn`

.

`ClimaCore.Operators.column_thomas_solve!`

— Method`column_thomas_solve!(A, b)`

Solves the linear system `A * x = b`

, where `A`

is a tri-diagonal matrix (represented by a `Field`

of tri-diagonal matrix rows), and where `b`

is a vector (represented by a `Field`

of numbers). The data in `b`

is overwritten with the solution `x`

, and the upper diagonal of `A`

is also overwritten with intermediate values used to compute `x`

. The algorithm is described here: https://en.wikipedia.org/wiki/Tridiagonal*matrix*algorithm.

`ClimaCore.Operators.copyto_slab!`

— Method`copyto_slab!(out, bc, slabidx)`

Copy the slab indexed by `slabidx`

from `bc`

to `out`

.

`ClimaCore.Operators.is_valid_index`

— Method`is_valid_index(space, ij, slabidx)::Bool`

Returns `true`

if the node indices `ij`

and slab indices `slabidx`

are valid for `space`

.

`ClimaCore.Operators.left_interior_idx`

— Method`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.left_interior_window_idx`

— Method`left_interior_window_idx(arg, space, loc)`

Compute the index of the leftmost point which uses only the interior stencil of the space.

`ClimaCore.Operators.linear_remap_op`

— Method`linear_remap_op(target::AbstractSpace, source::AbstractSpace)`

Computes linear remapping operator `R`

for remapping from `source`

to `target`

spaces.

Entry `R_{ij}`

gives the contribution weight to the target node `i`

from source node `j`

; nodes are indexed by their global position, determined by both element index and nodal order within the element.

`ClimaCore.Operators.local_weights`

— Method`local_weights(space::AbstractSpace)`

Each degree of freedom is associated with a local weight J*i. For finite volumes the local weight J*i would represent the geometric area of the associated region. For nodal finite elements, the local weight represents the value of the global Jacobian, or some global integral of the associated basis function.

See [Ullrich2015] section 2.

`ClimaCore.Operators.matrix_interpolate`

— Function`matrix_interpolate(field, quadrature)`

Computes the tensor product given a uniform quadrature `out = (M ⊗ M) * in`

on each element. Returns a 2D Matrix for plotting / visualizing 2D Fields.

`ClimaCore.Operators.matrix_interpolate`

— Method`matrix_interpolate(field, Nu::Integer)`

Computes the tensor product given a uniform quadrature degree of Nu on each element. Returns a 2D Matrix for plotting / visualizing 2D Fields.

`ClimaCore.Operators.operator_axes`

— Function`operator_axes(space)`

Return a tuple of the axis indicies a given field operator works over.

`ClimaCore.Operators.overlap`

— Method`overlap(target, source)`

Computes local weights of the overlap mesh for `source`

to `target`

spaces.

`ClimaCore.Operators.overlap_weights!`

— Method`overlap_weights!(J_ov, target, source, i, j, coords_t, coords_s)`

Computes the overlap weights for a pair of source and target elements.

The spatial overlap of element `i`

on the `target`

space and element `j`

of the `source`

space, computed for each nodal pair, approximating Ullrich2015 eq. 19 via quadrature.

`ClimaCore.Operators.reconstruct_placeholder_broadcasted`

— Method`reconstruct_placeholder_broadcasted(space, obj)`

Recurively reconstructs objects that have been stripped via `strip_space`

.

`ClimaCore.Operators.remap!`

— Method`remap!(target_field::Field, R::LinearRemap, source_field::Field)`

Applies the remapping operator `R`

to `source_field`

and stores the solution in `target_field`

.

`ClimaCore.Operators.remap`

— Method`remap(R::LinearRemap, source_field::Field)`

Applies `R`

to `source_field`

and outputs a new field on the target space.

`ClimaCore.Operators.resolve_operator`

— Method`resolve_operator(bc, slabidx)`

Recursively evaluate any operators in `bc`

at `slabidx`

, replacing any `SpectralBroadcasted`

objects.

- if
`bc`

is a regular`Broadcasted`

object, return a new`Broadcasted`

with`resolve_operator`

called on each`arg`

- if
`bc`

is a regular`SpectralBroadcasted`

object: - call
`resolve_operator`

called on each`arg`

- call
`apply_operator`

, returning the resulting "pseudo Field": a`Field`

with an

`IF`

/`IJF`

data object.

- if
`bc`

is a`Field`

, return that

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

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

`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.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_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.tensor_product!`

— Function```
tensor_product!(out, in, M)
tensor_product!(inout, M)
```

Computes the tensor product `out = (M ⊗ M) * in`

on each element.

`ClimaCore.Operators.unionall_type`

— Method`unionall_type(::Type{T})`

Extract the type of the input, and strip it of any type parameters.

`ClimaCore.Operators.x_overlap`

— Method`x_overlap(target, source)`

Computes 1D local weights of the overlap mesh for `source`

to `target`

spaces.

For 1D spaces, this returns the full local weight matrix of the overlap. In 2D, this returns the overlap weights in the first dimension.

`ClimaCore.Operators.y_overlap`

— Method`y_overlap(target, source)`

Computes 1D local weights of the overlap mesh for `source`

to `target`

spaces.

This returns the overlap weights in the second dimension for 2D spaces.

`ClimaCore.DataLayouts`

— Module`ClimaCore.DataLayouts`

Defines the following DataLayouts (see individual docs for more info):

TODO: Add links to these datalayouts

`IJKFVH`

`IJFH`

`IFH`

`DataF`

`IJF`

`IF`

`VF`

`VIJFH`

`VIFH`

`IH1JH2`

`IV1JH2`

Notation:

`i,j`

are horizontal node indices within an element`k`

is the vertical node index within an element`f`

is the field index (1 if field is scalar, >1 if it is a vector field)`v`

is the vertical element index in a stack`h`

is the element stack index

Data layout is specified by the order in which they appear, e.g. `IJKFVH`

indexes the underlying array as `[i,j,k,f,v,h]`

`ClimaCore.DataLayouts.Data0D`

— Type`Data0D{S}`

`ClimaCore.DataLayouts.Data1D`

— Type`Data1D{S,Ni}`

Abstract type for data storage for a 1D field made up of `Ni`

values of type `S`

.

Objects `data`

should define `slab(data, h)`

to return a `DataSlab2D{S,Nij}`

object.

`ClimaCore.DataLayouts.Data1DX`

— Type`Data1DX{S, Nv, Ni}`

Abstract type for data storage for a 1D field with extruded columns. The horizontal is made up of `Ni`

values of type `S`

.

Objects `data`

should define `slab(data, v, h)`

to return a `DataSlab1D{S,Ni}`

object, and a `column(data, i, h)`

to return a `DataColumn`

.

`ClimaCore.DataLayouts.Data2D`

— Type`Data2D{S,Nij}`

Abstract type for data storage for a 2D field made up of `Nij × Nij`

values of type `S`

.

Objects `data`

should define `slab(data, h)`

to return a `DataSlab2D{S,Nij}`

object.

`ClimaCore.DataLayouts.Data2DX`

— Type`Data2DX{S,Nv,Nij}`

Abstract type for data storage for a 2D field with extruded columns. The horizontal is made is made up of `Nij × Nij`

values of type `S`

.

Objects `data`

should define `slab(data, v, h)`

to return a `DataSlab2D{S,Nv,Nij}`

object, and a `column(data, i, j, h)`

to return a `DataColumn`

.

`ClimaCore.DataLayouts.Data3D`

— Type`Data3D{S,Nij,Nk}`

Abstract type for data storage for a 3D field made up of `Nij × Nij × Nk`

values of type `S`

.

`ClimaCore.DataLayouts.DataColumn`

— Type`DataColumn{S, Nv}`

Abstract type for data storage for a column. Objects `data`

should define a `data[k,v]`

, returning a value of type `S`

.

`ClimaCore.DataLayouts.DataF`

— Type`DataF{S, A} <: Data0D{S}`

Backing `DataLayout`

for 0D point data.

`ClimaCore.DataLayouts.DataSlab1D`

— Type`DataSlab1D{S,Ni}`

Abstract type for data storage for a slab of `Ni`

values of type `S`

. Objects `data`

should define a `data[i]`

, returning a value of type `S`

.

`ClimaCore.DataLayouts.DataSlab2D`

— Type`DataSlab2D{S,Nij}`

Abstract type for data storage for a slab of `Nij × Nij`

values of type `S`

. Objects `data`

should define a `data[i,j]`

, returning a value of type `S`

.

`ClimaCore.DataLayouts.IF`

— Type`IF{S, Ni, A} <: DataSlab1D{S, Ni}`

Backing `DataLayout`

for 1D spectral element slab data.

Nodal element data (I) are contiguous for each `S`

datatype struct field (F) for a single element slab.

A `DataSlab1D`

view can be returned from other `Data1D`

objects by calling `slab(data, idx...)`

.

`ClimaCore.DataLayouts.IFH`

— Type```
IFH{S, Ni, A} <: Data1D{S, Ni}
IFH{S,Ni}(ArrayType, nelements)
```

Backing `DataLayout`

for 1D spectral element slabs.

Element nodal point (I) data is contiguous for each datatype `S`

struct field (F), for each 1D mesh element (H).

The `ArrayType`

-constructor makes a IFH 1D Spectral DataLayout given the backing `ArrayType`

, quadrature degrees of freedom `Ni`

, and the number of mesh elements `nelements`

.

`ClimaCore.DataLayouts.IH1JH2`

— Type`IH1JH2{S, Nij}(data::AbstractMatrix{S})`

Stores a 2D field in a matrix using a column-major format. The primary use is for interpolation to a regular grid for ex. plotting / field output.

`ClimaCore.DataLayouts.IJF`

— Type`IJF{S, Nij, A} <: DataSlab2D{S, Nij}`

Backing `DataLayout`

for 2D spectral element slab data.

Nodal element data (I,J) are contiguous for each `S`

datatype struct field (F) for a single element slab.

A `DataSlab2D`

view can be returned from other `Data2D`

objects by calling `slab(data, idx...)`

.

`ClimaCore.DataLayouts.IJFH`

— Type```
IJFH{S, Nij, A} <: Data2D{S, Nij}
IJFH{S,Nij}(ArrayType, nelements)
```

Backing `DataLayout`

for 2D spectral element slabs.

Element nodal point (I,J) data is contiguous for each datatype `S`

struct field (F), for each 2D mesh element slab (H).

The `ArrayType`

-constructor constructs a IJFH 2D Spectral DataLayout given the backing `ArrayType`

, quadrature degrees of freedom `Nij × Nij`

, and the number of mesh elements `nelements`

.

`ClimaCore.DataLayouts.IJKFVH`

— Type`IJKFVH{S, Nij, Nk}(array::AbstractArray{T, 6}) <: Data3D{S, Nij, Nk}`

A 3D DataLayout. TODO: Add more docs

`ClimaCore.DataLayouts.IV1JH2`

— Type`IV1JH2{S, n1, Ni}(data::AbstractMatrix{S})`

Stores values from an extruded 1D spectral field in a matrix using a column-major format. The primary use is for interpolation to a regular grid for ex. plotting / field output.

`ClimaCore.DataLayouts.VF`

— Type`VF{S, A} <: DataColumn{S, Nv}`

Backing `DataLayout`

for 1D FV column data.

Column level data (V) are contiguous for each `S`

datatype struct field (F).

A `DataColumn`

view can be returned from other `Data1DX`

, `Data2DX`

objects by calling `column(data, idx...)`

.

`ClimaCore.DataLayouts.VIFH`

— Type`VIFH{S, Nv, Ni, A} <: Data1DX{S, Nv, Ni}`

Backing `DataLayout`

for 1D spectral element slab + extruded 1D FV column data.

Column levels (V) are contiguous for every element nodal point (I) for each datatype `S`

struct field (F), for each 1D mesh element slab (H).

`ClimaCore.DataLayouts.VIJFH`

— Type`VIJFH{S, Nij, A} <: Data2DX{S, Nij}`

Backing `DataLayout`

for 2D spectral element slab + extruded 1D FV column data.

Column levels (V) are contiguous for every element nodal point (I, J) for each `S`

datatype struct field (F), for each 2D mesh element slab (H).

`ClimaCore.DataLayouts.array2data`

— Function`array2data(array, ::AbstractData)`

Reshapes `array`

(of scalars) to fit into the given `DataLayout`

.

The dimensions of `array`

are assumed to be

`([number of vertical nodes], number of horizontal nodes)`

.

`ClimaCore.DataLayouts.bypass_constructor`

— Method`StructArrays.bypass_constructor(T, args)`

Create an instance of type `T`

from a tuple of field values `args`

, bypassing possible internal constructors. `T`

should be a concrete type.

`ClimaCore.DataLayouts.check_basetype`

— Method`check_basetype(::Type{T}, ::Type{S})`

Check whether the types `T`

and `S`

have well-defined sizes, and whether an object of type `S`

can be stored in an array with elements of type `T`

.

`ClimaCore.DataLayouts.data2array`

— Function`data2array(::AbstractData)`

Reshapes the DataLayout's parent array into a `Vector`

, or (for DataLayouts with vertical levels) `Nv x N`

matrix, where `Nv`

is the number of vertical levels and `N`

is the remaining dimensions.

The dimensions of the resulting array are

`([number of vertical nodes], number of horizontal nodes)`

.

Also, this assumes that `eltype(data) <: Real`

.

`ClimaCore.DataLayouts.fieldtypeoffset`

— Method`fieldtypeoffset(T,S,i)`

Similar to `fieldoffset(S,i)`

, but gives result in multiples of `sizeof(T)`

instead of bytes.

`ClimaCore.DataLayouts.get_struct`

— Method`get_struct(array, S, Val(D), start_index)`

Construct an object of type `S`

packed along the `D`

dimension, from the values of `array`

, starting at `start_index`

.

`ClimaCore.DataLayouts.is_valid_basetype`

— Method`is_valid_basetype(::Type{T}, ::Type{S})`

Determines whether an object of type `S`

can be stored in an array with elements of type `T`

by recursively checking whether the non-empty fields of `S`

can be stored in such an array. If `S`

is empty, this is always true.

`ClimaCore.DataLayouts.parent_array_type`

— Method`parent_array_type(::Type{<:AbstractArray})`

Returns the parent array type underlying any wrapper types, with all dimensionality information removed.

`ClimaCore.DataLayouts.promote_parent_array_type`

— Method`promote_parent_array_type(::Type{<:AbstractArray}, ::Type{<:AbstractArray})`

Given two parent array types (without any dimensionality information), promote both the element types and the array types themselves.

`ClimaCore.DataLayouts.replace_basetype`

— Method`replace_basetype(::Type{T}, ::Type{T′}, ::Type{S})`

Changes the type parameters of `S`

to produce a new type `S′`

such that, if `is_valid_basetype(T, S)`

is true, then so is `is_valid_basetype(T′, S′)`

.

`ClimaCore.DataLayouts.set_struct!`

— Method`set_struct!(array, val::S, Val(D), start_index)`

Store an object `val`

of type `S`

packed along the `D`

dimension, into `array`

, starting at `start_index`

.

`ClimaCore.DataLayouts.typesize`

— Method`typesize(T,S)`

Similar to `sizeof(S)`

, but gives the result in multiples of `sizeof(T)`

.

`ClimaCore.Remapping.Remapper`

— MethodRemapper(space, target*hcoords, target*zcoords; buffer*length = 1) Remapper(space, target*hcoords; buffer_length = 1)

Return a `Remapper`

responsible for interpolating any `Field`

defined on the given `space`

to the Cartesian product of `target_hcoords`

with `target_zcoords`

.

`target_zcoords`

can be `nothing`

for interpolation on horizontal spaces.

The `Remapper`

is designed to not be tied to any particular `Field`

. You can use the same `Remapper`

for any `Field`

as long as they are all defined on the same `topology`

.

`Remapper`

is the main argument to the `interpolate`

function.

**Keyword arguments**

`buffer_length`

is size of the internal buffer in the Remapper to store intermediate values for interpolation. Effectively, this controls how many fields can be remapped simultaneously in `interpolate`

. When more fields than `buffer_length`

are passed, the remapper will batch the work in sizes of `buffer_length`

.

`ClimaCore.Remapping._apply_mpi_bitmask!`

— Method`_apply_mpi_bitmask!(remapper::Remapper, num_fields::Int)`

Change to local (private) state of the `remapper`

by applying the MPI bitmask and reconstructing the correct shape for the interpolated values.

Internally, `remapper`

performs interpolation on a flat list of points, this function moves points around according to MPI-ownership and the expected output shape.

`num_fields`

is the number of fields that have been processed and have to be moved in the `interpolated_values`

. We assume that it is always the first `num_fields`

that have to be moved.

`ClimaCore.Remapping._collect_and_return_interpolated_values!`

— Method```
_collect_and_return_interpolated_values!(remapper::Remapper,
num_fields::Int)
```

Perform an MPI call to aggregate the interpolated points from all the MPI processes and save the result in the local state of the `remapper`

. Only the root process will return the interpolated data.

`_collect_and_return_interpolated_values!`

is type-unstable and allocates new return arrays.

`num_fields`

is the number of fields that have been interpolated in this batch.

`ClimaCore.Remapping._reset_interpolated_values!`

— Method`_reset_interpolated_values!(remapper::Remapper)`

Reset the local (private) state in `remapper`

. This function has to be called before performing interpolation.

`ClimaCore.Remapping._set_interpolated_values!`

— Method`_set_interpolated_values!(remapper, field)`

Change the local state of `remapper`

by performing interpolation of `fields`

on the vertical and horizontal points.

`ClimaCore.Remapping.batched_ranges`

— Method`batched_ranges(num_fields, buffer_length)`

Partition the indices from 1 to num*fields in such a way that no range is larger than buffer*length.

`ClimaCore.Remapping.containing_pid`

— Method`containing_pid(target_point, topology)`

Return the process id that contains the `target_point`

in the given `topology`

.

`ClimaCore.Remapping.interpolate`

— Method` interpolate(field::ClimaCore.Fields, target_hcoords, target_zcoords)`

Interpolate the given fields on the Cartesian product of `target_hcoords`

with `target_zcoords`

(if not empty).

Coordinates have to be `ClimaCore.Geometry.Points`

.

Note: do not use this method when performance is important. Instead, define a `Remapper`

and call `interpolate(remapper, fields)`

. Different `Field`

s defined on the same `Space`

can share a `Remapper`

, so that interpolation can be optimized.

**Example**

Given `field`

, a `Field`

defined on a cubed sphere.

```
longpts = range(-180.0, 180.0, 21)
latpts = range(-80.0, 80.0, 21)
zpts = range(0.0, 1000.0, 21)
hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
zcoords = [Geometry.ZPoint(z) for z in zpts]
interpolate(field, hcoords, zcoords)
```

`ClimaCore.Remapping.interpolate`

— Methodinterpolate(remapper::Remapper, fields) interpolate!(dest, remapper::Remapper, fields)

Interpolate the given `field`

(s) as prescribed by `remapper`

.

The optimal number of fields passed is the `buffer_length`

of the `remapper`

. If more fields are passed, the `remapper`

will batch work with size up to its `buffer_length`

.

This call mutates the internal (private) state of the `remapper`

.

Horizontally, interpolation is performed with the barycentric formula in Berrut2004, equation (3.2). Vertical interpolation is linear except in the boundary elements where it is 0th order.

`interpolate!`

writes the output to the given `dest`

iniation. `dest`

is expected to be defined on the root process and to be `nothing`

for the other processes.

Note: `interpolate`

allocates new arrays and has some internal type-instability, `interpolate!`

is non-allocating and type-stable.

When using `interpolate!`

, the `dest`

ination has to be the same array type as the device in use (e.g., `CuArray`

for CUDA runs).

**Example**

Given `field1`

,`field2`

, two `Field`

defined on a cubed sphere.

```
longpts = range(-180.0, 180.0, 21)
latpts = range(-80.0, 80.0, 21)
zpts = range(0.0, 1000.0, 21)
hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
zcoords = [Geometry.ZPoint(z) for z in zpts]
space = axes(field1)
remapper = Remapper(space, hcoords, zcoords)
int1 = interpolate(remapper, field1)
int2 = interpolate(remapper, field2)
# Or
int12 = interpolate(remapper, [field1, field2])
# With int1 = int12[1, :, :, :]
```

`ClimaCore.Remapping.interpolate_array`

— Function```
interpolate_array(field, xpts, ypts)
interpolate_array(field, xpts, ypts, zpts)
```

Interpolate a field to a regular array using pointwise interpolation.

This is primarily used for plotting and diagnostics.

**Examples**

```
longpts = range(Geometry.LongPoint(-180.0), Geometry.LongPoint(180.0), length = 21)
latpts = range(Geometry.LatPoint(-80.0), Geometry.LatPoint(80.0), length = 21)
zpts = range(Geometry.ZPoint(0.0), Geometry.ZPoint(1000.0), length = 21)
interpolate_array(field, longpts, latpts, zpts)
```

Hypsography is not currently handled correctly.

`ClimaCore.Remapping.interpolate_column`

— Method`interpolate_column(field, zpts, weights, gidx)`

Interpolate the given `field`

on the given points assuming the given interpolation_matrix and global index in the topology.

The coefficients `weights`

are computed with `Quadratures.interpolation_matrix`

. See also `interpolate_array`

.

**Keyword arguments**

`physical_z`

: When`true`

, the given`zpts`

are interpreted as "from mean sea level" (instead of "from surface") and hypsography is interpolated.`NaN`

s are returned for values that are below the surface.

`ClimaCore.Remapping.interpolate_slab!`

— Method`interpolate_slab!(output_array, field, slab_indices, weights)`

Interpolate horizontal field on the given `slab_indices`

using the given interpolation `weights`

.

`interpolate_slab!`

interpolates several values at a fixed `z`

coordinate. For this reason, it requires several slab indices and weights.

`ClimaCore.Remapping.interpolate_slab_level!`

— Method```
interpolate_slab_level(
field::Fields.Field,
h::Integer,
Is::Tuple,
zpts;
fill_value = eltype(field)(NaN)
)
```

Vertically interpolate the given `field`

on `zpts`

.

`interpolate_slab_level!`

interpolates several values at a fixed horizontal coordinate.

The field is linearly interpolated across two neighboring vertical elements.

For centered-valued fields, if `zcoord`

is in the top (bottom) half of a top (bottom) element in a column, no interpolation is performed and the value at the cell center is returned. Effectively, this means that the interpolation is first-order accurate across the column, but zeroth-order accurate close to the boundaries.

Return `fill_value`

when the vertical coordinate is negative.

`ClimaCore.Remapping.interpolation_weights`

— Function`interpolation_weights(horz_mesh, hcoord, quad_points)`

Return the weights (tuple of arrays) to interpolate fields onto `hcoord`

on the given mesh and quadrature points.

`ClimaCore.Remapping.target_hcoords_pid_bitmask`

— Method`target_hcoords_pid_bitmask(target_hcoords, topology, pid)`

Return a bitmask for points in `target_hcoords`

in the given (horizontal) `topology`

that belong to the process with process number `pid`

.

This mask can be used to extract the `target_hcoords`

relevant to the given `pid`

. The mask is the same shape and size as the input `target_hcoords`

, which makes it particularly useful.

`ClimaCore.Remapping.vertical_bounding_indices`

— Method`vertical_bounding_indices(space, zcoords)`

Return the vertical element indices needed to perform linear interpolation of `zcoords`

.

For centered-valued fields, if `zcoord`

is in the top (bottom) half of a top (bottom) element in a column, no interpolation is performed and the value at the cell center is returned. Effectively, this means that the interpolation is first-order accurate across the column, but zeroth-order accurate close to the boundaries.

`ClimaCore.Remapping.vertical_indices`

— Method`vertical_indices(space, zcoords)`

Return the vertical index of the element that contains `zcoords`

.

`zcoords`

is interpreted as "reference z coordinates".

`ClimaCore.Remapping.vertical_indices_ref_coordinate`

— Function`vertical_indices_ref_coordinate(space, zcoord)`

Return the vertical indices of the elements below and above `zcoord`

.

Return also the correct reference coordinate `zcoord`

for vertical interpolation.

`ClimaCore.Remapping.vertical_interpolation_weights`

— Method`vertical_interpolation_weights(space, zcoords)`

Compute the interpolation weights to vertically interpolate the `zcoords`

in the given `space`

.

This assumes a linear interpolation, where the first weight is to be multiplied with the "lower" element, and the second weight with the "higher" element in a stack.

That is, this function returns `A`

, `B`

such that `f(zcoord) = A f_lo + B f_hi`

, where `f_lo`

and `f_hi`

are the values on the neighboring elements of `zcoord`

.

`ClimaCore.Remapping.vertical_reference_coordinates`

— Method`vertical_reference_coordinates(space, zcoords)`

Return the reference coordinates of the element that contains `zcoords`

.

Reference coordinates (ξ) typically go from -1 to 1, but for center spaces this function remaps in such a way that they can directly be used for linear interpolation (so, if ξ is negative, it is remapped to (0,1), if positive to (-1, 0)). This is best used alongside with `vertical_bounding_indices`

.

`zcoords`

is interpreted as "reference z coordinates".