MatrixFields

MatrixFields

Module for defining and manipulating Fields that represent matrices. It adds the BandMatrixRow type, which stores the entries of one row 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,).
• Compatibility with generic data types, e.g., the entries of matrix_field can be iterators 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 is 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) converts each column of matrix_field into a BandedMatrix from BandedMatrices.jl.
• Custom printing, e.g., matrix_field is displayed as the BandedMatrix that corresponds to its first column.

This module also supports sparse block matrices of Fields through the FieldMatrix type (see FieldNameDict), 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.
• Solving linear equations with FieldMatrixSolver, a generalization of ldiv! that is designed to optimize solver performance.

BandMatrixRow{ld}(entries...)

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), with ld = 0.
• BidiagonalMatrixRow(entry_1, entry_2), with ld = -1 + half.
• TridiagonalMatrixRow(entry_1, entry_2, entry_3), with ld = -1.
• QuaddiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4), with ld = -2 + half.
• PentadiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4, entry_5), with ld = -2.

DiagonalMatrixRow{T}
DiagonalMatrixRow(entry)

Alias for BandMatrixRow{0, 1, T}: a row of a BandMatrixRow matrix field with a single entry on the main diagonal (index 0). Fields of DiagonalMatrixRows represent diagonal matrices, and a DiagonalMatrixRow can also appear directly as a scaling entry of a FieldMatrix.

BidiagonalMatrixRow{T}
BidiagonalMatrixRow(entry_1, entry_2)

Alias for BandMatrixRow{-1 + half, 2, T}: a row of a BandMatrixRow matrix field with entries on the diagonals -1/2 and +1/2. This is the row type of matrices that map between cell centers and cell faces, e.g. the matrices of two-point interpolation and difference operators.

TridiagonalMatrixRow{T}
TridiagonalMatrixRow(entry_1, entry_2, entry_3)

Alias for BandMatrixRow{-1, 3, T}: a row of a BandMatrixRow matrix field with entries on the diagonals -1, 0, and +1. This is the row type of square (center-to-center or face-to-face) matrices with a three-point stencil.

QuaddiagonalMatrixRow{T}
QuaddiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4)

Alias for BandMatrixRow{-2 + half, 4, T}: a row of a BandMatrixRow matrix field with entries on the diagonals -3/2, -1/2, +1/2, and +3/2. This is the row type of center-to-face or face-to-center matrices with a four-point stencil, e.g. the matrices of third-order upwinding operators.

PentadiagonalMatrixRow{T}
PentadiagonalMatrixRow(entry_1, entry_2, entry_3, entry_4, entry_5)

Alias for BandMatrixRow{-2, 5, T}: a row of a BandMatrixRow matrix field with entries on the diagonals -2, -1, 0, +1, and +2. This is the row type of square matrices with a five-point stencil, e.g. the product of two BidiagonalMatrixRow matrices and a TridiagonalMatrixRow matrix.

LowerDiagonalMatrixRow{T}
LowerDiagonalMatrixRow(entry)

Alias for BandMatrixRow{-1 + half, 1, T}: a row of a BandMatrixRow matrix field with a single entry on the diagonal -1/2. Together with UpperDiagonalMatrixRow, it is used for the one-sided rows of center-to-face and face-to-center operator matrices, e.g. the boundary rows generated by operator_matrix for interpolation and upwinding operators.

UpperDiagonalMatrixRow{T}
UpperDiagonalMatrixRow(entry)

Alias for BandMatrixRow{half, 1, T}: a row of a BandMatrixRow matrix field with a single entry on the diagonal +1/2. Together with LowerDiagonalMatrixRow, it is used for the one-sided rows of center-to-face and face-to-center operator matrices, e.g. the boundary rows generated by operator_matrix for interpolation and upwinding operators.

MultiplyColumnwiseBandMatrixField()

Operator that multiplies a ColumnwiseBandMatrixField by another Field, i.e., matrix-vector or matrix-matrix multiplication.

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

operator_matrix(op)

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

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

These identities do not hold as stated for gradient and divergence operators. A gradient operator matrix has vector-valued entries and a divergence operator matrix has covector-valued entries, so when ClimaCore itself rewrites a gradient or divergence broadcast into a matrix multiply, it compensates with an adjoint: on the argument for gradients and on the result for divergences. The explicit @. op_matrix() * arg form applies no such compensation, so for a divergence operator it evaluates to adjoint.(@. op(arg)) rather than @. op(arg). When the divergence's result is a scalar (e.g. the divergence of a vector field), the adjoint is a no-op and the identity holds exactly; when the argument is a higher-rank tensor field, the result holds the same components in transposed (row) form, and materializing it into a destination field with the operator's own element type throws a DimensionMismatch.

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 adds an argument to the expression op_matrix.(), from which the space and element type of the operator matrix are inferred.

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) prints a warning and zeros 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) throws an error.

FieldNameDict(keys, entries)
FieldNameDict{T}(key_entry_pairs...)

An AbstractDict with keys of type T, stored as a FieldNameSet{T}, and a tuple of entries. The two-argument constructor takes the key set and the entries; the parametric constructor takes key => entry pairs. T is either FieldName or FieldNamePair, which gives the two aliases:

• FieldMatrix = FieldNameDict{FieldNamePair}, which maps FieldMatrixKeys to ColumnwiseBandMatrixFields, DiagonalMatrixRows, or multiples of LinearAlgebra.I; this is the user-facing alias.
• FieldVectorView = FieldNameDict{FieldName}, which maps FieldVectorKeys to Fields; it is 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 is "lazy" when its entries include AbstractBroadcasted objects that become Fields on materialization. Internal operations produce lazy FieldNameDicts so that a chain of operations materializes once, in a single call to materialize!.

dict[key] returns the entry at key, including entries nested inside a stored entry (e.g. a component of a vector-valued entry); dict[set] returns a FieldNameDict with the entries of dict at the keys in the FieldNameSet set. For a FieldMatrix, one(dict) returns the identity matrix with the diagonal keys of dict.

Broadcasting over FieldNameDicts supports:

• addition, subtraction, and negation,
• multiplication and division by a single value (a Number or a Geometry.SingleValue wrapped in a Ref or a Tuple),
• multiplication, where the first argument is a FieldMatrix,
• inversion of a diagonal FieldMatrix, i.e. one in which every entry is a ColumnwiseBandMatrixField of DiagonalMatrixRows or a multiple of LinearAlgebra.I.

The result of materialize!(dest, bc) covers all keys of dest; entries of dest that are multiples of LinearAlgebra.I are immutable and must equal the corresponding entries of the result.

FieldMatrix

Alias for FieldNameDict{FieldNamePair} (see FieldNameDict): a sparse block matrix that maps (row_name, col_name) pairs of FieldNames, stored as FieldMatrixKeys, to ColumnwiseBandMatrixFields, DiagonalMatrixRows, or multiples of LinearAlgebra.I. Construct one with FieldMatrix(key => entry, ...), and wrap it in a FieldMatrixWithSolver to use it with ldiv!.

FieldVectorView

Alias for FieldNameDict{FieldName} (see FieldNameDict): a dictionary that maps FieldNames to Fields, with keys stored as FieldVectorKeys. It is the form in which a FieldVector participates in operations with a FieldMatrix, such as matrix-vector multiplications or FieldMatrixSolver inversions.

FieldMatrixKeys(values, [name_tree])

Alias for FieldNameSet{Tuple{FieldName, FieldName}}: the key set of a FieldMatrix, i.e. a set of (row_name, col_name) pairs of FieldNames such as ((@name(c.ρ), @name(c.ρ)), (@name(c.ρ), @name(f.u₃))), that serves as the analogue of a KeySet for a FieldNameDict.

FieldVectorKeys(values, [name_tree])

Alias for FieldNameSet{FieldName}: the key set of a FieldVectorView, i.e. a set of FieldNames such as (@name(c.ρ), @name(f.u₃)), that serves as the analogue of a KeySet for a FieldNameDict.

identity_field_matrix(x::Fields.FieldVector)

Construct the FieldMatrix that represents the identity operator on the FieldVector x. It has one diagonal key for every Field in x whose element type is a Geometry.SingleValue (a number or a vector): the entry for a number-valued field is UniformScaling(one(T)), and the entry for a vector-valued field is a DiagonalMatrixRow holding the identity tensor of the field's basis.

Unlike one(matrix), whose keys are the diagonal keys of matrix, the result has an entry for every single-valued field of x, so it covers all the entries needed to solve matrix * x = b for x when matrix is sparse.

field_vector_view(x, [name_tree])

Construct a FieldVectorView whose entries are the Fields in the FieldVector x, keyed by their names in x. name_tree is the FieldNameTree of the keys and defaults to FieldNameTree(x).

concrete_field_vector(vector)

Convert the FieldVectorView vector into a FieldVector, rebuilding the nesting given by the FieldNameTree of its keys.

FieldName(name_chain...)

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

A FieldName behaves like a scalar for broadcasting.

@name(expr)

Construct a FieldName from a chain of getproperty calls. For example:

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

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

FieldNameSet{T}(values, [name_tree])

AbstractSet that contains values of type T, serving as an analogue of a KeySet for a FieldNameDict. There are two aliases of FieldNameSet:

• FieldVectorKeys, for which T is FieldName.
• FieldMatrixKeys, for which T is 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 have no 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). To support this (and to support set complements), FieldNameSet stores a FieldNameTree name_tree, which it uses to infer child values. If name_tree is not specified, it defaults to nothing, which disables some FieldNameSet operations. For binary operations like union or setdiff, only one set needs to specify a name_tree; if both sets specify a name_tree, the name_trees must be identical.

has_field(x, name::FieldName)

Return whether get_field(x, name) is valid, i.e., whether each component of the name chain of name is one of the propertynames of the value selected by the preceding components. Every x has the empty field @name().

get_field(x, name::FieldName)

Return the field of x selected by name, by calling getproperty once for each component of the name chain; e.g., get_field(x, @name(a.b)) is x.a.b. The empty name @name() returns x itself.

is_child_name(child_name::FieldName, parent_name::FieldName)

Return whether the name chain of parent_name is a prefix of the name chain of child_name, so that child_name refers to parent_name or to a field nested inside of it. Every name is a child of itself and of the empty name @name().

append_internal_name(name::FieldName, internal_name::FieldName)

Return the FieldName whose name chain is the concatenation of the name chains of name and internal_name; e.g., append_internal_name(@name(a.b), @name(c)) == @name(a.b.c). This is the inverse of extract_internal_name.

top_level_names(x)

Return a tuple of single-component FieldNames, one for each of the propertynames of x; e.g., (@name(a), @name(b)) for a NamedTuple with keys a and b. The result is an empty tuple when x has no properties.

extract_first(name::FieldName)

Return the first component of the name chain of name, which is either a Symbol or an Integer; e.g., extract_first(@name(a.b.c)) == :a.

drop_first(name::FieldName)

Return the FieldName obtained by removing the first component of the name chain of name; e.g., drop_first(@name(a.b.c)) == @name(b.c).

filtered_names(f, x)

Return a tuple of the FieldNames of all fields of x (including x itself, as @name()) for which f(field) is true, searching the properties of x recursively. The recursion stops at any field that satisfies f, so no returned name is a child of another, and fields without properties that do not satisfy f are omitted.

replace_name_tree(dict::FieldNameDict, name_tree)
replace_name_tree(set::FieldNameSet, name_tree)

Return a copy of dict (or set) whose keys carry the FieldNameTree name_tree in place of their current one, leaving the values unchanged. A FieldMatrix is usually constructed without a name tree, so this is used to attach the name tree of a FieldVector to it before set operations that need to resolve overlapping names, e.g. inside FieldMatrixSolver.

FieldMatrixSolverAlgorithm

Abstract type for algorithms that 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 use the sparsity pattern of A. Subtypes: BlockDiagonalSolve, BlockLowerTriangularSolve, BlockArrowheadSolve, SchurComplementReductionSolve, and LazyFieldMatrixSolverAlgorithm.

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

• field_matrix_solver_cache.
• check_field_matrix_solver.
• run_field_matrix_solver!.

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 passed to this constructor must be similar to the ones passed to field_matrix_solve!, so that the cache is allocated correctly.

FieldMatrixWithSolver(A, b, alg = BlockDiagonalSolve())

Wrapper that combines a FieldMatrix A with a FieldMatrixSolver that can be used to solve the equation A * x = b for x, where x and b are both FieldVectors. Like a LinearAlgebra.Factorization, this wrapper can be passed to ldiv!, whereas a regular FieldMatrix cannot.

By default, the FieldMatrixSolverAlgorithm alg is a BlockDiagonalSolve, so a custom alg must be specified when A is not a block diagonal matrix.

matrix(A::FieldMatrixWithSolver)

Return the FieldMatrix wrapped by the FieldMatrixWithSolver A.

solver(A::FieldMatrixWithSolver)

Return the FieldMatrixSolver of the FieldMatrixWithSolver A, which is nothing if A was created by one(::FieldMatrixWithSolver).

field_matrix_solve!(solver, x, A, b)

Solve the equation A * x = b for x in place using the FieldMatrixSolver solver, and return x.

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.

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, on a CPU, the equation Aₙₙ * xₙ = bₙ is solved using Gaussian elimination (without pivoting), which makes two passes over the data. This is only implemented for tridiagonal and pentadiagonal 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).

On a GPU, tridiagonal systems are solved with the parallel cyclic reduction (PCR) method, which makes better use of GPU parallelism and shared memory. Since this solver launches one thread per row of the matrix, it is only used for systems with at most 512 rows, to stay within CUDA thread block limits. Above that size, the code falls back to the Gaussian elimination method used on the CPU, with degraded performance.

PCR recursively decomposes a tridiagonal system into two systems of half the size, until the systems have size 1 and can be solved directly.

BlockLowerTriangularSolve(names₁...; alg₁ = BlockDiagonalSolve(), alg₂ = BlockDiagonalSolve())

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.

BlockArrowheadSolve(names₁...; alg₂ = BlockDiagonalSolve())

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 also diagonal, so 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₂), can be computed directly with band matrix operations.

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

SchurComplementReductionSolve(names₁...; alg₁ = BlockDiagonalSolve(), 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₁₁) is generally 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 a LazyFieldMatrixSolverAlgorithm, which can evaluate the matrix-vector product (A₂₂ - A₂₁ * inv(A₁₁) * A₁₂) * x₂ without 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.

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.

StationaryIterativeSolve(; P_alg = nothing, n_iters = 1, correlated_solves = false, eigsolve_kwargs = (;), debug = nothing)

A LazyFieldMatrixSolverAlgorithm that solves A * x = b by setting x to some initial value x[0] (usually 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.

Keyword Arguments

• 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 the value of x that was generated during the previous 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: Whether to print debug information with @debug. By default, debug is true when the current logger's minimum level is Debug, e.g., when "error_norm" or "spectral_radius" is in ENV["JULIA_DEBUG"]; setting debug = true without enabling debug logging prints a warning.

Extended help

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 supports 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 not possible on GPUs.

ApproximateBlockArrowheadIterativeSolve(names₁...; P_alg₁ = MainDiagonalPreconditioner(), alg₁ = BlockDiagonalSolve(), alg₂ = BlockDiagonalSolve(), 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 similar to a StationaryIterativeSolve with a BlockArrowheadPreconditioner, but it usually converges more quickly, i.e., the spectral radius of its iteration matrix (I - inv(P) * A) tends to be 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 faster.

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.

PreconditionerAlgorithm

Abstract type for algorithms that approximate a FieldMatrix (or something similar, like a LazySchurComplement) with a preconditioner P for which P * x = b is inexpensive 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. Subtypes: MainDiagonalPreconditioner, BlockDiagonalPreconditioner, BlockArrowheadPreconditioner, BlockArrowheadSchurComplementPreconditioner, WeightedPreconditioner, and CustomPreconditioner.

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

• solver_algorithm.
• lazy_preconditioner.

MainDiagonalPreconditioner()

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

BlockDiagonalPreconditioner()

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

BlockArrowheadPreconditioner(names₁...; P_alg₁ = MainDiagonalPreconditioner(), alg₂ = BlockDiagonalSolve())

A PreconditionerAlgorithm for a 2×2 block matrix:

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

The FieldNames in names₁ correspond to the subscript ₁, while all other FieldNames correspond to the subscript ₂. The preconditioner P is set to the following matrix:

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

The internal preconditioner P₁₁ is generated by the PreconditionerAlgorithm P_alg₁, which is set to a MainDiagonalPreconditioner by default. The Schur complement of P₁₁ in P, A₂₂ - A₂₁ * inv(P₁₁) * A₁₂, is inverted using the FieldMatrixSolverAlgorithm alg₂, which is set to a BlockDiagonalSolve by default.

BlockArrowheadSchurComplementPreconditioner(; P_alg₁ = MainDiagonalPreconditioner(), alg₂ = BlockDiagonalSolve())

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.

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

CustomPreconditioner(M; alg = nothing)

A PreconditionerAlgorithm that sets P to the FieldMatrix M and inverts P using the FieldMatrixSolverAlgorithm alg. The default alg = nothing is only permitted when M is a diagonal matrix.

column_field2array(field)

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

column_field2array_view(field)

Convert a field like column_field2array(field), but return a view that does not copy the data stored in the field.

field2arrays(field)

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

field2arrays_view(field)

Convert a field like field2arrays(field), but call column_field2array_view instead of column_field2array and return a lazy iterator instead of a vector.

scalar_field_matrix(field_matrix::FieldMatrix)

Construct a FieldMatrix whose entries are the scalar blocks of field_matrix: every entry with vector- or tensor-valued blocks is split into one entry per scalar component, keyed by the component names (see get_scalar_keys). Entries that are Fields of scalars are views of the parent arrays of field_matrix.

Examples

e¹² = Geometry.Covariant12Vector(1.6, 0.7)
e₃ = Geometry.Contravariant3Vector(1.0)
e³ = Geometry.Covariant3Vector(1)
ᶜᶜmat3 = fill(TridiagonalMatrixRow(2.0, 3.2, 2.1), center_space)
ᶜᶠmat2 = fill(BidiagonalMatrixRow(4.3, 1.7), center_space)
ᶜᶜmat3_uₕ_scalar = ᶜᶜmat3 .* (e¹²,)
ρχ_unit = (; ρq_liq = 1.0, ρq_ice = 1.0)
ᶜᶠmat2_ρχ_u₃ = map(Base.Fix1(map, Base.Fix2(*, ρχ_unit * e₃')), ᶜᶠmat2)

A = MatrixFields.FieldMatrix(
    (@name(c.ρχ), @name(f.u₃)) => ᶜᶠmat2_ρχ_u₃,
    (@name(c.uₕ), @name(c.sgsʲs.:(1).ρa)) => ᶜᶜmat3_uₕ_scalar,
)

A_scalar = MatrixFields.scalar_field_matrix(A)
keys(A_scalar)
# Output:
# (@name(c.ρχ.ρq_liq), @name(f.u₃.:(1)))
# (@name(c.ρχ.ρq_ice), @name(f.u₃.:(1)))
# (@name(c.uₕ.:(1)), @name(c.sgsʲs.:(1).ρa))
# (@name(c.uₕ.:(2)), @name(c.sgsʲs.:(1).ρa))

band_matrix_info(field)

Return (n_rows, n_cols, matrix_ld, matrix_ud) for a ColumnwiseBandMatrixField: the size of the BandedMatrix that each column of field represents, and the Int indices of its lowest and highest diagonals. The diagonal indices of field's BandMatrixRows, which are half-integers for center-to-face and face-to-center matrices, are shifted by matrix_shape(field) to obtain matrix_ld and matrix_ud. Throws an error if the main diagonal lies outside the band, since BandedMatrices.jl does not support such matrices.