# System Solvers

## Non-linear solvers

`ClimateMachine.SystemSolvers.LSOnly`

— Type`LSOnly`

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 Q`Qdq::AT`

: container for Q + ϵΔQ`Fq::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`

— Type`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**

### Generalized Minimal Residual Method

`ClimateMachine.SystemSolvers.GeneralizedMinimalResidual`

— Type`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**

### Batched Generalized Minimal Residual Method

`ClimateMachine.SystemSolvers.BatchedGeneralizedMinimalResidual`

— Type```
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.

### Conjugate Gradient Solver Method

`ClimateMachine.SystemSolvers.ConjugateGradient`

— Type```
ConjugateGradient(
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 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**

- 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 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**

- ConjugateGradient struct

`ClimateMachine.SystemSolvers.initialize!`

— Function```
initialize!(
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! 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

JacobianFreeNewtonKrylovSolver initialize the residual

`ClimateMachine.SystemSolvers.doiteration!`

— Function```
doiteration!(
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! 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

```
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! 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

### LU Decomposition

`ClimateMachine.SystemSolvers.ManyColumnLU`

— Type`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.

`ClimateMachine.SystemSolvers.SingleColumnLU`

— Type`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.

## 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 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

...

## Shared components

`ClimateMachine.SystemSolvers.AbstractSystemSolver`

— Type`AbstractSystemSolver`

This is an abstract type representing a generic linear solver.

`ClimateMachine.SystemSolvers.AbstractNonlinearSolver`

— Type`AbstractNonlinearSolver`

This is an abstract type representing a generic nonlinear solver.

`ClimateMachine.SystemSolvers.AbstractIterativeSystemSolver`

— Type`AbstractIterativeSystemSolver`

This is an abstract type representing a generic iterative linear solver.

The available concrete implementations are:

`ClimateMachine.SystemSolvers.AbstractPreconditioner`

— Type`Abstract 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!`

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

`ClimateMachine.SystemSolvers.settolerance!`

— Function`settolerance!(solver::AbstractIterativeSystemSolver, tolerance, relative)`

Sets the relative or absolute tolerance of the iterative linear solver `solver`

to `tolerance`

.

`ClimateMachine.SystemSolvers.prefactorize`

— Function`prefactorize(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 information`dg`

: 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