MatrixFields
ClimaCore.MatrixFields
— ModuleMatrixFields
This module adds support for defining and manipulating Field
s that represent matrices. Specifically, it adds the BandMatrixRow
type, which can be used to store the entries of a band matrix. A Field
of BandMatrixRow
s on a FiniteDifferenceSpace
can be interpreted as a band matrix by vertically concatenating the BandMatrixRow
s. Similarly, a Field
of BandMatrixRow
s on an ExtrudedFiniteDifferenceSpace
can be interpreted as a collection of band matrices, one for each column of the Field
. Such Field
s are called ColumnwiseBandMatrixField
s, and this module adds the following functionality for them:
- Constructors, e.g.,
matrix_field = @. BidiagonalMatrixRow(field1, field2)
- Linear combinations, e.g.,
@. 3 * matrix_field1 + matrix_field2 / 3
- Matrix-vector multiplication, e.g.,
@. matrix_field ⋅ field
- Matrix-matrix multiplication, e.g.,
@. matrix_field1 ⋅ matrix_field2
- Compatibility with
LinearAlgebra.I
, e.g.,@. matrix_field = (4I,)
or@. matrix_field - (4I,)
- Integration with
RecursiveApply
, e.g., the entries ofmatrix_field
can beTuple
s orNamedTuple
s instead of single values, which allowsmatrix_field
to represent multiple band matrices at the same time - Integration with
Operators
, e.g., thematrix_field
that gets applied to the argument of anyFiniteDifferenceOperator
op
can be obtained using theFiniteDifferenceOperator
operator_matrix(op)
- Conversions to native array types, e.g.,
field2arrays(matrix_field)
can convert each column ofmatrix_field
into aBandedMatrix
fromBandedMatrices.jl
- Custom printing, e.g.,
matrix_field
gets displayed as aBandedMatrix
, specifically, as theBandedMatrix
that corresponds to its first column
This module also adds support for defining and manipulating sparse block matrices of Field
s. Specifically, it adds the FieldMatrix
type, which is a dictionary that maps pairs of FieldName
s to ColumnwiseBandMatrixField
s or multiples of LinearAlgebra.I
. This comes with the following functionality:
- Addition and subtraction, e.g.,
@. field_matrix1 + field_matrix2
- Matrix-vector multiplication, e.g.,
@. field_matrix * field_vector
- Matrix-matrix multiplication, e.g.,
@. field_matrix1 * field_matrix2
- Integration with
RecursiveApply
, e.g., the entries offield_matrix
can be specified either as matrixField
s ofTuple
s orNamedTuple
s, or as separate matrixField
s of single values - The ability to solve linear equations using
FieldMatrixSolver
, which is a generalization ofldiv!
that is designed to optimize solver performance
Matrix Field Element Type
ClimaCore.MatrixFields.BandMatrixRow
— TypeBandMatrixRow{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)
Matrix Field Multiplication
ClimaCore.MatrixFields.MultiplyColumnwiseBandMatrixField
— TypeMultiplyColumnwiseBandMatrixField()
An operator that multiplies a ColumnwiseBandMatrixField
by another Field
, i.e., matrix-vector or matrix-matrix multiplication. The ⋅
symbol is an alias for MultiplyColumnwiseBandMatrixField()
.
What follows is a derivation of the algorithm used by this operator with single-column Field
s. For Field
s on multiple columns, the same computation is done for each column.
In this derivation, we will use $M_1$ and $M_2$ to denote two ColumnwiseBandMatrixField
s, and we will use $V$ to denote a regular (vector-like) Field
. For both $M_1$ and $M_2$, we will use the array-like index notation $M[row, col]$ to denote $M[row][col-row]$, i.e., the entry in the BandMatrixRow
$M[row]$ located on the diagonal with index $col - row$. We will also use outer_indices
$($space
$)$ to denote the tuple $($left_idx
$($space
$),$right_idx
$($space
$))$.
1. Matrix-Vector Multiplication
From the definition of matrix-vector multiplication,
\[(M_1 ⋅ V)[i] = \sum_k M_1[i, k] * V[k].\]
To establish bounds on the values of $k$, let us define the following values:
- $li_1, ri_1 ={}$
outer_indices
$($column_axes
$(M_1))$ - $ld_1, ud_1 ={}$
outer_diagonals
$($eltype
$(M_1))$
Since $M_1[i, k]$ is only well-defined if $k$ is a valid column index and $k - i$ is a valid diagonal index, we know that
\[li_1 \leq k \leq ri_1 \quad \text{and} \quad ld_1 \leq k - i \leq ud_1.\]
Combining these into a single inequality gives us
\[\text{max}(li_1, i + ld_1) \leq k \leq \text{min}(ri_1, i + ud_1).\]
So, we can rewrite the expression for $(M_1 ⋅ V)[i]$ as
\[(M_1 ⋅ V)[i] = \sum_{k\ =\ \text{max}(li_1, i + ld_1)}^{\text{min}(ri_1, i + ud_1)} M_1[i, k] * V[k].\]
If we replace the variable $k$ with $d = k - i$ and switch from array-like indexing to Field
indexing, we find that
\[(M_1 ⋅ V)[i] = \sum_{d\ =\ \text{max}(li_1 - i, ld_1)}^{\text{min}(ri_1 - i, ud_1)} M_1[i][d] * V[i + d].\]
1.1 Interior vs. Boundary Indices
Now, suppose that the row index $i$ is such that
\[li_1 - ld_1 \leq i \leq ri_1 - ud_1.\]
If this is the case, then the bounds on $d$ can be simplified to
\[\text{max}(li_1 - i, ld_1) = ld_1 \quad \text{and} \quad \text{min}(ri_1 - i, ud_1) = ud_1.\]
The expression for $(M_1 ⋅ V)[i]$ then becomes
\[(M_1 ⋅ V)[i] = \sum_{d = ld_1}^{ud_1} M_1[i][d] * V[i + d].\]
The values of $i$ in this range are considered to be in the "interior" of the operator, while those not in this range (for which we cannot make the above simplification) are considered to be on the "boundary".
2. Matrix-Matrix Multiplication
From the definition of matrix-matrix multiplication,
\[(M_1 ⋅ M_2)[i, j] = \sum_k M_1[i, k] * M_2[k, j].\]
To establish bounds on the values of $j$ and $k$, let us define the following values:
- $li_1, ri_1 ={}$
outer_indices
$($column_axes
$(M_1))$ - $ld_1, ud_1 ={}$
outer_diagonals
$($eltype
$(M_1))$ - $li_2, ri_2 ={}$
outer_indices
$($column_axes
$(M_2))$ - $ld_2, ud_2 ={}$
outer_diagonals
$($eltype
$(M_2))$
In addition, let $ld_{prod}$ and $ud_{prod}$ denote the outer diagonal indices of the product matrix $M_1 ⋅ M_2$. We will derive the values of $ld_{prod}$ and $ud_{prod}$ in the last section.
Since $M_1[i, k]$ is only well-defined if $k$ is a valid column index and $k - i$ is a valid diagonal index, we know that
\[li_1 \leq k \leq ri_1 \quad \text{and} \quad ld_1 \leq k - i \leq ud_1.\]
Since $M_2[k, j]$ is only well-defined if $j$ is a valid column index and $j - k$ is a valid diagonal index, we also know that
\[li_2 \leq j \leq ri_2 \quad \text{and} \quad ld_2 \leq j - k \leq ud_2.\]
Finally, $(M_1 ⋅ M_2)[i, j]$ is only well-defined if $j - i$ is a valid diagonal index, so
\[ld_{prod} \leq j - i \leq ud_{prod}.\]
These inequalities can be combined to obtain
\[\begin{gather*} \text{max}(li_2, i + ld_{prod}) \leq j \leq \text{min}(ri_2, i + ud_{prod}) \\ \text{and} \\ \text{max}(li_1, i + ld_1, j - ud_2) \leq k \leq \text{min}(ri_1, i + ud_1, j - ld_2). \end{gather*}\]
So, we can rewrite the expression for $(M_1 ⋅ M_2)[i, j]$ as
\[\begin{gather*} (M_1 ⋅ M_2)[i, j] = \sum_{ k\ =\ \text{max}(li_1, i + ld_1, j - ud_2) }^{\text{min}(ri_1, i + ud_1, j - ld_2)} M_1[i, k] * M_2[k, j], \text{ where} \\[0.5em] \text{max}(li_2, i + ld_{prod}) \leq j \leq \text{min}(ri_2, i + ud_{prod}). \end{gather*}\]
If we replace the variable $k$ with $d = k - i$, replace the variable $j$ with $d_{prod} = j - i$, and switch from array-like indexing to Field
indexing, we find that
\[\begin{gather*} (M_1 ⋅ M_2)[i][d_{prod}] = \sum_{ d\ =\ \text{max}(li_1 - i, ld_1, d_{prod} - ud_2) }^{\text{min}(ri_1 - i, ud_1, d_{prod} - ld_2)} M_1[i][d] * M_2[i + d][d_{prod} - d], \text{ where} \\[0.5em] \text{max}(li_2 - i, ld_{prod}) \leq d_{prod} \leq \text{min}(ri_2 - i, ud_{prod}). \end{gather*}\]
2.1 Interior vs. Boundary Indices
Now, suppose that the row index $i$ is such that
\[\text{max}(li_1 - ld_1, li_2 - ld_{prod}) \leq i \leq \text{min}(ri_1 - ud_1, ri_2 - ud_{prod}).\]
If this is the case, then the bounds on $d_{prod}$ can be simplified to
\[\text{max}(li_2 - i, ld_{prod}) = ld_{prod} \quad \text{and} \quad \text{min}(ri_2 - i, ud_{prod}) = ud_{prod}.\]
Similarly, the bounds on $d$ can be simplified using the fact that
\[\text{max}(li_1 - i, ld_1) = ld_1 \quad \text{and} \quad \text{min}(ri_1 - i, ud_1) = ud_1.\]
The expression for $(M_1 ⋅ M_2)[i][d_{prod}]$ then becomes
\[\begin{gather*} (M_1 ⋅ M_2)[i][d_{prod}] = \sum_{ d\ =\ \text{max}(ld_1, d_{prod} - ud_2) }^{\text{min}(ud_1, d_{prod} - ld_2)} M_1[i][d] * M_2[i + d][d_{prod} - d], \text{ where} \\[0.5em] ld_{prod} \leq d_{prod} \leq ud_{prod}. \end{gather*}\]
The values of $i$ in this range are considered to be in the "interior" of the operator, while those not in this range (for which we cannot make these simplifications) are considered to be on the "boundary".
2.2 $ld_{prod}$ and $ud_{prod}$
We only need to compute $(M_1 ⋅ M_2)[i][d_{prod}]$ for values of $d_{prod}$ that correspond to a nonempty sum in the interior, i.e, those for which
\[\text{max}(ld_1, d_{prod} - ud_2) \leq \text{min}(ud_1, d_{prod} - ld_2).\]
This can be broken down into the four inequalities
\[ld_1 \leq ud_1, \qquad ld_1 \leq d_{prod} - ld_2, \qquad d_{prod} - ud_2 \leq ud_1, \quad \text{and} \quad d_{prod} - ud_2 \leq d_{prod} - ld_2.\]
By definition, $ld_1 \leq ud_1$ and $ld_2 \leq ud_2$, so the first and last inequality are always true. Rearranging the remaining two inequalities tells us that
\[ld_1 + ld_2 \leq d_{prod} \leq ud_1 + ud_2.\]
In other words, the outer diagonal indices of $M_1 ⋅ M_2$ are
\[ld_{prod} = ld_1 + ld_2 \quad \text{and} \quad ud_{prod} = ud_1 + ud_2.\]
This means that we can express the bounds on the interior values of $i$ as
\[\text{max}(li_1, li_2 - ld_2) - ld_1 \leq i \leq \text{min}(ri_1, ri_2 - ud_2) - ud_1.\]
Operator Matrices
ClimaCore.MatrixFields.operator_matrix
— Functionoperator_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)
.
Linear Solvers
ClimaCore.MatrixFields.FieldMatrixSolverAlgorithm
— TypeFieldMatrixSolverAlgorithm
Description of how to solve an equation of the form A * x = b
for x
, where A
is a FieldMatrix
and where x
and b
are both FieldVector
s. Different algorithms can be nested inside each other, enabling the construction of specialized linear solvers that fully utilize the sparsity pattern of A
.
ClimaCore.MatrixFields.FieldMatrixSolver
— TypeFieldMatrixSolver(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.
ClimaCore.MatrixFields.field_matrix_solve!
— Functionfield_matrix_solve!(solver, x, A, b)
Solves the equation A * x = b
for x
using the given FieldMatrixSolver
.
ClimaCore.MatrixFields.BlockDiagonalSolve
— TypeBlockDiagonalSolve()
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 settingxᵢ .= inv(λᵢ) .* bᵢ
. - If
Aᵢᵢ = Dᵢ
, whereDᵢ
is a diagonal matrix, the equation is solved by making a single pass over the data, setting eachxᵢ[n] = inv(Dᵢ[n]) * bᵢ[n]
. - If
Aᵢᵢ = Lᵢ * Dᵢ * Uᵢ
, whereDᵢ
is a diagonal matrix and whereLᵢ
andUᵢ
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 byinv(Lᵢ * Dᵢ)
, replacingAᵢᵢ
withUᵢ
andbᵢ
withUᵢxᵢ
, which is also referred to as puttingAᵢᵢ
into "reduced row echelon form". The second pass solvesUᵢ * xᵢ = Uᵢxᵢ
forxᵢ
using a unit upper triangular matrix solver, which is also referred to as "back substitution". Only tri-diagonal and penta-diagonal matricesAᵢᵢ
are currently supported. - The general case of
Aᵢᵢ = inv(Pᵢ) * Lᵢ * Uᵢ
, wherePᵢ
is a row permutation matrix (i.e., LU factorization with partial pivoting), is not currently supported.
ClimaCore.MatrixFields.BlockLowerTriangularSolve
— TypeBlockLowerTriangularSolve(names₁...; [alg₁], [alg₂])
A FieldMatrixSolverAlgorithm
for a block lower triangular matrix A
, which solves for x
by executing the following steps:
- Partition the entries in
A
,x
, andb
into the blocksA₁₁
,A₁₂
,A₂₁
,A₂₂
,x₁
,x₂
,b₁
, andb₂
, based on theFieldName
s innames₁
. In this notation, the subscript₁
corresponds toFieldName
s that are covered bynames₁
, while the subscript₂
corresponds to all otherFieldNames
. A subscript in the first position refers toFieldName
s that are used as row indices, while a subscript in the second position refers to column indices. This algorithm requires that the upper triangular blockA₁₂
be empty. (Any upper triangular solve can also be expressed as a lower triangular solve by swapping the subscripts₁
and₂
.) - Solve
A₁₁ * x₁ = b₁
forx₁
using the algorithmalg₁
, which is set toBlockDiagonalSolve()
by default. - Solve
A₂₂ * x₂ = b₂ - A₂₁ * x₁
forx₂
using the algorithmalg₂
, which is set toBlockDiagonalSolve()
by default.
ClimaCore.MatrixFields.SchurComplementSolve
— TypeSchurComplementSolve(names₁...; [alg₁])
A FieldMatrixSolverAlgorithm
for a block matrix A
, which solves for x
by executing the following steps:
- Partition the entries in
A
,x
, andb
into the blocksA₁₁
,A₁₂
,A₂₁
,A₂₂
,x₁
,x₂
,b₁
, andb₂
, based on theFieldName
s innames₁
. In this notation, the subscript₁
corresponds toFieldName
s that are covered bynames₁
, while the subscript₂
corresponds to all otherFieldNames
. A subscript in the first position refers toFieldName
s that are used as row indices, while a subscript in the second position refers to column indices. This algorithm requires that the blockA₂₂
be a diagonal matrix, which allows it to assume thatinv(A₂₂)
can be computed on the fly. - Solve
(A₁₁ - A₁₂ * inv(A₂₂) * A₂₁) * x₁ = b₁ - A₁₂ * inv(A₂₂) * b₂
forx₁
using the algorithmalg₁
, which is set toBlockDiagonalSolve()
by default. The matrixA₁₁ - A₁₂ * inv(A₂₂) * A₂₁
is called the "Schur complement" ofA₂₂
inA
. - Set
x₂
toinv(A₂₂) * (b₂ - A₂₁ * x₁)
.
ClimaCore.MatrixFields.ApproximateFactorizationSolve
— TypeApproximateFactorizationSolve(name_pairs₁...; [alg₁], [alg₂])
A FieldMatrixSolverAlgorithm
for a block matrix A
, which (approximately) solves for x
by executing the following steps:
- Use the entries in
A = M + I = M₁ + M₂ + I
to computeA₁ = M₁ + I
andA₂ = M₂ + I
, based on the pairs ofFieldName
s inname_pairs₁
. In this notation, the subscript₁
refers to pairs ofFieldName
s that are covered byname_pairs₁
, while the subscript₂
refers to all other pairs ofFieldNames
s. This algorithm approximates the matrixA
as the productA₁ * A₂
, which introduces an error that scales roughly with the norm ofA₁ * A₂ - A = M₁ * M₂
. (More precisely, the error introduced by this algorithm isx_exact - x_approx = inv(A) * b - inv(A₁ * A₂) * b
.) - Solve
A₁ * A₂x = b
forA₂x
using the algorithmalg₁
, which is set toBlockDiagonalSolve()
by default. - Solve
A₂ * x = A₂x
forx
using the algorithmalg₂
, which is set toBlockDiagonalSolve()
by default.
Internals
ClimaCore.MatrixFields.outer_diagonals
— Functionouter_diagonals(::Type{<:BandMatrixRow})
Gets the indices of the lower and upper diagonals, ld
and ud
, of the given subtype of BandMatrixRow
.
ClimaCore.MatrixFields.band_matrix_row_type
— Functionband_matrix_row_type(ld, ud, T)
A shorthand for getting the subtype of BandMatrixRow
that has entries of type T
on the diagonals with indices in the range ld:ud
.
ClimaCore.MatrixFields.mul_with_projection
— Functionmul_with_projection(x, y, lg)
Similar to x * y
, except that this version automatically projects y
to avoid DimensionMismatch
errors for AxisTensor
s. For example, if x
is a covector along the Covariant3Axis
(e.g., Covariant3Vector(1)'
), then y
will be projected onto the Contravariant3Axis
. In general, the first axis of y
will be projected onto the dual of the last axis of x
.
ClimaCore.MatrixFields.rmul_with_projection
— Functionrmul_with_projection(x, y, lg)
Similar to rmul(x, y)
, except that this version calls mul_with_projection
instead of *
.
ClimaCore.MatrixFields.mul_return_type
— Functionmul_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.rmul_return_type
— Functionrmul_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.matrix_shape
— Functionmatrix_shape(matrix_field, [matrix_space])
Returns either Square()
, FaceToCenter()
, or CenterToFace()
, depending on whether the diagonal indices of matrix_field
are Int
s or PlusHalf
s and whether matrix_space
is on cell centers or cell faces. By default, matrix_space
is set to axes(matrix_field)
.
ClimaCore.MatrixFields.column_axes
— Functioncolumn_axes(matrix_field, [matrix_space])
Returns the space that corresponds to the columns of matrix_field
, i.e., the axes
of the Field
s by which matrix_field
can be multiplied. The matrix_space
, on the other hand, is the space that corresponds to the rows of matrix_field
. By default, matrix_space
is set to axes(matrix_field)
.
ClimaCore.MatrixFields.AbstractLazyOperator
— TypeAbstractLazyOperator
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.replace_lazy_operator
— Functionreplace_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.FieldName
— TypeFieldName(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.@name
— Macro@name(expr)
Shorthand for constructing a FieldName
. Some examples include
name = @name()
, in which caseget_field(x, name)
returnsx
name = @name(a)
, in which caseget_field(x, name)
returnsx.a
name = @name(a.b.c)
, in which caseget_field(x, name)
returnsx.a.b.c
name = @name(a.b.c.:(1).d)
, in which caseget_field(x, name)
returnsx.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.MatrixFields.FieldNameTree
— TypeFieldNameTree(x)
Tree of FieldName
s that can be used to access x
with get_field(x, name)
. Check whether a name
is valid by calling is_valid_name(name, tree)
, and extract the children of name
by calling child_names(name, tree)
.
ClimaCore.MatrixFields.FieldNameSet
— TypeFieldNameSet{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 whichT
is set toFieldName
FieldMatrixKeys
, for whichT
is set toTuple{FieldName, FieldName}
; each tuple of typeT
represents a pair of row-column indices
Since FieldName
s are singleton types, the result of almost any FieldNameSet
operation can be inferred during compilation. So, with the exception of map
, foreach
, and set_string
, functions of FieldNameSet
s do not have any performance cost at runtime (as long as their arguments are inferrable).
Unlike other AbstractSet
s, FieldNameSet
has special behavior for overlapping values. For example, the FieldName
s @name(a.b)
and @name(a.b.c)
overlap, so any set operation needs to first decompose @name(a.b)
into its child values before combining it with @name(a.b.c)
. In order to support this (and also to support the ability to compute set complements), FieldNameSet
stores a FieldNameTree
name_tree
, which it uses to infer child values. If name_tree
is not specified, it gets set to nothing
by default, which causes some FieldNameSet
operations to become disabled. For binary operations like union
or setdiff
, only one set needs to specify a name_tree
; if two sets both specify a name_tree
, the name_tree
s must be identical.
ClimaCore.MatrixFields.FieldNameDict
— TypeFieldNameDict{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 ofFieldMatrixKeys
to eitherColumnwiseBandMatrixField
s or multiples ofLinearAlgebra.I
; this is the only user-facing subtype ofFieldNameDict
FieldVectorView
, which maps a set ofFieldVectorKeys
toField
s; this subtype is automatically generated when aFieldVector
is used in the same operation as aFieldMatrix
(e.g., when both appear in the same broadcast expression or are passed to aFieldMatrixSolver
)FieldMatrixBroadcasted
andFieldVectorViewBroadcasted
, which are the same asFieldMatrix
andFieldVectorView
, except that they can also store unevaluated broadcast expressions; these subtypes are automatically generated when aFieldMatrix
or aFieldVectorView
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 FieldNameDict
s, the following operations are supported:
- Addition and subtraction
- Multiplication, where the first argument must be a
FieldMatrix
(orFieldMatrixBroadcasted
) - Inversion, where the argument must be a diagonal
FieldMatrix
(orFieldMatrixBroadcasted
), i.e., one in which every entry is either aColumnwiseBandMatrixField
ofDiagonalMatrixRow
s or a multiple ofLinearAlgebra.I
ClimaCore.MatrixFields.field_vector_view
— Functionfield_vector_view(x)
Constructs a FieldVectorView
that contains all the top-level Field
s in the FieldVector
x
.
Utilities
ClimaCore.MatrixFields.column_field2array
— Functioncolumn_field2array(field)
Converts a field defined on a FiniteDifferenceSpace
into either a Vector
or a BandedMatrix
, depending on whether the elements of the field are single values or BandMatrixRow
s. This involves copying the data stored in the field. Because BandedMatrix
does not currently support operations with CuArray
s, all GPU data is copied to the CPU.
ClimaCore.MatrixFields.column_field2array_view
— Functioncolumn_field2array_view(field)
Similar to column_field2array(field)
, except that this version avoids copying the data stored in the field.
ClimaCore.MatrixFields.field2arrays
— Functionfield2arrays(field)
Converts a field defined on a FiniteDifferenceSpace
or on an ExtrudedFiniteDifferenceSpace
into a collection of arrays, each of which corresponds to a column of the field. This is done by calling column_field2array
on each of the field's columns.
ClimaCore.MatrixFields.field2arrays_view
— Functionfield2arrays_view(field)
Similar to field2arrays(field)
, except that this version calls column_field2array_view
instead of column_field2array
.