# System Solvers

## Non-linear solvers

ClimateMachine.SystemSolvers.JacobianActionType

mutable struct JacobianAction{FT, AT} rhs! ϵ::FT Q::AT Qdq::AT Fq::AT Fqdq::AT end

Solve for Frhs = F(q), the Jacobian action is computed

∂F(Q)      F(Q + eΔQ) - F(Q)
---- ΔQ ≈ -------------------
∂Q                e

...

Arguments

• rhs! : nonlinear operator F(Q)
• ϵ::FT : ϵ used for finite difference, e = e(Q, ϵ)
• Q::AT : cache for Q
• Qdq::AT : container for Q + ϵΔQ
• Fq::AT : cache for F(Q)
• Fqdq::AT : container for F(Q + ϵΔQ)

...

source

## Linear solvers

### Generalized Conjugate Residual Method

ClimateMachine.SystemSolvers.GeneralizedConjugateResidualType
GeneralizedConjugateResidual(K, Q; rtol, atol)

Conjugate Residual

This is an object for solving linear systems using an iterative Krylov method. The constructor parameter K is the number of steps after which the algorithm is restarted (if it has not converged), Q is a reference state used only to allocate the solver internal state, and tolerance specifies the convergence criterion based on the relative residual norm. The amount of memory required by the solver state is roughly (2K + 2) * size(Q). This object is intended to be passed to the linearsolve! command.

This uses the restarted Generalized Conjugate Residual method of Eisenstat (1983).

References

source

### Generalized Minimal Residual Method

ClimateMachine.SystemSolvers.GeneralizedMinimalResidualType
GeneralizedMinimalResidual(Q; M, rtol, atol)

GMRES

This is an object for solving linear systems using an iterative Krylov method. The constructor parameter M is the number of steps after which the algorithm is restarted (if it has not converged), Q is a reference state used only to allocate the solver internal state, and rtol specifies the convergence criterion based on the relative residual norm. The amount of memory required for the solver state is roughly (M + 1) * size(Q). This object is intended to be passed to the linearsolve! command.

This uses the restarted Generalized Minimal Residual method of Saad and Schultz (1986).

References

source

### Batched Generalized Minimal Residual Method

ClimateMachine.SystemSolvers.BatchedGeneralizedMinimalResidualType
BatchedGeneralizedMinimalResidual(
Q,
dofperbatch,
Nbatch;
M = min(20, length(Q)),
rtol = √eps(eltype(AT)),
atol = eps(eltype(AT)),
forward_reshape = size(Q),
forward_permute = Tuple(1:length(size(Q))),
)

BGMRES

This is an object for solving batched linear systems using the GMRES algorithm. The constructor parameter M is the number of steps after which the algorithm is restarted (if it has not converged), Q is a reference state used only to allocate the solver internal state, dofperbatch is the size of each batched system (assumed to be the same throughout), Nbatch is the total number of independent linear systems, and rtol specifies the convergence criterion based on the relative residual norm (max across all batched systems). The argument forward_reshape is a tuple of integers denoting the reshaping (if required) of the solution vectors for batching the Arnoldi routines. The argument forward_permute describes precisely which indices of the array Q to permute. This object is intended to be passed to the linearsolve! command.

This uses a batched-version of the restarted Generalized Minimal Residual method of Saad and Schultz (1986).

Note

Eventually, we'll want to do something like this:

i = @index(Global)
linearoperator!(Q[:, :, :, i], args...)

This will help stop the need for constantly reshaping the work arrays. It would also potentially save us some memory.

source

Q::AT;
rtol = eps(eltype(Q)),
atol = eps(eltype(Q)),
max_iter = length(Q),
dims = :,
reshape_tuple = size(Q),
) where {AT}

Description

• Outer constructor for the ConjugateGradient struct

Arguments

• Q:(array). The kind of object that linearoperator! acts on.

Keyword Arguments

• rtol: (float). relative tolerance
• atol: (float). absolute tolerance
• dims: (tuple or : ). the dimensions to reduce over
• reshape_tuple: (tuple). the dimensions that the conjugate gradient solver operators over

Comment

• The reshape tuple is necessary in case the linearoperator! is defined over vectors of a different size as compared to what plays nicely with the dimension reduction in the ConjugateGradient. It also allows the user to define preconditioners over arrays that are more convenienently shaped.

Return

source
Q::MPIStateArray;
rtol = eps(eltype(Q)),
atol = eps(eltype(Q)),
max_iter = length(Q),
dims = :,
reshape_tuple = size(Q),
)

Description

Outer constructor for the ConjugateGradient struct with MPIStateArrays. THIS IS A HACK DUE TO RESHAPE FUNCTIONALITY ON MPISTATEARRAYS.

Arguments

• Q:(array). The kind of object that linearoperator! acts on.

Keyword Arguments

• rtol: (float). relative tolerance
• atol: (float). absolute tolerance
• dims: (tuple or : ). the dimensions to reduce over
• reshape_tuple: (tuple). the dimensions that the conjugate gradient solver operators over

Comment

• The reshape tuple is necessary in case the linearoperator! is defined over vectors of a different size as compared to what plays nicely with the dimension reduction in the ConjugateGradient. It also allows the user to define preconditioners over arrays that are more convenienently shaped.

Return

source
ClimateMachine.SystemSolvers.initialize!Function
initialize!(
linearoperator!,
Q,
Qrhs,
args...,
)

Description

• This function initializes the iterative solver. It is called as part of the AbstractIterativeSystemSolver routine. SEE CODEREF for documentation on AbstractIterativeSystemSolver

Arguments

• linearoperator!: (function). This applies the predefined linear operator on an array. Applies a linear operator to object "y" and overwrites object "z". The function argument i s linearoperator!(z,y, args...) and it returns nothing.
• Q: (array). This is an object that linearoperator! outputs
• Qrhs: (array). This is an object that linearoperator! acts on
• solver: (struct). This is a scruct for dispatch, in this case for ColumnwisePreconditionedConjugateGradient
• args...: (arbitrary). This is optional arguments that can be passed into linearoperator! function.

Keyword Arguments

• There are no keyword arguments

Return

• converged: (bool). A boolean to say whether or not the iterative solver has converged.
• threshold: (float). The value of the residual for the first timestep

Comment

• This function does nothing for conjugate gradient
source

JacobianFreeNewtonKrylovSolver initialize the residual

source
ClimateMachine.SystemSolvers.doiteration!Function
doiteration!(
linearoperator!,
preconditioner,
Q,
Qrhs,
threshold,
args...;
applyPC! = (x, y) -> x .= y,
)

Description

• This function enacts the iterative solver. It is called as part of the AbstractIterativeSystemSolver routine. SEE CODEREF for documentation on AbstractIterativeSystemSolver

Arguments

• linearoperator!: (function). This applies the predefined linear operator on an array. Applies a linear operator to object "y" and overwrites object "z". It is a function with arguments linearoperator!(z,y, args...), where "z" gets overwritten by "y" and "args..." are additional arguments passed to the linear operator. The linear operator is assumed to return nothing.
• Q: (array). This is an object that linearoperator! overwrites
• Qrhs: (array). This is an object that linearoperator! acts on. This is the rhs to the linear system
• solver: (struct). This is a scruct for dispatch, in this case for ConjugateGradient
• threshold: (float). Either an absolute or relative tolerance
• applyPC!: (function). Applies a preconditioner to objecy "y" and overwrites object "z". applyPC!(z,y)
• args...: (arbitrary). This is necessary for the linearoperator! function which has a signature linearoperator!(b, x, args....)

Keyword Arguments

• There are no keyword arguments

Return

• converged: (bool). A boolean to say whether or not the iterative solver has converged.
• iteration: (int). Iteration number for the iterative solver
• threshold: (float). The value of the residual for the first timestep

Comment

• This function does conjugate gradient
source
doiteration!(
linearoperator!,
preconditioner,
Q::MPIStateArray,
Qrhs::MPIStateArray,
threshold,
args...;
applyPC! = (x, y) -> x .= y,
)

Description

This function enacts the iterative solver. It is called as part of the AbstractIterativeSystemSolver routine. SEE CODEREF for documentation on AbstractIterativeSystemSolver. THIS IS A HACK TO WORK WITH MPISTATEARRAYS. THE ISSUE IS WITH RESHAPE.

Arguments

• linearoperator!: (function). This applies the predefined linear operator on an array. Applies a linear operator to object "y" and overwrites object "z". It is a function with arguments linearoperator!(z,y, args...), where "z" gets overwritten by "y" and "args..." are additional arguments passed to the linear operator. The linear operator is assumed to return nothing.
• Q: (array). This is an object that linearoperator! overwrites
• Qrhs: (array). This is an object that linearoperator! acts on. This is the rhs to the linear system
• solver: (struct). This is a scruct for dispatch, in this case for ConjugateGradient
• threshold: (float). Either an absolute or relative tolerance
• applyPC!: (function). Applies a preconditioner to objecy "y" and overwrites object "z". applyPC!(z,y)
• args...: (arbitrary). This is necessary for the linearoperator! function which has a signature linearoperator!(b, x, args....)

Keyword Arguments

• There are no keyword arguments

Return

• converged: (bool). A boolean to say whether or not the iterative solver has converged.
• iteration: (int). Iteration number for the iterative solver
• threshold: (float). The value of the residual for the first timestep

Comment

• This function does conjugate gradient
source

### LU Decomposition

ClimateMachine.SystemSolvers.ManyColumnLUType
ManyColumnLU()

This solver is used for systems that are block diagonal where each block is associated with a column of the mesh. The systems are solved using a non-pivoted LU factorization.

source
ClimateMachine.SystemSolvers.SingleColumnLUType
SingleColumnLU()

This solver is used for systems that are block diagonal where each block is associated with a column of the mesh. Moreover, each block is assumed to be the same. The systems are solved using a non-pivoted LU factorization.

source

## Preconditioners

ClimateMachine.SystemSolvers.ColumnwiseLUPreconditionerType

mutable struct ColumnwiseLUPreconditioner{AT} A::DGColumnBandedMatrix Q::AT PQ::AT counter::Int update_freq::Int end

...

Arguments

• A: the lu factor of the precondition (approximated Jacobian), in the DGColumnBandedMatrix format
• Q: MPIArray container, used to update A
• PQ: MPIArray container, used to update A
• counter: count the number of Newton, when counter > update_freq or counter < 0, update precondition
• update_freq: preconditioner update frequency

...

source

## Shared components

ClimateMachine.SystemSolvers.linearsolve!Function
linearsolve!(linearoperator!, solver::AbstractIterativeSystemSolver, Q, Qrhs, args...)

Solves a linear problem defined by the linearoperator! function and the state Qrhs, i.e,

$$$L(Q) = Q_{rhs}$$$

using the solver and the initial guess Q. After the call Q contains the solution. The arguments args is passed to linearoperator! when it is called.

source
ClimateMachine.SystemSolvers.preconditioner_update!Function

Do nothing, when there is no preconditioner, preconditioner = Nothing

source

update the DGColumnBandedMatrix by the finite difference approximation ...

Arguments

• op: operator used to compute the finte difference information
• dg: the DG model, use only the grid information

...

source