Linear Solvers¶
Introduction¶
Any physics solver relying on standard finite element and finite volume techniques requires the solution of algebraic linear systems, which are obtained upon linearization and discretization of the governing equations, of the form:
with a a square sparse matrix, the solution vector, and the righthand side. For example, in a classical linear elastostatics problem is the stiffness matrix, and and are the displacement and nodal force vectors, respectively.
This solution stage represents the most computationally expensive portion of a typical simulation. Solution algorithms generally belong to two families of methods: direct methods and iterative methods. In GEOSX both options are made available wrapping around wellestablished opensource linear algebra libraries, namely HYPRE, PETSc, SuperLU, and Trilinos.
Direct methods¶
The major advantages are their reliability, robustness, and ease of use. However, they have large memory requirements and exhibit poor scalability. Direct methods should be used in a prototyping stage, for example when developing a new formulation or algorithm, when the dimension of the problem, namely the size of matrix , is small. Irrespective of the selected direct solver implementation, three stages can be idenitified:
 Setup Stage: the matrix is first analyzed and then factorized
 Solve Stage: the solution to the linear systems involving the factorized matrix is computed
 Finalize Stage: the systems involving the factorized matrix have been solved and the direct solver lifetime ends
The default option in GEOSX relies on SuperLU, a general purpose library for the direct solution of large, sparse, nonsymmetric systems of linear equations, that is called taking advantage of the interface provided in HYPRE.
Iterative methods¶
As the problem size (number of computational cells) increases, global iterative solution strategies are the method of choice—typically nonsymmetric Krylov solvers. Because of the possible poor conditioning of , preconditioning is essential to solve such systems efficiently. ‘’Preconditioning is simply a means of transforming the original linear system into one which has the same solution, but which is likely to be easier to solve with an iterative solver’’ [Saad (2003)].
The design of a robust and efficient preconditioner is based on a tradeoff between two competing objectives:
 Robustness: reducing the number of iterations needed by the preconditioned solver to achieve convergence;
 Efficiency: limiting the time required to construct and apply the preconditioner.
Assuming a preconditioning matrix is available, three standard approaches are used to apply the preconditioner:
 Left preconditioning: the preconditioned system is
 Right preconditioning: the preconditioned system is , with
 Split preconditioning: the preconditioned system is , with
Summary¶
The following table summarizes the available input parameters for the linear solver.
Name  Type  Default  Description 

amgAggresiveCoarseningLevels  integer  0  AMG number levels for aggressive coarsening
Available options are: TODO

amgCoarseSolver  geosx_LinearSolverParameters_AMG_CoarseType  direct  AMG coarsest level solver/smoother type. Available options are: default\jacobi\l1jacobi\fgs\sgs\l1sgs\chebyshev\direct\bgs 
amgCoarseningType  string  HMIS  AMG coarsening algorithm
Available options are: TODO

amgInterpolationType  integer  6  AMG interpolation algorithm
Available options are: TODO

amgNullSpaceType  geosx_LinearSolverParameters_AMG_NullSpaceType  constantModes  AMG near null space approximation. Available options are:constantModes\rigidBodyModes 
amgNumFunctions  integer  1  AMG number of functions
Available options are: TODO

amgNumSweeps  integer  2  AMG smoother sweeps 
amgSmootherType  geosx_LinearSolverParameters_AMG_SmootherType  fgs  AMG smoother type. Available options are: default\jacobi\l1jacobi\fgs\bgs\sgs\l1sgs\chebyshev\ilu0\ilut\ic0\ict 
amgThreshold  real64  0  AMG strengthofconnection threshold 
directCheckResidual  integer  0  Whether to check the linear system solution residual 
directColPerm  geosx_LinearSolverParameters_Direct_ColPerm  metis  How to permute the columns. Available options are: none\MMD_AtplusA\MMD_AtA\colAMD\metis\parmetis 
directEquil  integer  1  Whether to scale the rows and columns of the matrix 
directIterRef  integer  1  Whether to perform iterative refinement 
directParallel  integer  1  Whether to use a parallel solver (instead of a serial one) 
directReplTinyPivot  integer  1  Whether to replace tiny pivots by sqrt(epsilon)*norm(A) 
directRowPerm  geosx_LinearSolverParameters_Direct_RowPerm  mc64  How to permute the rows. Available options are: none\mc64 
iluFill  integer  0  ILU(K) fill factor 
iluThreshold  real64  0  ILU(T) threshold factor 
krylovAdaptiveTol  integer  0  Use EisenstatWalker adaptive linear tolerance 
krylovMaxIter  integer  200  Maximum iterations allowed for an iterative solver 
krylovMaxRestart  integer  200  Maximum iterations before restart (GMRES only) 
krylovTol  real64  1e06  Relative convergence tolerance of the iterative method
If the method converges, the iterative solution is such that
the relative residual norm satisfies:
<
krylovTol * 
krylovWeakestTol  real64  0.001  Weakestallowed tolerance for adaptive method 
logLevel  integer  0  Log level 
preconditionerType  geosx_LinearSolverParameters_PreconditionerType  iluk  Preconditioner type. Available options are: none\jacobi\l1jacobi\fgs\sgs\l1sgs\chebyshev\iluk\ilut\icc\ict\amg\mgr\block\direct\bgs 
solverType  geosx_LinearSolverParameters_SolverType  direct  Linear solver type. Available options are: direct\cg\gmres\fgmres\bicgstab\preconditioner 
stopIfError  integer  1  Whether to stop the simulation if the linear solver reports an error 
Preconditioner descriptions¶
This section provides a brief description of the available preconditioners.
 None: no preconditioning is used, i.e., .
 Jacobi: diagonal scaling preconditioning, with , with the matrix diagonal. Further details can be found in:
 ILUK: incomplete LU factorization with fill level k of the original matrix: . Further details can be found in:
 ILUT: a dual threshold incomplete LU factorization: .
Further details can be found in:
 HYPRE documentation,
 not yet available through PETSc interface,
 Trilinos documentation.
 ICC: incomplete Cholesky factorization of a symmetric positive definite matrix: .
Further details can be found in:
 not yet available through hypre interface,
 PETSc documentation,
 Trilinos documentation.
 AMG: algebraic multigrid (can be classical or aggregationbased according to the specific package). Further details can be found in:
 MGR: multigrid reduction. Available through hypre interface only. Specific documentation coming soon. Further details can be found in MGR documentation.
 Block: custom preconditioner designed for a 2 x 2 block matrix.
HYPRE MGR Preconditioner¶
MGR stands for multigrid reduction, a multigrid method that uses the interpolation, restriction operators, and the Galerkin triple product, to reduce a linear system to a smaller one, similar to a Schur complement approach. As such, it is designed to target block linear systems resulting from discretizations of multiphysics problems. GEOSX uses MGR through an implementation in HYPRE. More information regarding MGR can be found here. Currently, MGR strategies are implemented for hydraulic fracturing, poroelastic, compositional flow with and without wells. More multiphysics solvers with MGR will be enabled in the future.
To use MGR for a specific block system, several components need to be specified.
 The number of reduction levels and the coarse points (corresponding to fields) for each level. For example, for singlephase hydraulic fracturing, there are two fields, i.e. displacement and fluid pressure, a twolevel MGR strategy can be used with the fluid pressure being the coarse degrees of freedom.
 Interpolation/restriction operators and the coarsegrid computation strategy. A simple but effective strategy is to use Jacobi diagonalization for interpolation and injection for restriction. For most cases, a Galerkin coarse grid strategy can be used, but for special cases such as poroelastic, a nonGalerkin approach is preferable.
 Global smoother. Depending on the problem, a global relaxation step could be beneficial. Some options include ILU(k), (block) Jacobi, (block) GaussSeidel.
 Solvers for Frelaxation and coarsegrid correction. These solvers should be chosen carefully for MGR to be effective. The choice of these solvers should correspond to the properties of the blocks specified by the C and Fpoints. For example, if the block is hyperbolic, a Jacobi smoother is sufficient while for an elliptic operator an AMG Vcycle might be required. For the singlephase hydraulic fracturing case, an AMG Vcycle is needed for both Frelaxation and coarsegrid correction.
Note that these are only general guidelines for designing a working MGR recipe. For complicated multiphysics problems, experimentation with different numbers of levels, choices of C and Fpoints, and smoothers/solvers, etc., is typically needed to find the best strategy. Currently, these options are only available to developers. We are working on exposing these functionalities to the users in future releases.
Block preconditioner¶
This framework allows the user to design a block preconditioner for a 2 x 2 block matrix. The key component is the Schur complement computation, that requires an approximation of the leading block. Currently, available options for are:
 diagonal with diagonal values (essentially, a Jacobi preconditioner);
 diagonal with row sums as values (e.g., used for CPRlike preconditioners).
Once the Schur complement is computed, to properly define the block preconditioner we need:
 the preconditioner for (any of the above listed singlematrix preconditioner);
 the preconditioner for (any of the above listed singlematrix preconditioner);
 the application strategy. This can be:
 diagonal: none of the coupling terms is used;
 upper triangular: only the upper triangular coupling term is used;
 lowerupper triangular: both coupling terms are used.
Moreover, a block scaling is available. Feasible options are:
 none: keep the original scaling;
 Frobenius norm: equilibrate Frobenius norm of the diagonal blocks;
 user provided.