ClimaCore.Hypsography.SLEVEAdaptionType
SLEVEAdaption(surface::Field, ηₕ::FT, s::FT)

Locate vertical levels using an exponential function between the surface field and the top of the domain, using the method of Schar2002. This method is modified such no warping is applied above some user defined parameter 0 ≤ ηₕ < 1.0, where the lower and upper bounds represent the domain bottom and top respectively. s governs the decay rate. If the decay-scale is poorly specified (i.e., s * zₜ is lower than the maximum surface elevation), a warning is thrown and s is adjusted such that it szₜ > maximum(z_surface).

ClimaCore.Hypsography.diffuse_surface_elevation!Method
diffuse_surface_elevation!(f::Field; κ::T, iter::Int, dt::T)

Option for 2nd order diffusive smoothing of generated terrain. Mutate (smooth) a given elevation profile f before assigning the surface elevation to the HypsographyAdaption type. A spectral second-order diffusion operator is applied with forward-Euler updates to generate profiles for each new iteration. Steps to generate smoothed terrain ( represented as a ClimaCore Field) are as follows:

  • Compute discrete elevation profile f
  • Compute 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⁻¹).

ClimaCore.Hypsography.physical_z_to_ref_zMethod
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.

ClimaCore.Hypsography.ref_z_to_physical_zMethod
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.

ClimaCore.Limiters.QuasiMonotoneLimiterType
QuasiMonotoneLimiter

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

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

  • ρq: tracer density Field, where q denotes tracer concentration per unit mass. This can be a scalar field, or a struct-valued field.
  • ρ: fluid density Field (scalar).

Constructor

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

Creates a limiter instance for the field ρq with relative tolerance rtol.

Usage

Call compute_bounds! on the input fields:

compute_bounds!(limiter, ρq, ρ)

Then call apply_limiter! on the output fields:

apply_limiter!(ρq, ρ, limiter)
ClimaCore.Limiters.apply_limit_slab!Method
apply_limit_slab!(slab_ρq, slab_ρ, slab_WJ, slab_q_bounds, rtol)

Apply the computed bounds of the tracer concentration (slab_q_bounds) in the limiter to slab_ρq, given the total mass slab_ρ, metric terms slab_WJ, and relative tolerance rtol. Return whether the tolerance condition could be satisfied.

ClimaCore.Limiters.apply_limiter!Method
apply_limiter!(ρq, ρ, limiter::QuasiMonotoneLimiter)

Apply the limiter on the tracer density ρq, using the computed desired bounds on the concentration q and density ρ as an optimal weight. This iterates over each element, calling apply_limit_slab!. If the limiter fails to converge for any element, a warning is issued.

ClimaCore.Grids.ColumnGridType
ColumnGrid(
    full_grid :: ExtrudedFiniteDifferenceGrid, 
    colidx    :: ColumnIndex,
)

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

ClimaCore.Grids.ColumnIndexType
ColumnIndex(ij,h)

An index into a column of a field. This can be used as an argument to getindex of a Field, to return a field on that column.

Example

colidx = ColumnIndex((1,1),1)
field[colidx]
ClimaCore.Grids.ExtrudedFiniteDifferenceGridType
ExtrudedFiniteDifferenceGrid(
    horizontal_space::AbstractSpace,
    vertical_space::FiniteDifferenceSpace,
    hypsography::HypsographyAdaption = Flat(),
)

Construct an ExtrudedFiniteDifferenceGrid from the horizontal and vertical spaces.

ClimaCore.Grids.FiniteDifferenceGridType
FiniteDifferenceGrid(topology::Topologies.IntervalTopology)
FiniteDifferenceGrid(mesh::Meshes.IntervalMesh)

Construct a FiniteDifferenceGrid from an IntervalTopology (or an IntervalMesh).

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

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

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

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

ClimaCore.Grids.SpectralElementGrid2DMethod
SpectralElementSpace2D(topology, quadrature_style; enable_bubble)

Construct a SpectralElementSpace2D instance given a topology and quadrature. The flag enable_bubble enables the bubble correction for more accurate element areas.

Input arguments:

  • topology: Topology2D
  • quadrature_style: QuadratureStyle
  • enable_bubble: Bool

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

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

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

and the case of non linear elements, by

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

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

Note: This is accurate only for cubed-spheres of the Meshes.EquiangularCubedSphere and Meshes.EquidistantCubedSphere type, not for Meshes.ConformalCubedSphere.

ClimaCore.Grids.local_geometry_dataFunction
Grids.local_geometry_data(
    grid       :: AbstractGrid,
    staggering :: Union{Staggering, Nothing},
)

Get the DataLayout object containing the local geometry data information of the grid with staggering staggering.

If the grid is not staggered, staggering should be nothing.

ClimaCore.RecursiveApplyModule
RecursiveApply

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

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

ClimaCore.RecursiveApply.rmaptypeMethod
rmaptype(fn, T)
rmaptype(fn, T1, T2)

Recursively apply fn to each type parameter of the type T, or to each type parameter of the types T1 and T2, where fn returns a type.

ClimaCore.RecursiveApply.rmatmul1Method
rmatmul1(W, S, i, j)

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

mapreduce(⊠, ⊞, W[i,:], S[:,j])
ClimaCore.RecursiveApply.rmatmul2Method
rmatmul2(W, S, i, j)

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

mapreduce(⊠, ⊞, W[j,:], S[i, :])
ClimaCore.Geometry.AbstractPointType
AbstractPoint

Represents a point in space.

The following types are supported:

  • XPoint(x)
  • YPoint(y)
  • ZPoint(z)
  • XYPoint(x, y)
  • XZPoint(x, z)
  • XYZPoint(x, y, z)
  • LatPoint(lat)
  • LongPoint(long)
  • LatLongPoint(lat, long)
  • LatLongZPoint(lat, long, z)
  • Cartesian1Point(x1)
  • Cartesian2Point(x2)
  • Cartesian3Point(x3)
  • Cartesian12Point(x1, x2)
  • Cartesian13Point(x1, x3)
  • Cartesian123Point(x1, x2, x3)
ClimaCore.Geometry.AxisTensorType
AxisTensor(axes, components)

An AxisTensor is a wrapper around a StaticArray, where each dimension is labelled with an AbstractAxis. These axes must be consistent for operations such as addition or subtraction, or be dual to each other for operations such as multiplication.

See also

components to obtain the underlying array.

ClimaCore.Geometry.CartesianGlobalGeometryType
CartesianGlobalGeometry()

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

ClimaCore.Geometry.SphericalGlobalGeometryType
SphericalGlobalGeometry(radius)

Specifies that the local coordinates are specified in reference to a sphere of radius radius. The x1 axis is aligned with the zero longitude line.

The local vector basis has u in the zonal direction (with east being positive), v in the meridonal (north positive), and w in the radial direction (outward positive). For a point located at the pole, we take the limit along the zero longitude line:

  • at the north pole, this corresponds to u being aligned with the x2 direction, v being aligned with the negative x1 direction, and w being aligned with the x3 direction.
  • at the south pole, this corresponds to u being aligned with the x2 direction, v being aligned with the x1 direction, and w being aligned with the negative x3 direction.
ClimaCore.Geometry.bilinear_interpolateMethod
bilinear_interpolate(coords::NTuple{4}, ξ1, ξ2)

Bilinear interpolate between coords by parameters ξ1, ξ2, each in the interval [-1,1].

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

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

        ξ1-->
ClimaCore.Geometry.curl_result_typeMethod
curl_result_type(Val(I), Val(L), V)

The return type when taking the curl along dimensions I of a field of eltype V, defined on dimensions L

Required for statically infering the result type of the curl operation for AxisVector subtypes. Curl is only defined for CovariantVector` field input types.

Input VectorOperator directionCurl output vector
Covariant12Vector(1,2)Contravariant3Vector
Covariant3Vector(1,2)Contravariant12Vector
Covariant123Vector(1,2)Contravariant123Vector
Covariant1Vector(1,)Contravariant1Vector
Covariant2Vector(1,)Contravariant3Vector
Covariant3Vector(1,)Contravariant2Vector
Covariant12Vector(3,)Contravariant12Vector
Covariant1Vector(3,)Contravariant2Vector
Covariant2Vector(3,)Contravariant1Vector
Covariant3Vector(3,)Contravariant3Vector
ClimaCore.Geometry.divergence_result_typeMethod
divergence_result_type(V)

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

Required for statically infering the result type of the divergence operation for AxisVector subtypes.

ClimaCore.Geometry.dualFunction
dual(ax::AbstractAxis)

The dual axis to ax.

julia> using ClimaCore.Geometry

julia> Geometry.dual(Geometry.Covariant12Axis())
ClimaCore.Geometry.ContravariantAxis{(1, 2)}()

julia> Geometry.dual(Geometry.Cartesian123Axis())
ClimaCore.Geometry.CartesianAxis{(1, 2, 3)}()
ClimaCore.Geometry.gradient_result_typeMethod
gradient_result_type(Val(I), V)

The return type when taking the gradient along dimension I of a field V.

Required for statically infering the result type of the gradient operation for AxisVector subtypes.

ClimaCore.Geometry.great_circle_distanceMethod
great_circle_distance(pt1::LatLongPoint, pt2::LatLongPoint, global_geometry::SphericalGlobalGeometry)

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

ClimaCore.Geometry.great_circle_distanceMethod
great_circle_distance(pt1::LatLongZPoint, pt2::LatLongZPoint, global_geometry::SphericalGlobalGeometry)

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

ClimaCore.Geometry.linear_interpolateMethod
interpolate(coords::NTuple{2}, ξ1)

Interpolate between coords by parameters ξ1 in the interval [-1,1]. The type of interpolation is determined by the element type of coords

ClimaCore.Geometry.outerFunction
outer(x, y)
x ⊗ y

Compute the outer product of x and y. Typically x will be a vector, and y can be either a number, vector or tuple/named tuple.

# vector ⊗ scalar = vector
julia> [1.0,2.0] ⊗ 2.0
2-element Vector{Float64}:
 2.0
 4.0

# vector ⊗ vector = matrix
julia> [1.0,2.0] ⊗ [1.0,3.0]
2×2 Matrix{Float64}:
 1.0  3.0
 2.0  6.0

# vector ⊗ tuple = recursion
julia> [1.0,2.0] ⊗ (1.0, (a=2.0, b=3.0))
([1.0, 2.0], (a = [2.0, 4.0], b = [3.0, 6.0]))
ClimaCore.Geometry.projectFunction
project(axis, V[, local_geometry])

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

This is equivalent to transform, but will not throw an error if the conversion is not exact.

ClimaCore.Geometry.transformFunction
transform(axis, V[, local_geometry])

Transform the first axis of the vector or tensor V to axis. This will throw an error if the conversion is not exact.

The conversion rules are defined as:

  • v::Covariant => Local: ∂ξ∂x' * v
  • v::Local => Contravariant: ∂ξ∂x * v
  • v::Contravariant => Local: ∂x∂ξ * v
  • v::Local => Covariant: ∂x∂ξ' * v
  • v::Covariant => Contravariant: ∂ξ∂x * (∂ξ∂x' * v) = gⁱʲ * v
  • v::Contravariant => Covariant: ∂x∂ξ' * ∂x∂ξ * v = gᵢⱼ * v

Example

Consider the conversion from a Covariant12Vector to a Contravariant12Axis. Mathematically, we can write this as

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

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

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

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

[ v¹ ]   [g¹¹ v₁  + g¹² v₂ ]
[ v² ] = [g²¹ v₁  + g²² v₂ ]
[ 0  ]   [<asserts g²³ v₁  + g²³ v₂ == 0>]
ClimaCore.Geometry.@pointtypeMacro
@pointtype name fieldname1 ...

Define a subtype name of AbstractPoint with appropriate conversion functions.

ClimaCore.Utilities.CacheModule

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.

ClimaCore.Utilities.Cache.clean_cache!Method
Utilities.Cache.clean_cache!(object)

Remove object from the cache of created objects.

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

ClimaCore.Utilities.Cache.clean_cache!Method
Utilities.Cache.clean_cache!()

Remove all objects from the cache of created objects.

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

ClimaCore.MatrixFieldsModule
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
ClimaCore.MatrixFields.AbstractLazyOperatorType
AbstractLazyOperator

Supertype for "lazy operators", i.e., operators that can be called without any arguments by users, as long as they appear in broadcast expressions that contain at least one Field. If lazy_op is an AbstractLazyOperator, the expression lazy_op.() will internally be translated to non_lazy_op.(fields...), as long as it appears in a broadcast expression with at least one Field. This translation is done by the function replace_lazy_operator(space, lazy_op), which must be implemented by every subtype of AbstractLazyOperator.

ClimaCore.MatrixFields.BandMatrixRowType
BandMatrixRow{ld}(entries...)

Stores the nonzero entries in a row of a band matrix, starting with the lowest diagonal, which has index ld. Supported operations include accessing the entry on the diagonal with index d by calling row[d], taking linear combinations with other band matrix rows (and with LinearAlgebra.I), and checking for equality with other band matrix rows (and with LinearAlgebra.I). There are several aliases for commonly used subtypes of BandMatrixRow:

  • DiagonalMatrixRow(entry_1)
  • BidiagonalMatrixRow(entry_1, entry_2)
  • TridiagonalMatrixRow(entry_1, entry_2, entry_3)
  • QuaddiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4)
  • PentadiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4, entry_5)
ClimaCore.MatrixFields.BlockArrowheadPreconditionerType
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.

ClimaCore.MatrixFields.BlockArrowheadSchurComplementPreconditionerType
BlockArrowheadSchurComplementPreconditioner(; [P_alg₁], [alg₂])

A PreconditionerAlgorithm that is equivalent to a BlockArrowheadPreconditioner, but only applied to the Schur complement of A₁₁ in A, A₂₂ - A₂₁ * inv(A₁₁) * A₁₂, which is represented by a LazySchurComplement. Specifically, the preconditioner this generates is the Schur complement of P₁₁ in P, A₂₂ - A₂₁ * inv(P₁₁) * A₁₂, where P₁₁ is generated by P_alg₁. Unlike the BlockArrowheadPreconditioner constructor, this constructor does not require names₁ because the block structure of A can be inferred from the LazySchurComplement.

ClimaCore.MatrixFields.BlockArrowheadSolveType
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.

ClimaCore.MatrixFields.BlockDiagonalSolveType
BlockDiagonalSolve()

A FieldMatrixSolverAlgorithm for a block diagonal matrix:

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

This algorithm solves the N block equations Aₙₙ * xₙ = bₙ in sequence (though we might want to parallelize it in the future).

If Aₙₙ is a diagonal matrix, the equation Aₙₙ * xₙ = bₙ is solved by making a single pass over the data, setting each xₙ[i] to inv(Aₙₙ[i, i]) * bₙ[i].

Otherwise, the equation Aₙₙ * xₙ = bₙ is solved using Gaussian elimination (without pivoting), which makes two passes over the data. This is currently only implemented for tri-diagonal and penta-diagonal matrices Aₙₙ. In Gaussian elimination, Aₙₙ is effectively factorized into the product Lₙ * Dₙ * Uₙ, where Dₙ is a diagonal matrix, and where Lₙ and Uₙ are unit lower and upper triangular matrices, respectively. The first pass multiplies both sides of the equation by inv(Lₙ * Dₙ), replacing Aₙₙ with Uₙ and bₙ with Uₙxₙ, which is referred to as putting Aₙₙ into "reduced row echelon form". The second pass solves Uₙ * xₙ = Uₙxₙ for xₙ with a unit upper triangular matrix solver, which is referred to as "back substitution". These operations can become numerically unstable when Aₙₙ has entries with large disparities in magnitude, but avoiding this would require swapping the rows of Aₙₙ (i.e., replacing Dₙ with a partial pivoting matrix).

ClimaCore.MatrixFields.BlockLowerTriangularSolveType
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.
ClimaCore.MatrixFields.FieldMatrixSolverType
FieldMatrixSolver(alg, A, b)

Combination of a FieldMatrixSolverAlgorithm alg and the cache it requires to solve the equation A * x = b for x. The values of A and b that get passed to this constructor should be similar to the ones that get passed to field_matrix_solve! in order to ensure that the cache gets allocated correctly.

ClimaCore.MatrixFields.FieldMatrixSolverAlgorithmType
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:

ClimaCore.MatrixFields.FieldNameType
FieldName(name_chain...)

A singleton type that represents a chain of getproperty calls, which can be used to access a property or sub-property of an object x using the function get_field(x, name). The entire object x can also be accessed with the empty FieldName().

ClimaCore.MatrixFields.FieldNameDictType
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
ClimaCore.MatrixFields.FieldNameSetType
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.

ClimaCore.MatrixFields.FieldNameTreeType
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).

ClimaCore.MatrixFields.LazyFieldMatrixSolverAlgorithmType
LazyFieldMatrixSolverAlgorithm

A FieldMatrixSolverAlgorithm that does not require A to be a FieldMatrix, i.e., a "matrix-free" algorithm. Internally, a FieldMatrixSolverAlgorithm (for example, SchurComplementReductionSolve) might run a LazyFieldMatrixSolverAlgorithm on a "lazy" representation of a FieldMatrix (like a LazySchurComplement).

The only operations used by a LazyFieldMatrixSolverAlgorithm that depend on A are lazy_mul and, when required, lazy_preconditioner. These and other lazy operations are used to minimize the number of calls to Base.materialize!, since each call comes with a small performance penalty.

ClimaCore.MatrixFields.LazySchurComplementType
LazySchurComplement(A₁₁, A₁₂, A₂₁, A₂₂, [alg₁, cache₁, A₁₂_x₂, invA₁₁_A₁₂_x₂])

An analogue of a FieldMatrix that represents the Schur complement of A₁₁ in A, A₂₂ - A₂₁ * inv(A₁₁) * A₁₂. Since inv(A₁₁) will generally be a dense matrix, it would not be efficient to directly compute the Schur complement. So, this object only supports the "lazy" functions lazy_mul, which allows it to be multiplied by the vector x₂, and lazy_preconditioner, which allows it to be approximated with a FieldMatrix.

The values alg₁, cache₁, A₁₂_x₂, and invA₁₁_A₁₂_x₂ need to be specified in order for lazy_mul to be able to compute inv(A₁₁) * A₁₂ * x₂. When a LazySchurComplement is not passed to lazy_mul, these values can be omitted.

ClimaCore.MatrixFields.MultiplyColumnwiseBandMatrixFieldType
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.\]

ClimaCore.MatrixFields.PreconditionerAlgorithmType
PreconditionerAlgorithm

Description of how to approximate a FieldMatrix or something similar like a LazySchurComplement with a preconditioner P for which P * x = b is easy to solve for x. If P is a diagonal matrix, then x can be computed as @. inv(P) * b; otherwise, the PreconditionerAlgorithm must specify a FieldMatrixSolverAlgorithm that can be used to solve P * x = b for x.

Interface

Every subtype of PreconditionerAlgorithm must implement methods for the following functions:

ClimaCore.MatrixFields.SchurComplementReductionSolveType
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.

ClimaCore.MatrixFields.StationaryIterativeSolveType
StationaryIterativeSolve(; [kwargs...])

A LazyFieldMatrixSolverAlgorithm that solves A * x = b by setting x to some initial value x[0] (usually just the zero vector, $\mathbf{0}$) and then iteratively updating it to

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

The matrix P is called a "left preconditioner" for A. In general, this algorithm converges more quickly when P is a close approximation of A, although more complicated approximations often come with a performance penalty.

Background

Let x' denote the value of x for which A * x = b. Replacing b with A * x' in the formula for x[n] tells us that

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

In other words, the error on iteration n, x[n] - x', can be expressed in terms of the error on the previous iteration, x[n - 1] - x', as

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

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

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

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

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

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

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

The spectral radius has the property that

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

So, as the value of n increases, the norm of the error becomes bounded by

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

This indicates that x[n] will converge to x' (i.e., that the norm of the error will converge to 0) when ρ(I - inv(P) * A) < 1, and that the convergence rate is roughly bounded by ρ(I - inv(P) * A) for large values of n. More precisely, it can be shown that x[n] will converge to x' if and only if ρ(I - inv(P) * A) < 1. In practice, though, the convergence eventually stops due to the limits of floating point precision.

Also, if we assume that x[n] ≈ x', we can use the formula for x[n] to approximate the error on the previous iteration as

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

Debugging

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

  • error_norm: prints ||x[n] - x'||₂ on every iteration, approximating the error x[n] - x' as described above
  • spectral_radius: prints ρ(I - inv(P) * A), approximating this value with the eigsolve function from KrylovKit.jl

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

Keyword Arguments

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

  • P_alg = nothing: a PreconditionerAlgorithm that specifies how to compute P and solve P * x = b for x, or nothing if preconditioning is not required (in which case P is effectively set to one(A))
  • n_iters = 1: the number of iterations
  • correlated_solves = false: whether to set x[0] to a value of x that was generated during an earlier call to field_matrix_solve!, instead of setting it to $\mathbf{0}$ (it is always set to $\mathbf{0}$ on the first call to field_matrix_solve!)
  • eigsolve_kwargs = (;): keyword arguments for the eigsolve function that can be used to tune its accuracy and speed (only applicable when debugging the spectral radius)
  • debug = nothing: keyword argument for printing debug information to @debug. By default, debug is set to true if "error_norm" or "spectral_radius" is in ENV["JULIA_DEBUG"], which must be set by users.
ClimaCore.MatrixFields.WeightedPreconditionerType
WeightedPreconditioner(M, unweighted_P_alg)

A PreconditionerAlgorithm that sets P to M * P′, where M is a diagonal FieldMatrix and P′ is the preconditioner generated by unweighted_P_alg. When the entries of M are larger than 1, this is called "relaxation" or "damping"; when the entries are smaller than 1, this is called "extrapolation".

ClimaCore.MatrixFields.ApproximateBlockArrowheadIterativeSolveMethod
ApproximateBlockArrowheadIterativeSolve(names₁...; [P_alg₁], [alg₁], [alg₂], [kwargs...])

Shorthand for constructing a SchurComplementReductionSolve whose alg₂ is set to a StationaryIterativeSolve with a BlockArrowheadSchurComplementPreconditioner. The keyword argument alg₁ is passed to the constructor for SchurComplementReductionSolve, the keyword arguments P_alg₁ and alg₂ are passed to the constructor for BlockArrowheadSchurComplementPreconditioner, and all other keyword arguments are passed to the constructor for StationaryIterativeSolve.

This algorithm is somewhat similar to a StationaryIterativeSolve with a BlockArrowheadPreconditioner, but it usually converges much more quickly, i.e., the spectral radius of its iteration matrix (I - inv(P) * A) tends to be significantly smaller. Roughly speaking, this is because it runs the iterative solver on an equation with fewer variables (the Schur complement equation, (A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂ = b₂′), which means that, on each iteration, it accumulates less error due to coupling between variables. However, even though it converges more quickly, its iterations take longer because they involve using alg₁ to invert A₁₁. So, when only a few iterations are needed, a StationaryIterativeSolve with a BlockArrowheadPreconditioner might be more performant.

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

ClimaCore.MatrixFields.apply_preconditionerMethod
apply_preconditioner(P_alg, P_cache, P, lazy_b)

Constructs a lazy FieldMatrix (or a concrete one when possible) that represents the product @. inv(P) * b. Here, lazy_b is a (possibly lazy) FieldVectorView that represents b.

ClimaCore.MatrixFields.check_field_matrix_solverFunction
check_field_matrix_solver(alg, cache, A, b)

Checks that the sparsity structure of A is supported by the FieldMatrixSolverAlgorithm alg, and that A is compatible with b in the equation A * x = b.

ClimaCore.MatrixFields.check_preconditionerMethod
check_preconditioner(P_alg, P_cache, A, b)

Checks that P is compatible with b in the equation P * x = b, where P is the preconditioner generated by the PreconditionerAlgorithm P_alg for A. If P_alg requires a FieldMatrixSolverAlgorithm alg to solve the equation, this also calls check_field_matrix_solver on alg.

ClimaCore.MatrixFields.column_axesFunction
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).

ClimaCore.MatrixFields.column_field2arrayMethod
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.

ClimaCore.MatrixFields.field2arraysMethod
field2arrays(field)

Converts a field defined on a FiniteDifferenceSpace or on an ExtrudedFiniteDifferenceSpace into a collection of arrays, each of which corresponds to a column of the field. This is done by calling column_field2array on each of the field's columns.

ClimaCore.MatrixFields.field_vector_viewFunction
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.

ClimaCore.MatrixFields.lazy_mulMethod
lazy_mul(A, args...)

Constructs a lazy FieldMatrix that represents the product @. *(A, args...). This involves regular broadcasting when A is a FieldMatrix, but it has more complex behavior for other objects like the LazySchurComplement.

ClimaCore.MatrixFields.lazy_preconditionerMethod
lazy_preconditioner(P_alg, A)

Constructs a lazy FieldMatrix (or a concrete one when possible) that approximates A according to the PreconditionerAlgorithm P_alg. If P_alg is nothing instead of a PreconditionerAlgorithm, this returns one(A).

ClimaCore.MatrixFields.matrix_shapeFunction
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).

ClimaCore.MatrixFields.mul_return_typeMethod
mul_return_type(X, Y)

Computes the return type of mul_with_projection(x, y, lg), where x isa X and y isa Y. This can also be used to obtain the return type of x * y, although x * y will throw an error when projection is necessary.

Note that this is equivalent to calling the internal function _return_type: Base._return_type(mul_with_projection, Tuple{X, Y, LG}), where lg isa LG.

ClimaCore.MatrixFields.mul_with_projectionMethod
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.

ClimaCore.MatrixFields.operator_matrixMethod
operator_matrix(op)

Constructs a new operator (or operator-like object) that generates the matrix applied by op to its final argument. If op_matrix = operator_matrix(op), we can use the following identities:

  • When op takes one argument, @. op(arg) == @. op_matrix() ⋅ arg.
  • When op takes multiple arguments, @. op(args..., arg) == @. op_matrix(args...) ⋅ arg.

When op takes more than one argument, operator_matrix(op) constructs a FiniteDifferenceOperator that generates the operator matrix. When op only takes one argument, it instead constructs an AbstractLazyOperator, which is internally converted into a FiniteDifferenceOperator when used in a broadcast expression. Implementing op_matrix as a lazy operator allows us to add an argument to the expression op_matrix.(), and we then use this argument to infer the space and element type of the operator matrix.

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

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

The GradientF2C() operator applies a similar matrix, but with different entries:

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

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

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

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

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

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

In general, when op has nonzero boundary conditions that make it apply an affine transformation, operator_matrix(op) will print out a warning and zero out the boundary conditions before computing the operator matrix.

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

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

When op is such an operator, operator_matrix(op) will throw an error. In the future, we may want to modify operator_matrix(op) so that it will instead return $\textrm{op}'(\textbf{0})$, where $\textbf{0} ={}$zero.(arg).

ClimaCore.MatrixFields.preconditioner_cacheMethod
preconditioner_cache(P_alg, A, b)

Allocates the cache required to solve the equation P * x = b, where P is the preconditioner generated by the PreconditionerAlgorithm P_alg for A.

ClimaCore.MatrixFields.replace_lazy_operatorMethod
replace_lazy_operator(space, lazy_op)

Generates an instance of Base.AbstractBroadcasted that corresponds to the expression lazy_op.(), where the broadcast in which this expression appears is being evaluated on the given space. Note that the staggering (CellCenter or CellFace) of this space depends on the specifics of the broadcast and is not predetermined.

ClimaCore.MatrixFields.rmul_return_typeMethod
rmul_return_type(X, Y)

Computes the return type of rmul_with_projection(x, y, lg), where x isa X and y isa Y. This can also be used to obtain the return type of rmul(x, y), although rmul(x, y) will throw an error when projection is necessary.

Note that this is equivalent to calling the internal function _return_type: Base._return_type(rmul_with_projection, Tuple{X, Y, LG}), where lg isa LG.

ClimaCore.MatrixFields.solver_algorithmMethod
solver_algorithm(P_alg)

A FieldMatrixSolverAlgorithm that can be used to solve P * x = b for x, where P is the preconditioner generated by the PreconditionerAlgorithm P_alg. If P_alg is nothing instead of a PreconditionerAlgorithm, or if P is a diagonal matrix (and no solver is required to invert it), this returns nothing.

ClimaCore.MatrixFields.@nameMacro
@name(expr)

Shorthand for constructing a FieldName. Some examples include

  • name = @name(), in which case get_field(x, name) returns x
  • name = @name(a), in which case get_field(x, name) returns x.a
  • name = @name(a.b.c), in which case get_field(x, name) returns x.a.b.c
  • name = @name(a.b.c.:(1).d), in which case get_field(x, name) returns x.a.b.c.:(1).d

This macro is preferred over the FieldName constructor because it checks whether expr is a syntactically valid chain of getproperty calls before calling the constructor.

ClimaCore.Utilities.UnrolledFunctionsModule
UnrolledFunctions

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

The functions exported by this module are

  • unrolled_map(f, values, [values2]): alternative to map
  • unrolled_any(f, values): alternative to any
  • unrolled_all(f, values): alternative to all
  • unrolled_filter(f, values): alternative to filter
  • unrolled_foreach(f, values): alternative to foreach
  • unrolled_in(value, values): alternative to in
  • unrolled_unique(values): alternative to unique
  • unrolled_flatten(values): alternative to Iterators.flatten
  • unrolled_flatmap(f, values): alternative to Iterators.flatmap
  • unrolled_product(values1, values2): alternative to Iterators.product
  • unrolled_findonly(f, values): checks that only one value satisfies f, and then returns that value
  • unrolled_split(f, values): returns a tuple that contains the result of calling unrolled_filter with f and the result of calling it with !f
  • unrolled_take(values, ::Val{N}): alternative to Iterators.take, but with an Int wrapped in a Val as the second argument instead of a regular Int; this usually compiles more quickly than values[1:N]
  • unrolled_drop(values, ::Val{N}): alternative to Iterators.drop, but with an Int wrapped in a Val as the second argument instead of a regular Int; this usually compiles more quickly than values[(end - N + 1):end]
ClimaCore.SpacesModule
Meshes
  • domain
  • topology
  • coordinates
  • metric terms (inverse partial derivatives)
  • quadrature rules and weights

References / notes

ClimaCore.Spaces.areaMethod
Spaces.area(space::Spaces.AbstractSpace)

The length/area/volume of space. This is computed as the sum of the quadrature weights $W_i$ multiplied by the Jacobian determinants $J_i$:

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

If space is distributed, this uses a ClimaComms.allreduce operation.

ClimaCore.Spaces.unique_nodesMethod
unique_nodes(space::SpectralElementSpace2D)

An iterator over the unique nodes of space. Each node is represented by the first ((i,j), e) triple.

This function is experimental, and may change in future.

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

1). Finish communications.

2). Call Spaces.load_from_recv_buffer! After the communication is complete, this adds data from the recv buffer to the corresponding location in perimeter_data. For ghost vertices, this data is added only to the representative vertices. The values are then scattered to other local vertices corresponding to each unique ghost vertex in dss_local_ghost.

3). Call Spaces.dss_untransform! on all local elements. This transforms the DSS'd local vectors back to Covariant12 vectors, and copies the DSS'd data from the perimeter_data to data.

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

1). Apply Spaces.dss_transform! on interior elements. Local elements are split into interior and perimeter elements to facilitate overlapping of communication with computation.

2). Probe communication

3). Spaces.dss_local! computes the weighted DSS on local vertices and faces.

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

It comprises of the following steps:

1). Apply Spaces.dss_transform! on perimeter elements. This weights and tranforms vector fields to physical basis if needed. Scalar fields are weighted. The transformed and/or weighted perimeter data is stored in perimeter_data.

2). Apply Spaces.dss_local_ghost! This computes partial weighted DSS on ghost vertices, using only the information from local vertices.

3). Spaces.fill_send_buffer! Loads the send buffer from perimeter_data. For unique ghost vertices, only data from the representative ghost vertices which store result of "ghost local" DSS are loaded.

4). Start DSS communication with neighboring processes

ClimaCore.Topologies.create_dss_bufferMethod
create_dss_buffer(
    data::Union{DataLayouts.IJFH{S, Nij}, DataLayouts.VIJFH{S, Nij}},
    hspace::AbstractSpectralElementSpace,
) where {S, Nij}

Creates a DSSBuffer for the field data corresponding to data

ClimaCore.Quadratures.differentiation_matrixMethod
differentiation_matrix(r::SVector{Nq, T}) where {Nq, T}

The spectral differentiation matrix for the Lagrange polynomial of degree Nq-1 interpolating at points r.

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

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

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

ClimaCore.Quadratures.interpolation_matrixMethod
interpolation_matrix(x::SVector, r::SVector{Nq})

The matrix which interpolates the Lagrange polynomial of degree Nq-1 through the points r, to points x. The matrix coefficients are computed using the Barycentric formula of Berrut2004, section 4:

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

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

ClimaCore.Quadratures.orthonormal_polyMethod
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.

ClimaCore.Domains.IntervalDomainMethod
IntervalDomain(coord⁻, coord⁺; periodic=true)
IntervalDomain(coord⁻, coord⁺; boundary_names::Tuple{Symbol,Symbol})

Construct a IntervalDomain, the closed interval is given by coord⁻, coord⁺ coordinate arguments.

Either a periodic or boundary_names keyword argument is required.

ClimaCore.Domains.RectangleDomainMethod
RectangleDomain(x1::ClosedInterval, x2::ClosedInterval;
    x1boundary::Tuple{Symbol,Symbol},
    x2boundary::Tuple{Symbol,Symbol},
    x1periodic = false,
    x2periodic = false,
)

Construct a RectangularDomain in the horizontal. If a given x1 or x2 boundary is not periodic, then x1boundary or x2boundary boundary name keyword arguments must be supplied.

ClimaCore.Domains.boundary_namesFunction
boundary_names(obj::Union{AbstractDomain, AbstractMesh, AbstractTopology})

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

ClimaCore.Fields.FieldVectorType
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, ....

ClimaCore.Fields.ScalarWrapperType
Fields.ScalarWrapper(val) <: AbstractArray{T,0}

This is a wrapper around scalar values that allows them to be mutated as part of a FieldVector. A call getproperty on a FieldVector with this component will return a scalar, instead of the boxed object.

Base.fill!Method
fill!(field::Field, value)

Fill field with value.

Base.fillMethod
fill(value, space::AbstractSpace)

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

Base.maximumMethod
maximum([f=identity,]v::Field)

Approximate maximum of v or f.(v) over the domain.

If v is a distributed field, this uses a ClimaComms.allreduce operation.

Base.onesMethod
ones(space::AbstractSpace)

Construct a field on space that is one everywhere.

Base.sumMethod
sum([f=identity,] v::Field)

Approximate integration of v or f.(v) over the domain. In an AbstractSpectralElementSpace, an integral over the entire space is computed by summation over the elements of the integrand multiplied by the Jacobian determinants and the quadrature weights at each node within an element. Hence, sum is computed by summation of the field values multiplied by the Jacobian determinants and quadrature weights:

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

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

If v is a distributed field, this uses a ClimaComms.allreduce operation.

Base.zerosMethod
zeros(space::AbstractSpace)

Construct a field on space that is zero everywhere.

ClimaCore.Fields.backing_arrayMethod
backing_array(x)

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

ClimaCore.Fields.bycolumnMethod
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.

Note

On GPUs this will simply evaluate f once with colidx=: (i.e. it doesn't perform evaluation by columns). This may change in future.

Example

∇ = GradientF2C()
div = DivergenceC2F()

bycolumn(axes(f)) do colidx
    @. ∇f[colidx] = ∇(f[colidx])
    @. df[colidx] = div(∇f[colidx])
end
ClimaCore.Fields.field_iteratorMethod
field_iterator(::Union{Field, FieldVector})

Returns an iterable field, that recursively calls getproperty for all of the propertynames, and returns a Tuple of

  • the individual scalar field and
  • the tuple chain used to reach the scalar field
ClimaCore.Fields.local_sumMethod
Fields.local_sum(v::Field)

Compute the approximate integral of v over the domain local to the current context.

See sum for the integral over the full domain.

ClimaCore.Fields.property_chainsMethod
property_chains

An array of "property chains", used to recursively getproperty until a single scalar field is reached.

A property chain may be, for example (:model, :submodel, :temperature) where model.submodel.temperature is a scalar field.

ClimaCore.Fields.set!Function
set!(f::Function, field::Field, args = ())

Apply function f to populate values in field field. f must have a function signature with signature f(::LocalGeometry[, args...]). Additional arguments may be passed to f with args.

ClimaCore.Fields.unwrapMethod
Fields.unwrap(x::T)

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

ClimaCore.Fields.wrapMethod
Fields.wrap(x)

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

ClimaCore.Fields.Δz_fieldMethod
Δz_field(field::Field)
Δz_field(space::AbstractSpace)

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

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

Apply weighted direct stiffness summation (DSS) to f. This operates in-place (i.e. it modifies the f). ghost_buffer contains the necessary information for communication in a distributed setting, see Spaces.create_dss_buffer.

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

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

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

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

Then the problem can be written as

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

which reduces to

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

LinearAlgebra.normFunction
norm(v::Field, p=2; normalize=true)

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

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

where $|v|$ is defined to be the absolute value if $v$ is a scalar-valued Field, or the 2-norm if it is a vector-valued Field or composite Field (see LinearAlgebra.norm). Similar to sum and mean, in an AbstractSpectralElementSpace, this is computed by summation of the field values multiplied by the Jacobian determinants and quadrature weights. If normalize=true (the default), then internally the discrete norm is divided by the sum of the Jacobian determinants and quadrature weights:

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

If p=Inf, then the norm is the maximum of the absolute values

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

Consequently all norms should have the same units for all $p$ (being the same as calling norm on a single value).

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

Statistics.meanMethod
mean([f=identity, ]v::Field)

The mean value of field or f.(field) over the domain, weighted by area. Similar to sum, in an AbstractSpectralElementSpace, this is computed by summation of the field values multiplied by the Jacobian determinants and quadrature weights:

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

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

If v is a distributed field, this uses a ClimaComms.allreduce operation.

ClimaCore.Meshes.AbstractCubedSphereType
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

ClimaCore.Meshes.AbstractMeshType
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:

The following types/methods are provided by AbstractMesh:

ClimaCore.Meshes.ConformalCubedSphereType
ConformalCubedSphere <: AbstractCubedSphere

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

Constructors

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

Constuct a ConformalCubedSphere on domain with ne elements across each panel.

ClimaCore.Meshes.EquiangularCubedSphereType
EquiangularCubedSphere <: AbstractCubedSphere

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

Constructors

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

Constuct an EquiangularCubedSphere on domain with ne elements across each panel.

ClimaCore.Meshes.ExponentialStretchingType
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)
ClimaCore.Meshes.GeneralizedExponentialStretchingType
GeneralizedExponentialStretching(dz_bottom::FT, dz_top::FT)

Apply a generalized form of exponential stretching to the domain when constructing elements. dz_bottom and dz_top are target element grid spacings at the bottom and at the top of the vertical column domain (m). In typical atmosphere configurations, dz_bottom is the smallest grid spacing and dz_top the largest one. On the other hand, for typical land configurations, dz_bottom is the largest grid spacing and dz_top the smallest one.

For land configurations, use reverse_mode = true (default value false).

Then, the user can define a generalized stretched mesh via

ClimaCore.Meshes.IntervalMesh(interval_domain, GeneralizedExponentialStretching(dz_bottom, dz_top); nelems::Int, reverse_mode = false)
ClimaCore.Meshes.RectilinearMeshType
RectilinearMesh <: AbstractMesh2D

Constructors

RectilinearMesh(domain::RectangleDomain, n1, n2)

Construct a RectilinearMesh of equally-spaced n1 by n2 elements on domain.

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

Construct the product mesh of intervalmesh1 and intervalmesh2.

ClimaCore.Meshes.boundary_face_nameFunction
Meshes.boundary_face_name(mesh::AbstractMesh, elem, face::Int)::Union{Symbol,Nothing}

The name of the boundary facing face of element elem, or nothing if it is not on the boundary.

ClimaCore.Meshes.containing_elementFunction
elem = Meshes.containing_element(mesh::AbstractMesh, coord)

The element elem in mesh containing the coordinate coord. If the coordinate falls on the boundary between two or more elements, an arbitrary element is chosen.

ClimaCore.Meshes.containing_panelMethod
panel = cubedspherepanel(coord::Geometry.Cartesian123Point)

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

ClimaCore.Meshes.coordinatesMethod
Meshes.coordinates(mesh, elem, vert::Int)
Meshes.coordinates(mesh, elem, ξ::SVector)

Return the physical coordinates of a point in an element elem of mesh. The position of the point can either be a vertex number vert or the coordinates ξ in the reference element.

ClimaCore.Meshes.elementsFunction
Meshes.elements(mesh::AbstractMesh)

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

ClimaCore.Meshes.face_connectivity_matrixFunction
M = Meshes.face_connectivity_matrix(mesh, elemorder = elements(mesh))

Construct a Bool-valued SparseCSCMatrix containing the face connections of mesh. Elements are indexed according to elemorder.

Note that M[i,i] == true only if two distinct faces of element i are connected.

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

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

ClimaCore.Meshes.linearindicesMethod
Meshes.linearindices(elemorder)

Given a data structure elemorder[i] = elem that orders elements, construct the inverse map from orderindex = linearindices(elemorder) such that orderindex[elem] = i.

This will try to use the most efficient structure available.

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

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

ClimaCore.Meshes.reference_coordinatesFunction
ξ = Meshes.reference_coordinates(mesh::AbstractMesh, elem, coord)

An SVector of coordinates in the reference element such that

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

This can be used for interpolation to a specific point.

ClimaCore.Meshes.refindexMethod
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]`.

ClimaCore.Meshes.to_panelMethod
to_panel(panel, coord1::Geometry.Cartesian123Point)

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

ClimaCore.Meshes.truncate_meshMethod
truncate_mesh(
    parent_mesh::AbstractMesh,
    trunc_domain::IntervalDomain{CT},
)

Constructs an IntervalMesh, truncating the given parent_mesh defined on a truncated trunc_domain. The truncation preserves the number of degrees of freedom covering the space from the trunc_domain's z_bottom to z_top, adjusting the stretching.

ClimaCore.Meshes.vertex_connectivity_matrixFunction
M = Meshes.vertex_connectivity_matrix(mesh, elemorder = elements(mesh))

Construct a Bool-valued SparseCSCMatrix containing the vertex connections of mesh. Elements are indexed according to elemorder.

Note that M[i,i] == true only if two distinct vertices of element i are connected.

ClimaCore.Topologies.DSSBufferType
DSSBuffer{G, D, A, B}

Fields

  • graph_context: ClimaComms graph context for communication

  • perimeter_data: Perimeter DataLayout object: typically a VIFH{TT,Np}, where TT is the transformed type, and Np is the length of the perimeter

  • send_data: send buffer AbstractVector{FT}

  • recv_data: recv buffer AbstractVector{FT}

  • send_buf_idx: indexing array for loading send buffer from perimeter_data

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

  • scalarfidx: field id for all scalar fields stored in the data array

  • covariant12fidx: field id for all covariant12vector fields stored in the data array

  • contravariant12fidx: field id for all contravariant12vector fields stored in the data array

  • internal_elems: internal local elements (lidx)

  • perimeter_elems: local elements (lidx) located on process boundary

ClimaCore.Topologies.Perimeter2DMethod
Perimeter2D(Nq)

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

ClimaCore.Topologies.Topology2DType
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
ClimaCore.Topologies.VertexIteratorType
Topologies.VertexIterator

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

ClimaCore.Topologies.boundary_facesFunction
boundary_faces(topology, boundarytag)

An iterator over the faces of topology which face the boundary with tag boundarytag. Each element of the iterator is an (elem, face) pair.

ClimaCore.Topologies.boundary_tagFunction
boundary_tag(topology, name::Symbol)

The boundary tag of the topology for boundary name name. A boundary tag is an integer that uniquely identifies a boundary.

ClimaCore.Topologies.boundary_tagsFunction
boundary_tags(topology)

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

ClimaCore.Topologies.create_dss_bufferMethod
create_dss_buffer(
    data::Union{DataLayouts.IJFH{S, Nij}, DataLayouts.VIJFH{S, Nij}},
    topology::Topology2D,
    local_geometry = nothing,
    local_weights = nothing,
) where {S, Nij}

Creates a DSSBuffer for the field data corresponding to data

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

Sets the value for all local vertices of each unique ghost vertex, in perimeter_data, to that of the representative ghost vertex.

Part of ClimaCore.Spaces.weighted_dss!.

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

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

Part of ClimaCore.Spaces.weighted_dss!.

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

Transforms vectors from Covariant axes to physical (local axis), weights the data at perimeter nodes, and stores result in the perimeter_data array. This function calls the appropriate version of dss_transform! based on the data layout of the input arguments.

Arguments:

  • dss_buffer: DSSBuffer generated by create_dss_buffer function for field data
  • data: field data
  • local_geometry: local metric information defined at each node
  • weight: local dss weights for horizontal space
  • perimeter: perimeter iterator
  • localelems: list of local elements to perform transformation operations on

Part of ClimaCore.Spaces.weighted_dss!.

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

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

Arguments:

  • perimeter_data: contains the perimeter field data, represented on the physical axis, corresponding to the full field data in data
  • data: field data
  • ∂ξ∂x: partial derivatives of the map from x to ξ: ∂ξ∂x[i,j] is ∂ξⁱ/∂xʲ
  • weight: local dss weights for horizontal space
  • perimeter: perimeter iterator
  • scalarfidx: field index for scalar fields in the data layout
  • covariant12fidx: field index for Covariant12 vector fields in the data layout
  • localelems: list of local elements to perform transformation operations on

Part of ClimaCore.Spaces.weighted_dss!.

ClimaCore.Topologies.dss_transformMethod
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!.

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

Transforms the DSS'd local vectors back to Covariant12 vectors, and copies the DSS'd data from the perimeter_data to data. This function calls the appropriate version of dss_transform! function based on the data layout of the input arguments.

Arguments:

  • dss_buffer: DSSBuffer generated by create_dss_buffer function for field data
  • data: field data
  • local_geometry: local metric information defined at each node
  • perimeter: perimeter iterator
  • localelems: list of local elements to perform transformation operations on

Part of ClimaCore.Spaces.weighted_dss!.

ClimaCore.Topologies.face_node_indexFunction
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.

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

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

Part of ClimaCore.Spaces.weighted_dss!.

ClimaCore.Topologies.ghost_facesMethod
ghost_faces(topology::AbstractTopology)

An iterator over the ghost faces of topology. Each element of the iterator is a 5-tuple the form

(elem1, face1, elem2, face2, reversed)

where elemX, faceX are the element and face numbers, and reversed indicates whether they have opposing orientations.

ClimaCore.Topologies.ghost_neighboring_elementsFunction
ghost_neighboring_elements(topology::AbstractTopology, ridx::Integer)

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

ClimaCore.Topologies.ghost_verticesFunction
ghost_vertices(topology)

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

ClimaCore.Topologies.interior_facesMethod
interior_faces(topology::AbstractTopology)

An iterator over the interior faces of topology. Each element of the iterator is a 5-tuple the form

(elem1, face1, elem2, face2, reversed)

where elemX, faceX are the element and face numbers, and reversed indicates whether they have opposing orientations.

ClimaCore.Topologies.load_from_recv_buffer!Method
load_from_recv_buffer!(::ClimaComms.AbstractCPUDevice, dss_buffer::DSSBuffer)

Adds data from the recv buffer to the corresponding location in perimeter_data. For ghost vertices, this data is added only to the representative vertices. The values are then scattered to other local vertices corresponding to each unique ghost vertex in dss_local_ghost.

Part of ClimaCore.Spaces.weighted_dss!.

ClimaCore.Topologies.local_neighboring_elementsFunction
local_neighboring_elements(topology::AbstractTopology, lidx::Integer)

An iterator of the local element indices (lidx) of the local elements which are neighbors of the local element lidx in topology (excluding lidx itself).

ClimaCore.Topologies.local_verticesFunction
local_vertices(topology)

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

ClimaCore.Topologies.opposing_faceFunction
(opelem, opface, reversed) = opposing_face(topology, elem, face)

The opposing face of face number face of element elem in topology.

  • opelem is the opposing element number, 0 for a boundary, negative for a ghost element
  • opface is the opposite face number, or boundary face number if a boundary
  • reversed indicates whether the opposing face has the opposite orientation.
ClimaCore.Topologies.spacefillingcurveMethod
spacefillingcurve(mesh::Meshes.AbstractCubedSphere)

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

ClimaCore.Utilities.PlusHalfType
PlusHalf(i)

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

Supports +, - and inequalities.

See also half.

ClimaCore.InputOutput.HDF5ReaderType
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

Usage

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

To explore the contents of the reader, use either

julia> reader |> propertynames

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

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

Once "unpacked" as shown above, ClimaCorePlots or ClimaCoreMakie can be used to visualise fields. ClimaCoreTempestRemap supports interpolation onto user-specified grids if necessary.

ClimaCore.InputOutput.HDF5WriterType
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.

Note

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

Interface

write!

Usage

writer = InputOutput.HDF5Writer(filename)
InputOutput.write!(writer, Y, "Y")
close(writer)
ClimaCore.InputOutput.read_attributesMethod
read_attributes(reader::AbstractReader, name::AbstractString, data::Dict)

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

ClimaCore.InputOutput.read_domainMethod
read_domain(reader::AbstractReader, name)

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

ClimaCore.InputOutput.read_fieldMethod
read_field(reader, name)

Reads a Field or FieldVector named name from reader. Fields are not cached, so that reading the same field multiple times will create multiple distinct objects.

ClimaCore.InputOutput.read_gridMethod
read_grid(reader::AbstractReader, name)

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

ClimaCore.InputOutput.read_meshMethod
read_mesh(reader::AbstractReader, name)

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

ClimaCore.InputOutput.read_spaceMethod
read_space(reader::AbstractReader, name)

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

ClimaCore.InputOutput.read_topologyMethod
read_topology(reader::AbstractReader, name)

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

ClimaCore.InputOutput.write!Function
write!(writer::AbstractWriter, obj[, preferredname])

Write the object obj using writer. An optional preferredname can be provided, otherwise defaultname will be used to generate a name. The name of the object will be returned.

A cache of domains, meshes, topologies and spaces is kept: if one of these objects has already been written, then the file will not be modified: instead the name under which the object was first written will be returned. Note that Fields and FieldVectors are not cached, and so can be written multiple times.

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

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

ClimaCore.InputOutput.write_attributes!Method
write_attributes!(writer::AbstractWriter, name::AbstractString, data::Dict)

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

ClimaCore.DeviceSideContextType
DeviceSideContext()

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

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

ClimaCore.DeviceSideDeviceType
DeviceSideDevice()

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

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

ClimaCore.columnFunction
column(data::AbstractData, i::Integer)

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

ClimaCore.enable_threadingFunction
enable_threading()

By default returns false signifying threading is disabled. Enable the threading runtime by redefining the method at the toplevel your experiment file:

import ClimaCore: enable_threading
enable_threading() = true

and running julia with julia --nthreads=N ...

This function is deprecated in version 0.10.42. Please use the ClimaComms context for threading.

ClimaCore.slabFunction
slab(data::AbstractData, h::Integer)

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

ClimaCore.Operators.AbstractOperatorType
AbstractOperator

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

ClimaCore.Operators.AdvectionC2CType
A = AdvectionC2C(;boundaries)
A.(v, θ)

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

It uses the following stencil

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

Supported boundary conditions:

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

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

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

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

ClimaCore.Operators.AdvectionF2FType
A = AdvectionF2F(;boundaries)
A.(v, θ)

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

It uses the following stencil

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

No boundary conditions are currently supported.

ClimaCore.Operators.CUDASpectralStyleType
CUDASpectralStyle()

Applies spectral-element operations by using threads for each node, and synchronizing when they occur. This is used for GPU kernels.

ClimaCore.Operators.CurlType
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

ClimaCore.Operators.CurlC2FType
C = CurlC2F(;boundaryname=boundarycondition...)
C.(v)

Compute the vertical-derivative contribution to the curl of a center-valued covariant vector field v. It acts on the horizontal covariant components of v (that is it only depends on $v₁$ and $v₂$), and will return a face-valued horizontal contravariant vector field (that is $C(v)³ = 0$).

Specifically it approximates:

\[\begin{align*} C(v)^1 &= -\frac{1}{J} \frac{\partial v_2}{\partial \xi^3} \\ C(v)^2 &= \frac{1}{J} \frac{\partial v_1}{\partial \xi^3} \\ \end{align*}\]

using the stencils

\[\begin{align*} C(v)[i]^1 &= - \frac{1}{J[i]} (v₂[i+\tfrac{1}{2}] - v₂[i-\tfrac{1}{2}]) \\ C(v)[i]^2 &= \frac{1}{J[i]} (v₁[i+\tfrac{1}{2}] - v₁[i-\tfrac{1}{2}]) \end{align*}\]

where $v₁$ and $v₂$ are the 1st and 2nd covariant components of $v$, and $J$ is the Jacobian determinant.

The following boundary conditions are supported:

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

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

  • SetCurl(v⁰): enforce the curl operator output at the boundary to be the contravariant vector v⁰.
ClimaCore.Operators.DivergenceType
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

ClimaCore.Operators.DivergenceC2FType
D = DivergenceC2F(;boundaryname=boundarycondition...)
D.(v)

Compute the vertical contribution to the divergence of a center-valued field vector v, returning a face-valued scalar field, using the stencil

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

where Jv³ is the Jacobian multiplied by the third contravariant component of v.

The following boundary conditions are supported:

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

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

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

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

ClimaCore.Operators.DivergenceF2CType
D = DivergenceF2C(;boundaryname=boundarycondition...)
D.(v)

Compute the vertical contribution to the divergence of a face-valued field vector v, returning a center-valued scalar field, using the stencil

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

where Jv³ is the Jacobian multiplied by the third contravariant component of v.

The following boundary conditions are supported:

  • by default, the value of v at the boundary face will be used.
  • SetValue(v₀): calculate the divergence assuming the value at the boundary is v₀. For the left boundary, this becomes:

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

  • Extrapolate(): set the value at the center closest to the boundary to be the same as the neighbouring interior value. For the left boundary, this becomes:

\[D(v)[1]³ = D(v)[2]³\]

ClimaCore.Operators.FCTBorisBookType
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.
Note

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

ClimaCore.Operators.FCTZalesakType
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.
Note

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

ClimaCore.Operators.GradientType
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

ClimaCore.Operators.GradientC2FType
G = GradientC2F(;boundaryname=boundarycondition...)
G.(x)

Compute the gradient of a center-valued field x, returning a face-valued Covariant3 vector field, using the stencil:

\[G(x)[i]^3 = x[i+\tfrac{1}{2}] - x[i-\tfrac{1}{2}]\]

The following boundary conditions are supported:

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

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

  • SetGradient(v₀): set the value of the gradient at the boundary to be v₀. For the left boundary, this becomes:

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

ClimaCore.Operators.GradientF2CType
G = GradientF2C(;boundaryname=boundarycondition...)
G.(x)

Compute the gradient of a face-valued field x, returning a center-valued Covariant3 vector field, using the stencil:

\[G(x)[i]^3 = x[i+\tfrac{1}{2}] - x[i-\tfrac{1}{2}]\]

We note that the usual division factor $1 / \Delta z$ that appears in a first-order finite difference operator is accounted for in the LocalVector basis. Hence, users need to cast the output of the GradientF2C to a UVector, VVector or WVector, according to the type of domain on which the operator is defined.

The following boundary conditions are supported:

  • by default, the value of x at the boundary face will be used.
  • SetValue(x₀): calculate the gradient assuming the value at the boundary is x₀. For the left boundary, this becomes:

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

  • Extrapolate(): set the value at the center closest to the boundary

to be the same as the neighbouring interior value. For the left boundary, this becomes:

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

ClimaCore.Operators.InterpolateType
i = Interpolate(space)
i.(f)

Interpolates f to the space. If space has equal or higher polynomial degree as the space of f, this is exact, otherwise it will be lossy.

In matrix form, it is the linear operator

\[I = \bigotimes_i I_i\]

where $I_i$ is the barycentric interpolation matrix in the $i$th dimension.

See also Restrict.

ClimaCore.Operators.InterpolateC2FType
I = InterpolateC2F(;boundaries..)
I.(x)

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

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

Supported boundary conditions are:

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

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

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

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

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

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

ClimaCore.Operators.LeftBiased3rdOrderC2FType
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:

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

ClimaCore.Operators.LeftBiased3rdOrderF2CType
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:

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

ClimaCore.Operators.LeftBiasedC2FType
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:

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

ClimaCore.Operators.LeftBiasedF2CType
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:

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

ClimaCore.Operators.RestrictType
r = Restrict(space)
r.(f)

Computes the projection of a field f on $\mathcal{V}_0$ to a lower degree polynomial space space ($\mathcal{V}_0^*$). space must be on the same topology as the space of f, but have a lower polynomial degree.

It is defined as the field $\theta \in \mathcal{V}_0^*$ such that for all $\phi \in \mathcal{V}_0^*$

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

In matrix form, this is

\[\phi^\top W^* J^* \theta = (I \phi)^\top WJ f\]

where $W^*$ and $J^*$ are the quadrature weights and Jacobian determinant of $\mathcal{V}_0^*$, and $I$ is the interpolation operator (see Interpolate) from $\mathcal{V}_0^*$ to $\mathcal{V}_0$. This reduces to

\[\theta = (W^* J^*)^{-1} I^\top WJ f\]

ClimaCore.Operators.RightBiased3rdOrderC2FType
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:

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

ClimaCore.Operators.RightBiased3rdOrderF2CType
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:

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

ClimaCore.Operators.RightBiasedC2FType
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:

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

ClimaCore.Operators.RightBiasedF2CType
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:

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

ClimaCore.Operators.SetGradientType
SetGradient(val)

Set the gradient at the boundary to be val. In the case of gradient operators this will set the output value of the gradient.

ClimaCore.Operators.SetValueType
SetValue(val)

Set the value at the boundary to be val. In the case of gradient operators, this will set the input value from which the gradient is computed.

ClimaCore.Operators.SpectralBroadcastedType
SpectralBroadcasted{Style}(op, args[,axes[, work]])

This is similar to a Base.Broadcast.Broadcasted object, except it contains space for an intermediate work storage.

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

ClimaCore.Operators.StencilBroadcastedType
StencilBroadcasted{Style}(op, args[,axes[, work]])

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

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

ClimaCore.Operators.Upwind3rdOrderBiasedProductC2FType
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.

Note

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

ClimaCore.Operators.UpwindBiasedProductC2FType
U = UpwindBiasedProductC2F(;boundaries)
U.(v, x)

Compute the product of the face-valued vector field v and a center-valued field x at cell faces by upwinding x according to the direction of v.

More precisely, it is computed based on the sign of the 3rd contravariant component, and it returns a Contravariant3Vector:

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

where $\boldsymbol{e}_3$ is the 3rd covariant basis vector.

Supported boundary conditions are:

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

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

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

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

ClimaCore.Operators.WeakCurlType
wcurl = WeakCurl()
wcurl.(u)

Computes the "weak curl" on each element of a covariant vector field u.

Note: The vector field $u$ needs to be excliclty converted to a CovaraintVector, as then the WeakCurl is independent of the local metric tensor.

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

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

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

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

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

In matrix form, this becomes

\[{\phi_i}^\top W J \theta^i = (J^{-1} \epsilon^{kji} D_j \phi_i)^\top W J u_k\]

which, by using the anti-symmetry of the Levi-Civita symbol, reduces to

\[\theta^i = - \epsilon^{ijk} (WJ)^{-1} D_j^\top W u_k\]

ClimaCore.Operators.WeakDivergenceType
wdiv = WeakDivergence()
wdiv.(u)

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

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

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

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

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

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

It can be written in matrix form as

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

which reduces to

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

where

  • $J$ is the diagonal Jacobian matrix
  • $W$ is the diagonal matrix of quadrature weights
  • $D_i$ is the derivative matrix along the $i$th dimension
ClimaCore.Operators.WeakGradientType
wgrad = WeakGradient()
wgrad.(f)

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

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

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

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

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

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

In matrix form, this becomes

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

which reduces to

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

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

ClimaCore.Operators.WeightedInterpolateC2FType
WI = WeightedInterpolateC2F(; boundaries)
WI.(w, x)

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

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

Supported boundary conditions are:

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

These have the same stencil as in InterpolateC2F.

ClimaCore.Operators.WeightedInterpolateF2CType
WI = WeightedInterpolateF2C(; boundaries)
WI.(w, x)

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

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

No boundary conditions are required (or supported)

ClimaCore.Operators.add_numerical_flux_internal!Method
add_numerical_flux_internal!(fn, dydt, args...)

Add the numerical flux at the internal faces of the spectral space mesh.

The numerical flux is determined by evaluating

fn(normal, argvals⁻, argvals⁺)

where:

  • normal is the unit normal vector, pointing from the "minus" side to the "plus" side
  • argvals⁻ is the tuple of values of args on the "minus" side of the face
  • argvals⁺ is the tuple of values of args on the "plus" side of the face

and should return the net flux from the "minus" side to the "plus" side.

For consistency, it should satisfy the property that

fn(normal, argvals⁻, argvals⁺) == -fn(-normal, argvals⁺, argvals⁻)

See also:

ClimaCore.Operators.column_integral_definite!Method
column_integral_definite!(∫field::Field, ᶜfield::Field)

Sets ∫field${}= \int_0^{z_{max}}\,$ᶜfield$(z)\,dz$, where $z_{max}$ is the value of z at the top of the domain. The input ᶜfield must lie on a cell-center space, and the output ∫field must lie on the corresponding horizontal space.

ClimaCore.Operators.column_integral_indefinite!Method
column_integral_indefinite!(ᶠ∫field::Field, ᶜfield::Field)

Sets ᶠ∫field$(z) = \int_0^z\,$ᶜfield$(z')\,dz'$. The input ᶜfield must lie on a cell-center space, and the output ᶠ∫field must lie on the corresponding cell-face space.

column_integral_indefinite!(
    f::Function,
    ᶠ∫field::Fields.ColumnField,
    ϕ₀ = 0,
    average = (ϕ⁻, ϕ⁺) -> (ϕ⁻ + ϕ⁺) / 2,
)

The indefinite integral ᶠ∫field = ϕ(z) = ∫ f(ϕ,z) dz given:

  • f the integral integrand function (which may be a function)
  • ᶠ∫field the resulting (scalar) field ϕ(z)
  • ϕ₀ (optional) the boundary condition
  • average (optional) a function to compute the cell center average between two cell faces (ϕ⁻, ϕ⁺).
ClimaCore.Operators.column_mapreduce!Method
column_mapreduce!(fn, op, reduced_field::Field, fields::Field...)

Applies mapreduce along the vertical direction. The input fields must all lie on the same space, and the output reduced_field must lie on the corresponding horizontal space. The function fn is mapped over every input, and the function op is used to reduce the outputs of fn.

ClimaCore.Operators.column_thomas_solve!Method
column_thomas_solve!(A, b)

Solves the linear system A * x = b, where A is a tri-diagonal matrix (represented by a Field of tri-diagonal matrix rows), and where b is a vector (represented by a Field of numbers). The data in b is overwritten with the solution x, and the upper diagonal of A is also overwritten with intermediate values used to compute x. The algorithm is described here: https://en.wikipedia.org/wiki/Tridiagonalmatrixalgorithm.

ClimaCore.Operators.left_interior_idxMethod
left_interior_idx(space::AbstractSpace, op::FiniteDifferenceOperator, bc::AbstractBoundaryCondition, args..)

The index of the left-most interior point of the operator op with boundary bc when used with arguments args.... By default, this is

left_idx(space) + boundary_width(op, bc)

but can be overwritten for specific stencil types (e.g. if the stencil is assymetric).

ClimaCore.Operators.linear_remap_opMethod
linear_remap_op(target::AbstractSpace, source::AbstractSpace)

Computes linear remapping operator R for remapping from source to target spaces.

Entry R_{ij} gives the contribution weight to the target node i from source node j; nodes are indexed by their global position, determined by both element index and nodal order within the element.

ClimaCore.Operators.local_weightsMethod
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.

ClimaCore.Operators.matrix_interpolateFunction
matrix_interpolate(field, quadrature)

Computes the tensor product given a uniform quadrature out = (M ⊗ M) * in on each element. Returns a 2D Matrix for plotting / visualizing 2D Fields.

ClimaCore.Operators.matrix_interpolateMethod
matrix_interpolate(field, Nu::Integer)

Computes the tensor product given a uniform quadrature degree of Nu on each element. Returns a 2D Matrix for plotting / visualizing 2D Fields.

ClimaCore.Operators.overlap_weights!Method
overlap_weights!(J_ov, target, source, i, j, coords_t, coords_s)

Computes the overlap weights for a pair of source and target elements.

The spatial overlap of element i on the target space and element j of the source space, computed for each nodal pair, approximating Ullrich2015 eq. 19 via quadrature.

ClimaCore.Operators.remap!Method
remap!(target_field::Field, R::LinearRemap, source_field::Field)

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

ClimaCore.Operators.remapMethod
remap(R::LinearRemap, source_field::Field)

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

ClimaCore.Operators.resolve_operatorMethod
resolve_operator(bc, slabidx)

Recursively evaluate any operators in bc at slabidx, replacing any SpectralBroadcasted objects.

  • if bc is a regular Broadcasted object, return a new Broadcasted with resolve_operator called on each arg
  • if bc is a regular SpectralBroadcasted object:
  • call resolve_operator called on each arg
  • call apply_operator, returning the resulting "pseudo Field": a Field with an

IF/IJF data object.

  • if bc is a Field, return that
ClimaCore.Operators.resolve_shmem!Method
resolve_shmem!(obj, ij, slabidx)

Recursively stores the arguments to all operators into shared memory, at the given indices (if they are valid).

As this calls sync_threads(), it should be called collectively on all threads at the same time.

ClimaCore.Operators.right_interior_idxMethod
right_interior_idx(space::AbstractSpace, op::FiniteDifferenceOperator, bc::AbstractBoundaryCondition, args..)

The index of the right-most interior point of the operator op with boundary bc when used with arguments args.... By default, this is

right_idx(space) + boundary_width(op, bc)

but can be overwritten for specific stencil types (e.g. if the stencil is assymetric).

ClimaCore.Operators.stencil_interiorFunction
stencil_interior(::Op, loc, space, idx, args...)

Defines the stencil of the operator Op in the interior of the domain at idx; args are the input arguments.

ClimaCore.Operators.stencil_interior_widthFunction
stencil_interior_width(::Op, args...)

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

Example

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

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

ClimaCore.Operators.x_overlapMethod
x_overlap(target, source)

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

For 1D spaces, this returns the full local weight matrix of the overlap. In 2D, this returns the overlap weights in the first dimension.

ClimaCore.Operators.y_overlapMethod
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.DataLayoutsModule
ClimaCore.DataLayouts

Notation:

  • i,j are horizontal node indices within an element
  • k is the vertical node index within an element
  • f is the field index (1 if field is scalar, >1 if it is a vector field)
  • v is the vertical element index in a stack
  • h is the element stack index

Data layout is specified by the order in which they appear, e.g. IJKFVH indexes the underlying array as [i,j,k,f,v,h]

ClimaCore.DataLayouts.Data1DType
Data1D{S,Ni}

Abstract type for data storage for a 1D field made up of Ni values of type S.

Objects data should define slab(data, h) to return a DataSlab2D{S,Nij} object.

ClimaCore.DataLayouts.Data1DXType
Data1DX{S,Ni}

Abstract type for data storage for a 1D field with extruded columns. The horizontal is made up of Ni values of type S.

Objects data should define slab(data, v, h) to return a DataSlab1D{S,Ni} object, and a column(data, i, h) to return a DataColumn.

ClimaCore.DataLayouts.Data2DType
Data2D{S,Nij}

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

Objects data should define slab(data, h) to return a DataSlab2D{S,Nij} object.

ClimaCore.DataLayouts.Data2DXType
Data2DX{S,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,Nij} object, and a column(data, i, j, h) to return a DataColumn.

ClimaCore.DataLayouts.Data3DType
Data3D{S,Nij,Nk}

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

ClimaCore.DataLayouts.DataColumnType
DataColumn{S}

Abstract type for data storage for a column. Objects data should define a data[k,v], returning a value of type S.

ClimaCore.DataLayouts.DataSlab1DType
DataSlab1D{S,Ni}

Abstract type for data storage for a slab of Ni values of type S. Objects data should define a data[i], returning a value of type S.

ClimaCore.DataLayouts.DataSlab2DType
DataSlab2D{S,Nij}

Abstract type for data storage for a slab of Nij × Nij values of type S. Objects data should define a data[i,j], returning a value of type S.

ClimaCore.DataLayouts.IFType
IF{S, Ni, A} <: DataSlab1D{S, Ni}

Backing DataLayout for 1D spectral element slab data.

Nodal element data (I) are contiguous for each S datatype struct field (F) for a single element slab.

A DataSlab1D view can be returned from other Data1D objects by calling slab(data, idx...).

ClimaCore.DataLayouts.IFHType
IFH{S, Ni, A} <: Data1D{S, Ni}

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

ClimaCore.DataLayouts.IFHMethod
IFH{S,Ni}(ArrayType, nelements)

Construct a IFH 1D Spectral DataLayout given the backing ArrayType, quadrature degrees of freedom Ni , and the number of mesh elements nelements.

ClimaCore.DataLayouts.IH1JH2Method
IH1JH2{S, Nij}(data::AbstractMatrix{S})

Stores a 2D field in a matrix using a column-major format. The primary use is for interpolation to a regular grid for ex. plotting / field output.

ClimaCore.DataLayouts.IJFType
IJF{S, Nij, A} <: DataSlab2D{S, Nij}

Backing DataLayout for 2D spectral element slab data.

Nodal element data (I,J) are contiguous for each S datatype struct field (F) for a single element slab.

A DataSlab2D view can be returned from other Data2D objects by calling slab(data, idx...).

ClimaCore.DataLayouts.IJFHType
IJFH{S, Nij, A} <: Data2D{S, Nij}

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

ClimaCore.DataLayouts.IJFHMethod
IJFH{S,Nij}(ArrayType, nelements)

Construct a IJFH 2D Spectral DataLayout given the backing ArrayType, quadrature degrees of freedom Nij × Nij , and the number of mesh elements nelements.

ClimaCore.DataLayouts.IV1JH2Method
IV1JH2{S, Ni}(data::AbstractMatrix{S})

Stores values from an extruded 1D spectral field in a matrix using a column-major format. The primary use is for interpolation to a regular grid for ex. plotting / field output.

ClimaCore.DataLayouts.VFType
VF{S, A} <: DataColumn{S}

Backing DataLayout for 1D FV column data.

Column level data (V) are contiguous for each S datatype struct field (F).

A DataColumn view can be returned from other Data1DX, Data2DX objects by calling column(data, idx...).

ClimaCore.DataLayouts.VIFHType
VIFH{S, Ni, A} <: Data1DX{S, Ni}

Backing DataLayout for 1D spectral element slab + extruded 1D FV column data.

Column levels (V) are contiguous for every element nodal point (I) for each datatype S struct field (F), for each 1D mesh element slab (H).

ClimaCore.DataLayouts.VIJFHType
VIJFH{S, Nij, A} <: Data2DX{S, Nij}

Backing DataLayout for 2D spectral element slab + extruded 1D FV column data.

Column levels (V) are contiguous for every element nodal point (I, J) for each S datatype struct field (F), for each 2D mesh element slab (H).

ClimaCore.DataLayouts.bypass_constructorMethod
StructArrays.bypass_constructor(T, args)

Create an instance of type T from a tuple of field values args, bypassing possible internal constructors. T should be a concrete type.

ClimaCore.DataLayouts.check_basetypeMethod
check_basetype(::Type{T}, ::Type{S})

Check whether the types T and S have well-defined sizes, and whether an object of type S can be stored in an array with elements of type T.

ClimaCore.DataLayouts.get_structMethod
get_struct(array, S, Val(D), start_index)

Construct an object of type S packed along the D dimension, from the values of array, starting at start_index.

ClimaCore.DataLayouts.is_valid_basetypeMethod
is_valid_basetype(::Type{T}, ::Type{S})

Determines whether an object of type S can be stored in an array with elements of type T by recursively checking whether the non-empty fields of S can be stored in such an array. If S is empty, this is always true.

ClimaCore.DataLayouts.parent_array_typeMethod
parent_array_type(::Type{<:AbstractArray})

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

ClimaCore.DataLayouts.promote_parent_array_typeMethod
promote_parent_array_type(::Type{<:AbstractArray}, ::Type{<:AbstractArray})

Given two parent array types (without any dimensionality information), promote both the element types and the array types themselves.

ClimaCore.DataLayouts.replace_basetypeMethod
replace_basetype(::Type{T}, ::Type{T′}, ::Type{S})

Changes the type parameters of S to produce a new type S′ such that, if is_valid_basetype(T, S) is true, then so is is_valid_basetype(T′, S′).

ClimaCore.DataLayouts.set_struct!Method
set_struct!(array, val::S, Val(D), start_index)

Store an object val of type S packed along the D dimension, into array, starting at start_index.

ClimaCore.Remapping.RemapperMethod

Remapper(targethcoords, targetzcoords, space)

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.

ClimaCore.Remapping.interpolateMethod

interpolate(remapper, field; physical_z = false)

Interpolate the given field as prescribed by remapper.

Keyword arguments

  • physical_z: When true, interpolate to the physical z coordinates, taking into account hypsography. If false (the default), interpolation will be based on the reference z coordinate, i.e. will correspond to constant model levels. NaNs are returned for values that are below the surface.

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(hcoords, zcoords, space)

int1 = interpolate(remapper, field1)
int2 = interpolate(remapper, field2)
ClimaCore.Remapping.interpolate_arrayFunction
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)
Note

Hypsography is not currently handled correctly.

ClimaCore.Remapping.interpolate_columnMethod
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.
ClimaCore.Remapping.interpolate_slab!Method
interpolate_slab!(output_array, field, slab_indices, weights)

Interpolate horizontal field on the given slab_indices using the given interpolation weights.

interpolate_slab! interpolates several values at a fixed z coordinate. For this reason, it requires several slab indices and weights.

ClimaCore.Remapping.interpolate_slab_level!Method
interpolate_slab_level(
                       field::Fields.Field,
                       h::Integer,
                       Is::Tuple,
                       zpts;
                       fill_value = eltype(field)(NaN)
                       )

Vertically interpolate the given field on zpts.

interpolate_slab_level! interpolates several values at a fixed horizontal coordinate.

The field is linearly interpolated across two neighboring vertical elements.

For centered-valued fields, if zcoord is in the top (bottom) half of a top (bottom) element in a column, no interpolation is performed and the value at the cell center is returned. Effectively, this means that the interpolation is first-order accurate across the column, but zeroth-order accurate close to the boundaries.

Return fill_value when the vertical coordinate is negative.

ClimaCore.Remapping.interpolation_weightsFunction
interpolation_weights(horz_mesh, hcoord, quad_points)

Return the weights (tuple of arrays) to interpolate fields onto hcoord on the given mesh and quadrature points.

ClimaCore.Remapping.target_hcoords_pid_bitmaskMethod
target_hcoords_pid_bitmask(target_hcoords, topology, pid)

Return a bitmask for points in target_hcoords in the given (horizontal) topology that belong to the process with process number pid.

This mask can be used to extract the target_hcoords relevant to the given pid. The mask is the same shape and size as the input target_hcoords, which makes it particularly useful.

ClimaCore.Remapping.vertical_indicesFunction
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.