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:

\mathsf{A} \mathsf{x} = \mathsf{b}

with a \mathsf{A} a square sparse matrix, \mathsf{x} the solution vector, and \mathsf{b} the right-hand side. For example, in a classical linear elastostatics problem \mathsf{A} is the stiffness matrix, and \mathsf{x} and \mathsf{b} 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 GEOS both options are made available wrapping around well-established open-source 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 \mathsf{A}, is small. Irrespective of the selected direct solver implementation, three stages can be idenitified:

  1. Setup Stage: the matrix is first analyzed and then factorized

  2. Solve Stage: the solution to the linear systems involving the factorized matrix is computed

  3. Finalize Stage: the systems involving the factorized matrix have been solved and the direct solver lifetime ends

The default option in GEOS 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 \mathsf{A}, 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 trade-off 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 \mathsf{M} is available, three standard approaches are used to apply the preconditioner:

  1. Left preconditioning: the preconditioned system is \mathsf{M}^{-1} \mathsf{A} \mathsf{x} = \mathsf{M}^{-1} \mathsf{b}

  2. Right preconditioning: the preconditioned system is \mathsf{A} \mathsf{M}^{-1} \mathsf{y} = \mathsf{b}, with \mathsf{x} = \mathsf{M}^{-1} \mathsf{y}

  3. Split preconditioning: the preconditioned system is \mathsf{M}^{-1}_L \mathsf{A} \mathsf{M}^{-1}_R \mathsf{y} = \mathsf{M}^{-1}_L \mathsf{b}, with \mathsf{x} = \mathsf{M}^{-1}_R \mathsf{y}

Summary

The following table summarizes the available input parameters for the linear solver.

Name

Type

Default

Description

amgAggressiveCoarseningLevels

integer

0

AMG number of levels for aggressive coarsening

amgAggressiveCoarseningPaths

integer

1

AMG number of paths for aggressive coarsening

amgAggressiveInterpType

geos_LinearSolverParameters_AMG_AggInterpType

multipass

AMG aggressive interpolation algorithm. Available options are: default\|extendedIStage2\|standardStage2\|extendedStage2\|multipass\|modifiedExtended\|modifiedExtendedI\|modifiedExtendedE\|modifiedMultipass

amgCoarseSolver

geos_LinearSolverParameters_AMG_CoarseType

direct

AMG coarsest level solver/smoother type. Available options are: default\|jacobi\|l1jacobi\|fgs\|sgs\|l1sgs\|chebyshev\|direct\|bgs

amgCoarseningType

geos_LinearSolverParameters_AMG_CoarseningType

HMIS

AMG coarsening algorithm. Available options are: default\|CLJP\|RugeStueben\|Falgout\|PMIS\|HMIS

amgInterpolationMaxNonZeros

integer

4

AMG interpolation maximum number of nonzeros per row

amgInterpolationType

geos_LinearSolverParameters_AMG_InterpType

extendedI

AMG interpolation algorithm. Available options are: default\|modifiedClassical\|direct\|multipass\|extendedI\|standard\|extended\|directBAMG\|modifiedExtended\|modifiedExtendedI\|modifiedExtendedE

amgNullSpaceType

geos_LinearSolverParameters_AMG_NullSpaceType

constantModes

AMG near null space approximation. Available options are:constantModes\|rigidBodyModes

amgNumFunctions

integer

1

AMG number of functions

amgNumSweeps

integer

1

AMG smoother sweeps

amgRelaxWeight

real64

1

AMG relaxation factor for the smoother

amgSeparateComponents

integer

0

AMG apply separate component filter for multi-variable problems

amgSmootherType

geos_LinearSolverParameters_AMG_SmootherType

l1sgs

AMG smoother type. Available options are: default\|jacobi\|l1jacobi\|fgs\|bgs\|sgs\|l1sgs\|chebyshev\|ilu0\|ilut\|ic0\|ict

amgThreshold

real64

0

AMG strength-of-connection threshold

directCheckResidual

integer

0

Whether to check the linear system solution residual

directColPerm

geos_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

geos_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 Eisenstat-Walker adaptive linear tolerance

krylovMaxIter

integer

200

Maximum iterations allowed for an iterative solver

krylovMaxRestart

integer

200

Maximum iterations before restart (GMRES only)

krylovTol

real64

1e-06

Relative convergence tolerance of the iterative method
If the method converges, the iterative solution \mathsf{x}_k is such that
the relative residual norm satisfies:
\left\lVert \mathsf{b} - \mathsf{A} \mathsf{x}_k \right\rVert_2 < krylovTol * \left\lVert\mathsf{b}\right\rVert_2

krylovWeakestTol

real64

0.001

Weakest-allowed tolerance for adaptive method

logLevel

integer

0

Log level

preconditionerType

geos_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

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

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

  1. The number of reduction levels and the coarse points (corresponding to fields) for each level. For example, for single-phase hydraulic fracturing, there are two fields, i.e. displacement and fluid pressure, a two-level MGR strategy can be used with the fluid pressure being the coarse degrees of freedom.

  2. Interpolation/restriction operators and the coarse-grid 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 non-Galerkin approach is preferable.

  3. Global smoother. Depending on the problem, a global relaxation step could be beneficial. Some options include ILU(k), (block) Jacobi, (block) Gauss-Seidel.

  4. Solvers for F-relaxation and coarse-grid 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 F-points. For example, if the \mathsf{A}_{FF} block is hyperbolic, a Jacobi smoother is sufficient while for an elliptic operator an AMG V-cycle might be required. For the single-phase hydraulic fracturing case, an AMG V-cycle is needed for both F-relaxation and coarse-grid 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 F-points, 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 \mathsf{S} = \mathsf{A}_{11} - \mathsf{A}_{10} \mathsf{\widetilde{A}}_{00}^{-1} \mathsf{A}_{01} computation, that requires an approximation of the leading block. Currently, available options for \mathsf{\widetilde{A}}_{00}^{-1} are:

  • diagonal with diagonal values (essentially, a Jacobi preconditioner);

  • diagonal with row sums as values (e.g., used for CPR-like preconditioners).

Once the Schur complement is computed, to properly define the block preconditioner we need:

  • the preconditioner for \mathsf{A}_{00} (any of the above listed single-matrix preconditioner);

  • the preconditioner for \mathsf{S} (any of the above listed single-matrix 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;

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