LinearAdaption(surface::Field)

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

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

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 diffusesurfaceelevation!(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⁻¹).

physical_z_to_ref_z(adaption::HypsographyAdaption, z_ref::ZPoint, z_surface::ZPoint, z_top::ZPoint) :: ZPoint

Convert physical zs to reference zs as prescribed by the given adaption.

This function has to be the inverse of ref_z_to_physical_z.

ref_z_to_physical_z(adaption::HypsographyAdaption, z_ref::ZPoint, z_surface::ZPoint, z_top::ZPoint) :: ZPoint

Convert reference zs to physical zs as prescribed by the given adaption.

This function has to be the inverse of physical_z_to_ref_z.

AbstractLimiter

Supertype for all limiters.

Interfaces

• apply_limiter!

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)

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.

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.

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

1. compute_element_bounds!
2. starts the ghost exchange (if distributed)
3. compute_neighbor_bounds_local!
4. completes the ghost exchange (if distributed)
5. compute_neighbor_bounds_ghost! (if distributed)

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

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

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

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

CellCenter()

Cell center location

CellFace()

Cell face location

ColumnGrid(
    full_grid :: ExtrudedFiniteDifferenceGrid, 
    colidx    :: ColumnIndex,
)

A view into a column of a ExtrudedFiniteDifferenceGrid. This can be used as an

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]

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

Construct an ExtrudedFiniteDifferenceGrid from the horizontal and vertical spaces.

FiniteDifferenceGrid(topology::Topologies.IntervalTopology)
FiniteDifferenceGrid(device::ClimaComms.AbstractDevice, 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.

Flat()

No surface hypsography.

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

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

SpectralElementSpace2D <: AbstractSpace

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

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.

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.

Grids.local_geometry_type(::Type)

Get the LocalGeometry type.

Grids.topology(grid::AbstractGrid)

Get the topology of a grid.

RecursiveApply

This module contains operators to recurse over nested Tuples or NamedTuples.

To extend to another type T, define RecursiveApply.rmap(fn, args::T...)

radd(X, Y)
X ⊞ Y

Recursively add elements of X and Y.

rconvert(T, X)

Identical to convert(T, X), but with improved type stability for nested types.

rdiv(X, Y)

Recursively divide each element of X by Y

rmap(fn, X...)

Recursively apply fn to each element of X

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.

rmatmul1(W, S, i, j)

Recursive matrix product along the 1st dimension of S. Equivalent to:

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

rmatmul2(W, S, i, j)

Recursive matrix product along the 2nd dimension S. Equivalent to:

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

rmul(X, Y)
X ⊠ Y

Recursively scale each element of X by Y.

rmuladd(w, X, Y)

Recursively add elements of w * X + Y.

rpromote_type(Ts...)

Recursively apply promote_type to the input types.

rsub(X, Y)
X ⊟ Y

Recursively subtract elements of Y from X.

rzero(T)

Recursively compute the zero value of type T.

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:

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

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

Utilities.Cache.cached_objects()

List all currently cached objects.

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.

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.

AbstractAxis

An axis of a AxisTensor.

AbstractGlobalGeometry

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

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)

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.

CartesianGlobalGeometry()

Specifies that the local coordinates align with the Cartesian coordinates, e.g. XYZPoint aligns with Cartesian123Point, and UVWVector aligns with Cartesian123Vector.

DeepSphericalGlobalGeometry(radius)

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

LocalGeometry

The necessary local metric information defined at each node.

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.

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.

SurfaceGeometry

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

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

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.

blockmat(m11, m22[, m12])

Construct an Axis2Tensor from sub-blocks

components(a::AxisTensor)

Returns a StaticArray containing the components of a in its stored basis.

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.

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.

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

euclidean_distance(pt1::XYPoint, pt2::XYPoint)

Compute the 2D or 3D Euclidean distance between pt1 and pt2.

float_type(T)

Return the floating point type backing T: T can either be an object or a type.

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.

great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint, global_geometry::SphericalGlobalGeometry)

Compute the great circle (spherical geodesic) distance between pt1 and pt2.

great_circle_distance(pt1::LatLongZPoint, pt2::LatLongZPoint, global_geometry::SphericalGlobalGeometry)

Compute the great circle (spherical geodesic) distance between pt1 and pt2.

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

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.

mul_with_projection(x, y, lg)

Similar to x * y, except that this version automatically projects y to avoid DimensionMismatch errors for AxisTensors. 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.

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

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.

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.

rmul_with_projection(x, y, lg)

Similar to rmul(x, y), except that this version calls mul_with_projection instead of *.

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

@pointtype name fieldname1 ...

Define a subtype name of AbstractPoint with appropriate conversion functions.

MatrixFields

This module adds support for defining and manipulating Fields that represent matrices. Specifically, it adds the BandMatrixRow type, which can be used to store the entries of a band matrix. A Field of BandMatrixRows on a FiniteDifferenceSpace can be interpreted as a band matrix by vertically concatenating the BandMatrixRows. Similarly, a Field of BandMatrixRows on an ExtrudedFiniteDifferenceSpace can be interpreted as a collection of band matrices, one for each column of the Field. Such Fields are called ColumnwiseBandMatrixFields, 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 Tuples or NamedTuples 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 Fields. Specifically, it adds the FieldMatrix type, which is a dictionary that maps pairs of FieldNames to ColumnwiseBandMatrixFields 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 Fields of Tuples or NamedTuples, or as separate matrix Fields of single values
• The ability to solve linear equations using FieldMatrixSolver, which is a generalization of ldiv! that is designed to optimize solver performance

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.

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)

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 FieldNames in names₁ correspond to the subscript ₁, while all other FieldNames 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.

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.

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 FieldNames in names₁ correspond to the subscript ₁, while all other FieldNames correspond to the subscript ₂. This algorithm has only 1 step:

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

BlockDiagonalPreconditioner()

A PreconditionerAlgorithm that sets P to the block diagonal entries of A. This is also called a "block Jacobi" preconditioner.

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

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 FieldNames in names₁ correspond to the subscript ₁, while all other FieldNames correspond to the subscript ₂. This algorithm has 2 steps:

1. Solve A₁₁ * x₁ = b₁ for x₁ using the algorithm alg₁, which is set to a BlockDiagonalSolve by default.
2. Solve A₂₂ * x₂ = b₂ - A₂₁ * x₁ for x₂ using the algorithm alg₂, which is also set to a BlockDiagonalSolve by default.

CustomPreconditioner(M, alg)

A PreconditionerAlgorithm that sets P to the FieldMatrix M and inverts P using the FieldMatrixSolverAlgorithm alg.

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.

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

• field_matrix_solver_cache
• check_field_matrix_solver
• run_field_matrix_solver!

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

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

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 ColumnwiseBandMatrixFields or multiples of LinearAlgebra.I; this is the only user-facing subtype of FieldNameDict
• FieldVectorView, which maps a set of FieldVectorKeys to Fields; 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 Fields when they are materialized. Many internal operations generate lazy FieldNameDicts 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 FieldNameDicts, 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 DiagonalMatrixRows or a multiple of LinearAlgebra.I

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 FieldNames 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 FieldNameSets do not have any performance cost at runtime (as long as their arguments are inferrable).

Unlike other AbstractSets, FieldNameSet has special behavior for overlapping values. For example, the FieldNames @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_trees must be identical.

FieldNameTree(x)

Tree of FieldNames 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).

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.

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.

MainDiagonalPreconditioner()

A PreconditionerAlgorithm that sets P to the main diagonal of A. This is also called a "Jacobi" preconditioner.

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 Fields. For Fields 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 ColumnwiseBandMatrixFields, 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.\]

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:

• solver_algorithm
• lazy_preconditioner

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 FieldNames in names₁ correspond to the subscript ₁, while all other FieldNames correspond to the subscript ₂. This algorithm has 3 steps:

1. Solve A₁₁ * x₁′ = b₁ for x₁′ using the algorithm alg₁, which is set to a BlockDiagonalSolve by default.
2. Solve (A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂ = b₂ - A₂₁ * x₁′ for x₂ using the algorithm alg₂.
3. 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.

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.

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

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.

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.

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.

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.

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.

column_axes(matrix_field, [matrix_space])

Returns the space that corresponds to the columns of matrix_field, i.e., the axes of the Fields 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).

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 BandMatrixRows. This involves copying the data stored in the field. Because BandedMatrix does not currently support operations with CuArrays, all GPU data is copied to the CPU.

column_field2array_view(field)

Similar to column_field2array(field), except that this version avoids copying the data stored in the field.

concrete_field_vector(vector)

Converts the FieldVectorView vector back into a FieldVector.

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.

field2arrays_view(field)

Similar to field2arrays(field), except that this version calls column_field2array_view instead of column_field2array.

field_matrix_solve!(solver, x, A, b)

Solves the equation A * x = b for x using the FieldMatrixSolver solver.

field_matrix_solver_cache(alg, A, b)

Allocates the cache required by the FieldMatrixSolverAlgorithm alg to solve the equation A * x = b.

field_vector_view(x, [name_tree])

Constructs a FieldVectorView that contains all of the Fields in the FieldVector x. The default name_tree is FieldNameTree(x), but this can be modified if needed.

is_lazy(dict)

Checks whether the FieldNameDict dict contains any un-materialized AbstractBroadcasted entries.

lazy_main_diagonal(matrix)

Constructs a lazy FieldMatrix that contains the main diagonal of matrix.

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.

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.

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

matrix_shape(matrix_field, [matrix_space])

Returns either Square(), FaceToCenter(), or CenterToFace(), depending on whether the diagonal indices of matrix_field are Ints or PlusHalfs and whether matrix_space is on cell centers or cell faces. By default, matrix_space is set to axes(matrix_field).

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

outer_diagonals(::Type{<:BandMatrixRow})

Gets the indices of the lower and upper diagonals, ld and ud, of the given subtype of BandMatrixRow.

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.

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.

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.

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.

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

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]

Meshes

• domain
• topology
• coordinates
• metric terms (inverse partial derivatives)
• quadrature rules and weights

References / notes

• ceed
• QA

AbstractSpace

Should define

• grid
• staggering

• space constructor

FiniteDifferenceSpace(
    grid::Grids.FiniteDifferenceGrid,
    staggering::Staggering, 
)

PointSpace <: AbstractSpace

A zero-dimensional space.

SpectralElementSpace1D(grid::SpectralElementGrid1D)

SpectralElementSpace2D(grid::SpectralElementGrid1D)

SpectralElementSpaceSlab <: AbstractSpace

A view into a SpectralElementSpace2D for a single slab.

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.

Spaces.local_area(space::Spaces.AbstractSpace)

The length/area/volume of space local to the current context. See Spaces.area

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.

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.

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!

2). Spaces.weighted_dss_internal!

3). Spaces.weighted_dss_ghost!

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.

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.

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

Δz_data(space::AbstractSpace)

A DataLayout containing the Δz on a given space space.

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

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

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

Creates a DSSBuffer for the field data corresponding to data

ClosedUniform{Nq}()

Uniformly-spaced quadrature including boundary.

GL{Nq}()

Gauss-Legendre quadrature using Nq quadrature points.

GLL{Nq}()

Gauss-Legendre-Lobatto quadrature using Nq quadrature points.

QuadratureStyle

Quadrature style supertype. See sub-types:

• GLL
• GL
• Uniform

Uniform{Nq}()

Uniformly-spaced quadrature.

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.

degrees_of_freedom(QuadratureStyle) -> Int

Returns the degreesoffreedom of the QuadratureStyle concrete type

differentiation_matrix(FT, quadstyle::QuadratureStyle)

The spectral differentiation matrix at the quadrature points of quadstyle, using floating point types FT.

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.

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.

V = orthonormal_poly(points, quad)

V_{ij} contains the j-1th Legendre polynomial evaluated at points[i]. i.e. it is the mapping from the modal to the nodal representation.

polynomial_degree(QuadratureStyle) -> Int

Returns the polynomial degree of the QuadratureStyle concrete type

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

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

AbstractDomain

A domain represents a region of space.

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.

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.

SphereDomain(radius)

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

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

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

AbstractFieldStyle

The supertype of all broadcasting-like operations on Fields.

Field(values, space)

A set of values defined at each point of a space.

FieldStyle{DS <: DataStyle}

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

FieldVector

A FieldVector is a wrapper around one or more Fields 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, ....

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.

fill!(field::Field, value)

Fill field with value.

fill(value, space::AbstractSpace)

Create a new Field on space and fill it with value.

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.

ones(space::AbstractSpace)

Construct a field on space that is one everywhere.

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.

zeros(space::AbstractSpace)

Construct a field on space that is zero everywhere.

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.

backing_array(x)

The AbstractArray that is backs an object x, allowing it to be treated as a component of a FieldVector.

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.

Example

∇ = GradientF2C()
div = DivergenceC2F()

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

coordinate_field(space::AbstractSpace)

Construct a Field of the coordinates of the space.

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.

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

local_geometry_field(space::AbstractSpace)

Construct a Field of the LocalGeometry of the space.

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.

ncolumns(::Field)
ncolumns(::Space)

Number of columns in a given space.

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.

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.

Fields.unwrap(x::T)

This is called when calling getproperty on a FieldVector property of element type T.

Fields.wrap(x)

Construct a mutable wrapper around x. This can be extended for new types (especially immutable ones).

Δz_field(field::Field)
Δz_field(space::AbstractSpace)

A Field containing the Δz values on the same space as the given field.

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

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

Call Spaces.weighted_dss! on multiple fields at once, overlapping communication as much as possible.

Spaces.create_dss_buffer(field::Field)

Create a buffer for communicating neighbour information of field.

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

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.

AbstractCubedSphere <: AbstractMesh2D

This is an abstract type of cubed-sphere meshes on SphereDomains. 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

• S2Geometry library
• MIT GCM exch2

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:

• SharedVertices(mesh, elem, vert)
• face_connectivity_matrix(mesh[,elemorder])
• vertex_connectivity_matrix(mesh[,elemorder])

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.

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.

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.

ExponentialStretching(H::FT)

Apply exponential stretching to the domain when constructing elements. H is the scale height (a typical atmospheric scale height H ≈ 7.5km).

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)

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)

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:

• Uniform()
• ExponentialStretching(H)
• GeneralizedExponentialStretching(dz_bottom, dz_top)

IntrinsicMap()

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

LocalElementMap

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

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.

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.

Meshes.SharedVertices(mesh, elem, vert)

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

Uniform()

Use uniformly-sized elements.

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.

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.

panel = cubedspherepanel(coord::Geometry.Cartesian123Point)

Given a point coord, return its panel number (an integer between 1 and 6).

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.

Meshes.domain(mesh::AbstractMesh)

The domain (a subtype of Domains.AbstractDomain) on which the mesh is defined.

Meshes.element_horizontal_length_scale(mesh::AbstractMesh)

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

Meshes.elements(mesh::AbstractMesh)

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

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.

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

Given a point coord and its panel number panel, return its coordinates in panel 1 (coord1).

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

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.

nelements(mesh::AbstractMesh)

The number of elements in the mesh.

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

The element and face (opelem, opface) that oppose face face of element elem.

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

i = refindex(ϕ, n)

Given a reference coordinate ϕ in the interval [-1, 1], divide the interval intonevenly spaced subintervals, and return the index of the subinterval (i), and the position in the subinterval (ϕs), normalized to normalized[-1, 1]`.

to_panel(panel, coord1::Geometry.Cartesian123Point)

Given a point at coord1 on panel 1 of a sphere, transform it to panel panel.

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.

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.

AbstractTopology

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

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

IntervalTopology([context::SingletonCommsContext,] mesh::IntervalMesh)

A sequential topology on an Meshes.IntervalMesh.

Perimeter2D(Nq)

Construct a perimeter iterator for a 2D spectral element with Nq nodes per dimension (i.e. polynomial degree Nq-1).

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

Topologies.VertexIterator

An iterator over the unique (shared) vertices of the topology topology. Each vertex returns a Vertex object, which is itself an iterator.

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.

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.

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.

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

dss!(data, topology)

Computed unweighted/pure DSS of data.

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

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

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

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

Apply dss to local vertices.

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

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

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.

See ClimaCore.Spaces.weighted_dss!.

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

dss_untransform(T, targ, local_geometry, I...)

Transform targ[I...] back to a value of type T after performing direct stiffness summation (DSS).

See ClimaCore.Spaces.weighted_dss!.

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

The node indices of the qth node on face face, where Nq is the number of face nodes in each direction.

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

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.

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.

ghost_vertices(topology)

An iterator over the ghost vertices of topology. Each vertex is an iterator over (isghost, lidx/ridx, vert) pairs.

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.

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

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

local_vertices(topology)

An iterator over the interior vertices of topology. Each vertex is an iterator over (lidx, vert) pairs.

localelemindex(topology, elem)

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

mesh(topology)

Returns the mesh underlying the topology

neighbors(topology)

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

nelems(topology)

The total number of elements in topology.

nghostelems(topology)

The number of ghost elements in topology.

nlocalelems(topology)

The number of local elements in topology.

nneighbors(topology)

The number of neighbors of this process in topology.

nsendelems(topology)

The number of elements to send to neighbors in topology.

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

spacefillingcurve(mesh::Meshes.AbstractCubedSphere)

Generate element ordering, elemorder, based on a space filling curve for a CubedSphere mesh.

spacefillingcurve(mesh::Meshes.RectilinearMesh)

Generate element ordering, elemorder, based on a space filling curve for a Rectilinear mesh.

(c1,c2,c3,c4) = vertex_coordinates(topology, elem)

The coordinates of the 4 vertices of element elem.

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.

const half = PlusHalf(0)

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.

cart_ind(n::NTuple, i::Integer)

Returns a CartesianIndex from the list CartesianIndices(map(x->Base.OneTo(x), n))[i] given size n and location i.

linear_ind(n::NTuple, ci::CartesianIndex)
linear_ind(n::NTuple, t::NTuple)

Returns a linear index from the list LinearIndices(map(x->Base.OneTo(x), n))[ci] given size n and cartesian index ci.

The linear_ind(n::NTuple, t::NTuple) wraps t in a Cartesian index and calls linear_ind(n::NTuple, ci::CartesianIndex).

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 Fields will be distributed using this context. As with HDF5Writer, this requires a HDF5 library with MPI support.

Interface

• read_domain
• read_mesh
• read_topology
• read_space
• read_field

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.

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.

Interface

write!

Usage

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

defaultname(obj)

Default name of object for InputOutput writers.

matrix_to_cartesianindices(elemorder_matrix)

Converts the elemorder_matrix to cartesian indices.

read_attributes(reader::AbstractReader, name::AbstractString, data::Dict)

Return the attributes associated to the object at name in the given HDF5 file.

read_domain(reader::AbstractReader, name)

Reads a domain named name from reader. Domain objects are cached in the reader to avoid creating duplicate objects.

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.

read_grid(reader::AbstractReader, name)

Reads a space named name from reader, or from the reader cache if it has already been read.

read_mesh(reader::AbstractReader, name)

Reads a mesh named name from reader, or from the reader cache if it has already been read.

read_space(reader::AbstractReader, name)

Reads a space named name from reader, or from the reader cache if it has already been read.

read_topology(reader::AbstractReader, name)

Reads a topology named name from reader, or from the reader cache if it has already been read.

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 Fields and FieldVectors are not cached, and so can be written multiple times.

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

Write one or more name => value pairs to the HDF5 file filename.

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

Write one or more name => value pairs to writer.

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

Write data as attributes to the object at name in the given HDF5 file.

write_new!(writer, domain, name)

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

write_new!(writer, topology, name)

Write IntervalTopology data to HDF5.

write_new!(writer, domain, name)

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

write_new!(writer, domain, name)

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

write_new!(writer, mesh, name)

Write CubedSphereMesh data to HDF5.

write_new!(writer, topology, name)

Write Topology2D data to HDF5.

write_new!(writer, space, name)

Write SpectralElementSpace1D data to HDF5.

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.

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.

column(data::AbstractData, i::Integer)

A contiguous "column" view into an underlying data layout data at nodal point index i.

slab(data::AbstractData, h::Integer)

A "pancake" view into an underlying data layout data at location h.

AbstractBoundaryCondition

An abstract type for boundary conditions for FiniteDifferenceOperators.

Subtypes should define:

• boundary_width
• stencil_left_boundary
• stencil_right_boundary

AbstractOperator

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

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

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.

CentralNumericalFlux(fluxfn)

Evaluates the central numerical flux using fluxfn.

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

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

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

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

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

Extrapolate()

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

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.

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.

FiniteDifferenceOperator

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

• return_eltype
• return_space
• stencil_interior_width
• stencil_interior

See also AbstractBoundaryCondition for how to define the boundaries.

FirstOrderOneSided()

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

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

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

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

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.

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

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

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

Supported boundary conditions are:

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

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

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

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

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

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

InterpolateF2C()

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

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

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

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

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

LinearRemap(target::AbstractSpace, source::AbstractSpace)

A remapping operator from the source space to the target space.

See Ullrich2015 eqs. 3 and 4.

NullBoundaryCondition()

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

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

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

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

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

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

RusanovNumericalFlux(fluxfn, wavespeedfn)

Evaluates the Rusanov numerical flux using fluxfn with wavespeed wavespeedfn

SetBoundaryOperator(;boundaries...)

This operator only modifies the values at the boundary:

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

SetCurl(val)

Set the curl at the boundary to be val.

SetDivergence(val)

Set the divergence at the boundary to be val.

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.

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.

SlabBlockSpectralStyle()

Applies spectral-element operations using by making use of intermediate temporaries for each operator. This is used for CPU kernels.

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

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

SpectralStyle()

Broadcasting requires use of spectral-element operations.

StencilBroadcasted{Style}(op, args[,axes[, work]])

This is similar to a Base.Broadcast.Broadcasted object.

This is returned by Base.Broadcast.broadcasted(op::FiniteDifferenceOperator).

ThirdOrderOneSided()

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

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.

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

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

wdiv = WeakDivergence()
wdiv.(u)

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

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

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

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

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

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

It can be written in matrix form as

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

which reduces to

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

where

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

wgrad = WeakGradient()
wgrad.(f)

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

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

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

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

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

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

In matrix form, this becomes

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

which reduces to

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

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

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

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

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

Supported boundary conditions are:

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

These have the same stencil as in InterpolateC2F.

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

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

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

No boundary conditions are required (or supported)

_resolve_operator(slabidx, args...)

Calls resolve_operator(arg, slabidx) for each arg in args

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:

• CentralNumericalFlux
• RusanovNumericalFlux

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.

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.

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 (ϕ⁻, ϕ⁺).

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.

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

copyto_slab!(out, bc, slabidx)

Copy the slab indexed by slabidx from bc to out.

is_valid_index(space, ij, slabidx)::Bool

Returns true if the node indices ij and slab indices slabidx are valid for space.

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

left_interior_window_idx(arg, space, loc)

Compute the index of the leftmost point which uses only the interior stencil of the space.

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.

local_weights(space::AbstractSpace)

Each degree of freedom is associated with a local weight Ji. For finite volumes the local weight Ji 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.

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.

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.

operator_axes(space)

Return a tuple of the axis indicies a given field operator works over.

overlap(target, source)

Computes local weights of the overlap mesh for source to target spaces.

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.

reconstruct_placeholder_broadcasted(space, obj)

Recurively reconstructs objects that have been stripped via strip_space.

remap!(target_field::Field, R::LinearRemap, source_field::Field)

Applies the remapping operator R to source_field and stores the solution in target_field.

remap(R::LinearRemap, source_field::Field)

Applies R to source_field and outputs a new field on the target space.

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

return_eltype(::Op, fields...)

Defines the element type of the result of operator Op

return_space(::Op, spaces...)

Defines the space upon which the operator Op returns given arguments on input spaces.

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

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.

stencil_interior_width(::Op, args...)

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

Example

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

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

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

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

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

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

tensor_product!(out, in, M)
tensor_product!(inout, M)

Computes the tensor product out = (M ⊗ M) * in on each element.

unionall_type(::Type{T})

Extract the type of the input, and strip it of any type parameters.

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.

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

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]

Data0D{S}

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.

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.

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.

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.

Data3D{S,Nij,Nk}

Abstract type for data storage for a 3D field made up of Nij × Nij × Nk values of type S.

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.

DataF{S, A} <: Data0D{S}

Backing DataLayout for 0D point data.

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.

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.

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

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.

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.

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

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.

IJKFVH{S, Nij, Nk}(array::AbstractArray{T, 6}) <: Data3D{S, Nij, Nk}

A 3D DataLayout. TODO: Add more docs

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.

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

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

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

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

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.

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.

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.

device_dispatch(data::AbstractData)

Returns an ToCPU or a ToCUDA for CPU and CUDA-backed arrays accordingly.

fieldtypeoffset(T,S,i)

Similar to fieldoffset(S,i), but gives result in multiples of sizeof(T) instead of bytes.

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.

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.

parent_array_type(::Type{<:AbstractArray})

Returns the parent array type underlying any wrapper types, with all dimensionality information removed.

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.

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

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.

typesize(T,S)

Similar to sizeof(S), but gives the result in multiples of sizeof(T).

(Ni, Nj, Nf, Nv, Nh) = universal_size(data::AbstractData)

Returns dimensions in a universal format for all data layouts:

• Ni number of spectral element nodal degrees of freedom in first horizontal direction
• Nj number of spectral element nodal degrees of freedom in second horizontal direction
• Nf number of field components
• Nv number of vertical degrees of freedom
• Nh number of horizontal elements

Note: this is similar to Base.size, except that universal_size does not return 1 for the number of field components.

Remapper(space, targethcoords, targetzcoords; bufferlength = 1) Remapper(space, targethcoords; 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.

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

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

_reset_interpolated_values!(remapper::Remapper)

Reset the local (private) state in remapper. This function has to be called before performing interpolation.

_set_interpolated_values!(remapper, field)

Change the local state of remapper by performing interpolation of fields on the vertical and horizontal points.

batched_ranges(num_fields, buffer_length)

Partition the indices from 1 to numfields in such a way that no range is larger than bufferlength.

containing_pid(target_point, topology)

Return the process id that contains the target_point in the given topology.

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

interpolate(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 destiniation. 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 destination 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, :, :, :]

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)

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. NaNs are returned for values that are below the surface.

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.

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.

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.

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.

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.

vertical_indices(space, zcoords)

Return the vertical index of the element that contains zcoords.

zcoords is interpreted as "reference z coordinates".

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.

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.

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