System Solvers
Non-linear solvers
ClimateMachine.SystemSolvers.LSOnly
— TypeLSOnly
Only applies the linear solver (no Newton solver)
ClimateMachine.SystemSolvers.JacobianAction
— Typemutable 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 QQdq::AT
: container for Q + ϵΔQFq::AT
: cache for F(Q)Fqdq::AT
: container for F(Q + ϵΔQ)
...
ClimateMachine.SystemSolvers.JacobianFreeNewtonKrylovSolver
— TypeSolve for Frhs = F(Q), by finite difference
∂F(Q) F(Q + eΔQ) - F(Q)
---- ΔQ ≈ -------------------
∂Q e
Q^n+1 = Q^n - dF/dQ(Q^{n})⁻¹ (F(Q^n) - Frhs)
set ΔQ = F(Q^n) - Frhs
Linear solvers
Generalized Conjugate Residual Method
ClimateMachine.SystemSolvers.GeneralizedConjugateResidual
— TypeGeneralizedConjugateResidual(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
Generalized Minimal Residual Method
ClimateMachine.SystemSolvers.GeneralizedMinimalResidual
— TypeGeneralizedMinimalResidual(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
Batched Generalized Minimal Residual Method
ClimateMachine.SystemSolvers.BatchedGeneralizedMinimalResidual
— TypeBatchedGeneralizedMinimalResidual(
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.
Conjugate Gradient Solver Method
ClimateMachine.SystemSolvers.ConjugateGradient
— TypeConjugateGradient(
Q::AT;
rtol = eps(eltype(Q)),
atol = eps(eltype(Q)),
max_iter = length(Q),
dims = :,
reshape_tuple = size(Q),
) where {AT}
ConjugateGradient
Description
- Outer constructor for the ConjugateGradient struct
Arguments
Q
:(array). The kind of object that linearoperator! acts on.
Keyword Arguments
rtol
: (float). relative toleranceatol
: (float). absolute tolerancedims
: (tuple or : ). the dimensions to reduce overreshape_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
- ConjugateGradient struct
ConjugateGradient(
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 toleranceatol
: (float). absolute tolerancedims
: (tuple or : ). the dimensions to reduce overreshape_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
- ConjugateGradient struct
ClimateMachine.SystemSolvers.initialize!
— Functioninitialize!(
linearoperator!,
Q,
Qrhs,
solver::ConjugateGradient,
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! outputsQrhs
: (array). This is an object that linearoperator! acts onsolver
: (struct). This is a scruct for dispatch, in this case for ColumnwisePreconditionedConjugateGradientargs...
: (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
JacobianFreeNewtonKrylovSolver initialize the residual
ClimateMachine.SystemSolvers.doiteration!
— Functiondoiteration!(
linearoperator!,
preconditioner,
Q,
Qrhs,
solver::ConjugateGradient,
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! overwritesQrhs
: (array). This is an object that linearoperator! acts on. This is the rhs to the linear systemsolver
: (struct). This is a scruct for dispatch, in this case for ConjugateGradientthreshold
: (float). Either an absolute or relative toleranceapplyPC!
: (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 solverthreshold
: (float). The value of the residual for the first timestep
Comment
- This function does conjugate gradient
doiteration!(
linearoperator!,
preconditioner,
Q::MPIStateArray,
Qrhs::MPIStateArray,
solver::ConjugateGradient,
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! overwritesQrhs
: (array). This is an object that linearoperator! acts on. This is the rhs to the linear systemsolver
: (struct). This is a scruct for dispatch, in this case for ConjugateGradientthreshold
: (float). Either an absolute or relative toleranceapplyPC!
: (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 solverthreshold
: (float). The value of the residual for the first timestep
Comment
- This function does conjugate gradient
LU Decomposition
ClimateMachine.SystemSolvers.ManyColumnLU
— TypeManyColumnLU()
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.
ClimateMachine.SystemSolvers.SingleColumnLU
— TypeSingleColumnLU()
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.
Preconditioners
ClimateMachine.SystemSolvers.NoPreconditioner
— Typemutable struct NoPreconditioner end
Do nothing
ClimateMachine.SystemSolvers.ColumnwiseLUPreconditioner
— Typemutable 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 formatQ
: MPIArray container, used to update APQ
: MPIArray container, used to update Acounter
: count the number of Newton, when counter > update_freq or counter < 0, update preconditionupdate_freq
: preconditioner update frequency
...
Shared components
ClimateMachine.SystemSolvers.AbstractSystemSolver
— TypeAbstractSystemSolver
This is an abstract type representing a generic linear solver.
ClimateMachine.SystemSolvers.AbstractNonlinearSolver
— TypeAbstractNonlinearSolver
This is an abstract type representing a generic nonlinear solver.
ClimateMachine.SystemSolvers.AbstractIterativeSystemSolver
— TypeAbstractIterativeSystemSolver
This is an abstract type representing a generic iterative linear solver.
The available concrete implementations are:
ClimateMachine.SystemSolvers.AbstractPreconditioner
— TypeAbstract base type for all preconditioners.
ClimateMachine.SystemSolvers.nonlinearsolve!
— FunctionSolving rhs!(Q) = Qrhs via Newton,
where F = rhs!(Q) - Qrhs
dF/dQ(Q^n) ΔQ ≈ jvp!(ΔQ; Q^n, F(Q^n))
preconditioner ≈ dF/dQ(Q)
ClimateMachine.SystemSolvers.linearsolve!
— Functionlinearsolve!(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.
ClimateMachine.SystemSolvers.settolerance!
— Functionsettolerance!(solver::AbstractIterativeSystemSolver, tolerance, relative)
Sets the relative or absolute tolerance of the iterative linear solver solver
to tolerance
.
ClimateMachine.SystemSolvers.prefactorize
— Functionprefactorize(linop!, linearsolver, args...)
Prefactorize the in-place linear operator linop!
for use with linearsolver
.
ClimateMachine.SystemSolvers.preconditioner_update!
— FunctionDo nothing, when there is no preconditioner, preconditioner = Nothing
update the DGColumnBandedMatrix by the finite difference approximation ...
Arguments
op
: operator used to compute the finte difference informationdg
: the DG model, use only the grid information
...
ClimateMachine.SystemSolvers.preconditioner_solve!
— FunctionDo nothing, when there is no preconditioner, preconditioner = Nothing
Inplace applying the preconditioner
Q = P⁻¹ * Q
ClimateMachine.SystemSolvers.preconditioner_counter_update!
— FunctionDo nothing, when there is no preconditioner, preconditioner = Nothing
Update the preconditioner counter, after each Newton iteration