MatrixFields

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
source

Matrix Field Element Type

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

Matrix Field Multiplication

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

source

Operator Matrices

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

source

Linear Solvers

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.

source
ClimaCore.MatrixFields.FieldMatrixSolverType
FieldMatrixSolver(alg, A, b)

Combination of a FieldMatrixSolverAlgorithm and the cache that 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.

source
ClimaCore.MatrixFields.BlockDiagonalSolveType
BlockDiagonalSolve()

A FieldMatrixSolverAlgorithm for a block diagonal matrix A, which solves each block's equation Aᵢᵢ * xᵢ = bᵢ in sequence. The equation for xᵢ is solved as follows:

  • If Aᵢᵢ = λᵢ * I, the equation is solved by setting xᵢ .= inv(λᵢ) .* bᵢ.
  • If Aᵢᵢ = Dᵢ, where Dᵢ is a diagonal matrix, the equation is solved by making a single pass over the data, setting each xᵢ[n] = inv(Dᵢ[n]) * bᵢ[n].
  • If Aᵢᵢ = Lᵢ * Dᵢ * Uᵢ, where Dᵢ is a diagonal matrix and where Lᵢ and Uᵢ are unit lower and upper triangular matrices, respectively, the equation is solved using Gauss-Jordan elimination, which makes two passes over the data. The first pass multiplies both sides of the equation by inv(Lᵢ * Dᵢ), replacing Aᵢᵢ with Uᵢ and bᵢ with Uᵢxᵢ, which is also referred to as putting Aᵢᵢ into "reduced row echelon form". The second pass solves Uᵢ * xᵢ = Uᵢxᵢ for xᵢ using a unit upper triangular matrix solver, which is also referred to as "back substitution". Only tri-diagonal and penta-diagonal matrices Aᵢᵢ are currently supported.
  • The general case of Aᵢᵢ = inv(Pᵢ) * Lᵢ * Uᵢ, where Pᵢ is a row permutation matrix (i.e., LU factorization with partial pivoting), is not currently supported.
source
ClimaCore.MatrixFields.BlockLowerTriangularSolveType
BlockLowerTriangularSolve(names₁...; [alg₁], [alg₂])

A FieldMatrixSolverAlgorithm for a block lower triangular matrix A, which solves for x by executing the following steps:

  1. Partition the entries in A, x, and b into the blocks A₁₁, A₁₂, A₂₁, A₂₂, x₁, x₂, b₁, and b₂, based on the FieldNames in names₁. In this notation, the subscript corresponds to FieldNames that are covered by names₁, while the subscript corresponds to all other FieldNames. A subscript in the first position refers to FieldNames that are used as row indices, while a subscript in the second position refers to column indices. This algorithm requires that the upper triangular block A₁₂ be empty. (Any upper triangular solve can also be expressed as a lower triangular solve by swapping the subscripts and .)
  2. Solve A₁₁ * x₁ = b₁ for x₁ using the algorithm alg₁, which is set to BlockDiagonalSolve() by default.
  3. Solve A₂₂ * x₂ = b₂ - A₂₁ * x₁ for x₂ using the algorithm alg₂, which is set to BlockDiagonalSolve() by default.
source
ClimaCore.MatrixFields.SchurComplementSolveType
SchurComplementSolve(names₁...; [alg₁])

A FieldMatrixSolverAlgorithm for a block matrix A, which solves for x by executing the following steps:

  1. Partition the entries in A, x, and b into the blocks A₁₁, A₁₂, A₂₁, A₂₂, x₁, x₂, b₁, and b₂, based on the FieldNames in names₁. In this notation, the subscript corresponds to FieldNames that are covered by names₁, while the subscript corresponds to all other FieldNames. A subscript in the first position refers to FieldNames that are used as row indices, while a subscript in the second position refers to column indices. This algorithm requires that the block A₂₂ be a diagonal matrix, which allows it to assume that inv(A₂₂) can be computed on the fly.
  2. Solve (A₁₁ - A₁₂ * inv(A₂₂) * A₂₁) * x₁ = b₁ - A₁₂ * inv(A₂₂) * b₂ for x₁ using the algorithm alg₁, which is set to BlockDiagonalSolve() by default. The matrix A₁₁ - A₁₂ * inv(A₂₂) * A₂₁ is called the "Schur complement" of A₂₂ in A.
  3. Set x₂ to inv(A₂₂) * (b₂ - A₂₁ * x₁).
source
ClimaCore.MatrixFields.ApproximateFactorizationSolveType
ApproximateFactorizationSolve(name_pairs₁...; [alg₁], [alg₂])

A FieldMatrixSolverAlgorithm for a block matrix A, which (approximately) solves for x by executing the following steps:

  1. Use the entries in A = M + I = M₁ + M₂ + I to compute A₁ = M₁ + I and A₂ = M₂ + I, based on the pairs of FieldNames in name_pairs₁. In this notation, the subscript refers to pairs of FieldNames that are covered by name_pairs₁, while the subscript refers to all other pairs of FieldNamess. This algorithm approximates the matrix A as the product A₁ * A₂, which introduces an error that scales roughly with the norm of A₁ * A₂ - A = M₁ * M₂. (More precisely, the error introduced by this algorithm is x_exact - x_approx = inv(A) * b - inv(A₁ * A₂) * b.)
  2. Solve A₁ * A₂x = b for A₂x using the algorithm alg₁, which is set to BlockDiagonalSolve() by default.
  3. Solve A₂ * x = A₂x for x using the algorithm alg₂, which is set to BlockDiagonalSolve() by default.
source

Internals

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

source
ClimaCore.MatrixFields.mul_return_typeFunction
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.

source
ClimaCore.MatrixFields.rmul_return_typeFunction
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.

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

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

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

source
ClimaCore.MatrixFields.replace_lazy_operatorFunction
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.

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

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

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

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

source
ClimaCore.MatrixFields.FieldNameDictType
FieldNameDict{T1, T2}(keys, entries)
FieldNameDict{T1, T2}(key_entry_pairs...)

An AbstractDict that contains keys of type T1 and entries of type T2, where the keys are stored as a FieldNameSet{T1}. There are four commonly used 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 are passed to a FieldMatrixSolver)
  • FieldMatrixBroadcasted and FieldVectorViewBroadcasted, which are the same as FieldMatrix and FieldVectorView, except that they can also store unevaluated broadcast expressions; these subtypes are automatically generated when a FieldMatrix or a FieldVectorView is used in a broadcast expression

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 (or FieldMatrixBroadcasted)
  • Inversion, where the argument must be a diagonal FieldMatrix (or FieldMatrixBroadcasted), i.e., one in which every entry is either a ColumnwiseBandMatrixField of DiagonalMatrixRows or a multiple of LinearAlgebra.I
source

Utilities

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

source
ClimaCore.MatrixFields.field2arraysFunction
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.

source