4.3. Linear Solvers

Starting with the 5.0 release the Tpetra based linear solvers are the only linear solvers supported.

Within the input file the linear solver command block may optionally appear at the Domain scope as illustrated below. Note that different linear solvers can be used for different equations. A solver command block can appear in the input but not be used, and if no solver command blocks are explicitly specified Fuego will use sensible default linear solvers for each equation.

Begin Sierra
  .
  .
  Begin Tpetra Equation Solver first_solver
    .
    Tpetra linear solver commands
    .
  End
  Begin Tpetra Equation Solver second_solver
    .
    Tpetra linear solver commands
    .
  End
  .
  Begin Procedure My_Fuego_Procedure
     .
     Begin Fuego Region My_Region
       .
       use equation solver first_solver for equation continuity
       use equation solver second_solver for equation X_Momentum
       use equation solver second_solver for equation Y_Momentum
       use equation solver second_solver for equation Z_Momentum
       .
     End

  End
  .
End Sierra

See Solver Selection Guidelines for suggestions on how to select a type of linear solver.

4.3.1. Tpetra Solver Overview

The syntax for specifying a Tpetra-based solver uses a nested command block structure. The outermost level at the Sierra scope is the BEGIN TPETRA EQUATION SOLVER block. This block must contain a single nested solver block corresponding to your desired solver type. Several types of direct and iterative linear solvers are provided through the Amesos2 and Belos packages of Trilinos. The supported solver blocks are listed in Tpetra Solvers.

If one of the iterative solvers is selected, it must also have a nested preconditioner block. A variety of preconditioners are supported through the Ifpack2 and MueLu Trilinos packages. See Tpetra Preconditioners for details on the supported preconditioners and their available options.

4.3.2. Solver Error Codes

If the linear solve fails it will report an error code. While these codes are not universally standardized across linear solvers, they can generally be interpreted with the following rules:

Error 1 (or -1)

The linear solver reached the maximum number of linear solver iterations and did not reach a solution. Check the final linear solver residual to see if this is serious or not. Even if the error is not serious (residuals are small) this usually means the solve is much more expensive than necessary and a better solver/preconditioner combination should be used.

Error 2 (or -2, 4, or -4)

The linear solver encountered a fatal internal error. This is often because of a singular matrix or a NaN or Infinite term. This is almost always a serious error.

Error 3 (or -3)

The linear solver had an internal loss of precision. This is often not a serious error, but may indicate the need for a better solver/preconditioner combination. As always, check the final linear solver residual reported in the log file to determine if it is serious or not.

4.3.3. Expert Syntax

The standard solver and preconditioner blocks only support a subset of the commonly used options that are supported by the underlying Trilinos packages. The subset was chosen to cover the vast majority of linear solver use cases within Fuego simulations, but an alternative approach is provided for expert users who would like to be able to pass parameters directly to the underlying Belos, Ifpack2, and MueLu packages. This expert interface is used by providing the path to an XML file containing the definition for a Teuchos::ParameterList that is supported by the underlying Trilinos package in one of the BEGIN EXPERT BELOS SOLVER, BEGIN EXPERT IFPACK2 PRECONDITIONER, or BEGIN MUELU PRECONDITIONER blocks as desired.

Sample XML for expert Belos solver

<ParameterList>
  <Parameter name="solution method" type="string" value="CG"/>
  <ParameterList name="Belos-Tpetra">
    <Parameter name="Convergence Tolerance" type="double" value="1.e-12"/>
    <Parameter name="Maximum Iterations" type="int" value="1000"/>
  </ParameterList>
</ParameterList>

The solution method parameter is required. Any additional parameters that are not specified in the Belos-Tpetra sublist will use their default values as defined by Belos.

Sample XML for expert Ifpack2 preconditioner

<ParameterList>
  <Parameter name="preconditioner type:" type="string" value="ILUT"/>
  <Parameter name="fact: ilut level-of-fill" type="double" value="1.0"/>
  <Parameter name="fact: drop tolerance" type="double" value="0.0"/>
</ParameterList>

The preconditioner type: parameter is required. Any additional parameters that are not specified will use their default values as defined by Ifpack2.

Sample XML for expert MueLu preconditioner

<ParameterList name="MueLu-Tpetra">
  <Parameter        name="verbosity"                        type="string"   value="high"/>
  <Parameter        name="coarse: type"                     type="string"   value="ILUT"/>
  <Parameter        name="coarse: max size"                 type="int"      value="1000"/>
  <Parameter        name="max levels"                       type="int"      value="4"/>
  <Parameter        name="smoother: type"                   type="string"   value="CHEBYSHEV"/>
  <ParameterList    name="smoother: params">
     <Parameter name="chebyshev: degree"                    type="int"      value="2"/>
     <Parameter name="chebyshev: ratio eigenvalue"          type="double"   value="20"/>
     <Parameter name="chebyshev: min eigenvalue"            type="double"   value="1.0"/>
     <Parameter name="chebyshev: zero starting solution"    type="bool"     value="true"/>
     <Parameter name="chebyshev: eigenvalue max iterations" type="int"      value="15"/>
  </ParameterList>
  <Parameter        name="aggregation: type"                type="string"   value="uncoupled"/>
  <Parameter        name="aggregation: drop tol"            type="double"   value="0.02"/>

  <Parameter        name="repartition: enable"              type="bool"     value="true"/>
  <Parameter        name="repartition: min rows per proc"   type="int"      value="1000"/>
  <Parameter        name="repartition: start level"         type="int"      value="2"/>
  <Parameter        name="repartition: max imbalance"       type="double"   value="1.327"/>
  <Parameter        name="repartition: partitioner"         type="string"   value="zoltan2"/>
</ParameterList>

Unspecified parameters will use their default values as defined by MueLu.

For further documentation of the parameters supported by each Trilinos package please see the corresponding documentation on trilinos.org.

4.3.4. Solver Selection Guidelines

Selecting an appropriate linear solver that is both robust and efficient for a particular Fuego simulation can be confusing. In general there is a tradeoff between solver speed and robustness, and this can significantly affect the total time to solution. The following guidelines for choosing a linear solver may be useful:

  1. As a starting point preset solvers are provided for several common types of problems. These are recommended as a starting point for any analysis, and Fuego will use an appropriate preset for each equation type if a specific solver is not specified.

  • For fluids problems the continuity preset is recommended for the continuity equation system, and the scalar_transport preset is recommended for the momentum equation system, as well as any additional scalar advection-diffusion equations. The continuity preset is a GMRES solver with an algebraic multigrid preconditioner provided through the MueLu Trilinos package. The scalar_transport preset is a GMRES solver with a SGS preconditioner.

  • The thermal preset solver type is recommended as the starting point for general thermal analysis problems. It uses a GMRES solver with Jacobi preconditioning and is a good balance between runtime per linear solver iteration and robustness for most thermal problems.

  • The thermal_symmetric and thermal_bicgstab presets can be slightly faster for thermal analysis problems that do not involve contact. They are CG and BiCGSTAB solvers respectively, each with Jacobi preconditioning. These trade some robustness for shorter runtime per linear iteration when compared with the GMRES solver used by the thermal preset.

  • For more general multiphysics problems involving multiple coupled equations in a single equation system the multiphysics preset solver is the recommended starting point. It is a GMRES solver with a DD-ILU preconditioner.

  • When using any of the preset solver types they will print an equivalent solver block to the log file that can be used as a starting point for modifying additional solver or preconditioner settings.

  1. If the suggested preset solver shows evidence of struggling to solve a particular problem in the log file (for example frequently approaching or reaching the maximum number of linear iterations) consider using a stronger solver and/or preconditioner. For the solvers move from CG to BiCGSTAB to GMRES. For the preconditioners move from Jacobi to SGS to DD-ILU or DD-ILUT. Initially leave all preconditioner-specific options to their default values.

  2. If GMRES/DD-ILU still shows evidence of struggling to solve the problem consider increasing the ILU fill and overlap levels. First trying overlap of 1, then fill of 1, then both is a good balance between the cost and robustness of the preconditioner. See Dd-Ilu Preconditioner. for setting the fill and overlap.

  3. If GMRES/DD-ILU with fill and overlap is still insufficient consider:

  • Trying one of the direct solvers (KLU2, SuperLU, or UMFPACK) if your problem is small. These are the most robust option available to Fuego, but significantly more expensive in both memory usage and runtime than the iterative solver options. They also do not scale in parallel effectively, and should not be used for problems requiring more than a single compute node on HPC systems.

  • For multiphysics problems involving many equations in a single monolithic equation system consider moving to a segregated solution scheme involving multiple smaller equation systems. This will make the linear solves easier, but can make nonlinear convergence of the overall problem more challenging. Loosely coupled equations are more likely to be successfully split into separate equation systems without causing nonlinear convergence problems.

  1. Multigrid preconditioning can be extremely effective and exhibit excellent parallel scalability, however setup of multigrid preconditioners is an extremely complex topic that is both problem type and problem size dependent. Fuego supports general algebraic multigrid preconditioning through the MueLu package from Trilinos, and can read MueLu XML files for defining the multigrid setup. Please consult the MueLu documentation for suggestions on selecting a multigrid preconditioner.

4.3.5. Tolerance Settings

In addition to selecting the solver and preconditioner type, selecting appropriate convergence tolerance and scaling settings is also important. The default settings (R0 scaling and AUTOMATIC tolerance) are a reasonable starting place, but may need to be adjusted depending on the problem and nonlinear solver settings. Typically setting the RESIDUAL SCALING to either NONE or R0 is most useful.

  • NONE applies no scaling to the residual for determining linear solver convergence. If the problem has a low residual to being with the solver will only take a few iterations, however if the initial residual is large this option may take more linear solver iterations than is ideal. The linear solver tolerance selected with this option should be at least an order of magnitude lower than the nonlinear residual tolerance used in the Newton solve or the nonlinear problem may have trouble converging.

  • R0 scales the residual by the initial residual. This is useful if you want a looser acceptance tolerance when the problem starts with high residual, however if the problem starts with a small residual this may cause the solver to take an unreasonably large number of iterations. If the initial residual is large and the linear solver tolerance is too loose this can result in the linear solve producing bad/poorly solved results which may cause convergence problems.

An autonomous linear convergence tolerance can be computed within Fuego, and is declared within the Tpetra solvers by setting the parameter CONVERGENCE TOLERANCE = AUTOMATIC in the Tpetra solver block. This option computes a linear convergence tolerance for each nonlinear iteration based on the initial linear residual and the user-specified nonlinear residual tolerance. The user can also specify an optional minimum linear residual tolerance (with a default value of 1\times10^{-6}) to ensure that the computed linear convergence tolerance within Fuego is not used if it’s too small. Below is an example of using this autonomous setting with a specified minimum

BEGIN TPETRA EQUATION SOLVER  GMRES_DDILU
  BEGIN GMRES SOLVER
    BEGIN DD-ILU PRECONDITIONER
    END
    MAXIMUM ITERATIONS = 1000
    RESIDUAL SCALING = NONE
    CONVERGENCE TOLERANCE = AUTOMATIC MINIMUM = 1.0E-8
  END
END TPETRA EQUATION SOLVER