4.7. Solving the problem

4.7.1. Nonlinear Solution Strategy

4.7.1.1. Equation System Overview

In coupled physics problems one is tasked with solving a set of nonlinear coupled equations and ideally one would solve these equations as a fully-coupled monolithic system. In practice solving the set of equations monolithically may prove impractical owing to large linear system memory requirements and difficulty in locating a small radius of convergence. Thus in some cases one may wish to perform a segregated solve of the equations where each set of equations can be solved in a slightly different manner. This ability to solve subsets of the entire equation set is supported through a feature called Equation System by which criteria for solving a set of nonlinear equations can be applied to a select group of equations.

By default a single set of coupled equations is treated internally as an Equation System without need of additional input syntax. Here the solution strategy for the Equation System is solely determined by a collection of nonlinear solution specifications. A typical setup for a single equation system might be

Begin Aria Region myRegion
  .
  Use Linear Solver solve_temperature
  .
  eq energy for temperature On block_1 using q1 with lumped_mass diff
  .
  nonlinear solution strategy = newton
  use dof averaged nonlinear residual
  maximum nonlinear iterations = 20
  nonlinear residual tolerance = 1.e-6
  .
  .
End Aria Region myRegion

Simulations in which different subsets of the coupled equations are solved in a segregated manner require that each subset be explicitly defined as an Equation System with a corresponding solution strategy which may differ for each equation system. When multiple Equation Systems are present, the overall segregated solution is governed by segregated iteration and convergence criteria at an outer level.

A typical setup for segregated solution of multiple equation systems might be

Begin Aria Region myRegion
  .
  Segregated global nonlinear residual tolerance = 6.0e-9
  .
  Begin Equation System temperature
    .
    eq energy for temperature On block_1 using q1 with lumped_mass diff
    Use Linear Solver solve_temperature
    .
    nonlinear solution strategy = newton
    use dof averaged nonlinear residual
    maximum nonlinear iterations = 20
    nonlinear residual tolerance = 1.e-6
    .
  End Equation System temperature
  .
  Begin Equation System species
    .
    eq species for species_0 On block_1 using q1 with diff
    .
    Use Linear Solver solve_species
    .
    nonlinear solution strategy = newton
    use dof averaged nonlinear residual
    maximum nonlinear iterations = 10
    nonlinear residual tolerance = 1.e-8
    .
  End Equation System species
  .
End Aria Region myRegion

4.7.1.2. Nonlinear Solution Specifications

The general solution strategy used in Aria is to solve the requested system of equations for solution variables $\symSlnVec$ incrementally. In the context of a Newton formulation then within the Newton step one will solve

$\symJac \Delta \symSlnVec = -R(\symSlnVec)$

where $\symJac$ is the Jacobian matrix, $\Delta \symSlnVec$ is the solution increment and $R(\symSlnVec)$ is the system residual for the governing equations. See the theory section for more information on the solution strategy. This chapter describes basic nonlinear solution command line settings that influence how a converged solution will be obtained. For transient solutions these settings are used in combination with time integration settings for advancing a solution to the next time step.

It is instructive here to review a typical log file output of the Aria solution summary in order to establish how a user might influence the solution procedure.

Equation System AriaRegion->main:
* Step  : Transient, Strategy: NEWTON, Time: 2.20e+03, Step: 5.29e+00
* Matrix: Solver: "Solve_Temperature", Unknowns: 1360, Nonzeros: 11578
* Mesh  : Processor 0 of 1: 1253 of 1253 elems, 1363 of 1363 nodes
* Computing View Factors for enclosure space
  N O N L I N E A R                L I N E A R
----------------------   ------------------------------------
Step  Resid    Delta     Itns Status  Resid    Asm/Slv Time
---- -------- --------   ---- ------ -------- ---------------
1    6.37e-02 1.10e+01   59    ok    8.23e-07 7.0e-03/2.0e-03
2    1.10e-04 8.18e-03   58    ok    1.19e-09 7.0e-03/2.0e-03
3    3.06e-07 8.81e-05   61    ok    3.68e-12 6.0e-03/2.0e-03
Termination reason: 8.80515e-05 < nonlinear_correction_tolerance(0.001)

In the log file output shown above, columns under the N O N L I N E A R, heading are related to the command lines described in this chapter. The columns Resid and Delta represent L2 norms of an aggregate residual $R^{*}(\symSlnVec)$ and solution increment vector $\Delta\symSlnVec$ while the column denoted Step denotes the iteration within this Newton step. Columns under the L I N E A R heading represent the number of linear solver iterations (Itns), the Resid column represents the linear system residual $R^{*}(\Delta \symSlnVec)$ . Here $R^{*}(\Delta\symSlnVec) = \symJac \Delta \symSlnVec + R(\symSlnVec)$ , the residual for the assembled system of equations, is solved for $\Delta\symSlnVec$ to the tolerance requested in the Linear Solver specifications. The column labeled Status indicates whether the linear solve completed successfully, ok or fail, within the allowable number of iterations defined in the Linear Solver specifications. The column labeled Asm/Slv Time contain assembly and solution times for the system of equations.

Given a current value of $\symSlnVec$ then within each Newton step the residual $R^{*}(\Delta\symSlnVec)$ is first evaluated using $\Delta\symSlnVec=0$ and tested for convergence. Here the convergence conditions are satisfied when either the NONLINEAR RESIDUAL TOLERANCE or NONLINEAR CORRECTION TOLERANCE are achieved. If the solution is not yet converged, the system of equations is turned over to the equation solver where it is solved for a new $\Delta\symSlnVec$ subject to the Linear Solver specifications. Then the Solver solution is used to update the nonlinear solution $\symSlnVec=\symSlnVec_{old}+\Delta\symSlnVec$ and the process begins anew subject to a user supplied MAXIMUM NONLINEAR ITERATIONS.

Within the Newton step a solution is said to FAIL if it is not converged within MAXIMUM NONLINEAR ITERATIONS and this occurrence will be noted in the log file. For transient simulations one will normally repeat the solution procedure with a smaller time step as described in Time Step Selection. In some cases the user may chose to override the FAIL condition using the ACCEPT SOLUTION AFTER MAXIMUM NONLINEAR ITERATIONS command line but one is cautioned that this option should be used with discretion as it may discredit any expectation of accuracy in transient results. For completeness, a summary of pass, fail, and fail but accepted steps are included at the end of the log file.

For most applications the default NONLINEAR RESIDUAL TOLERANCE (1e-6) is adequate when accompanied by a Linear Solver specification of RESIDUAL NORM TOLERANCE minimum of 1.0e-8. Inclusion of the command line USE DOF AVERAGED NONLINEAR RESIDUAL is highly recommended, particularly when performing uniform mesh refinement studies. A summary of these and other nonlinear solution related command lines is given in what follows.

4.7.2. Linear Solution Strategy

4.7.2.1. Overview

Aria supports several choices of linear solver packages. Starting with the 5.0 release the Tpetra based linear solvers are the only linear solvers supported.

Within the input file linear solver command block will appear at the Domain scope as illustrated below. Note that different linear solvers can be used for different equation systems. Furthermore, a solver command block can appear in the input but not be used.

Begin Sierra
  .
  .
  Begin Tpetra Equation Solver third_solver
    .
    Tpetra linear solver commands
    .
  End
  .
  Begin Procedure My_Aria_Procedure
     .
     Begin Aria Region My_Region
       .
       use linear solver first_solver
       .
     End

     Begin Aria Region My_Other_Region
       .
       use linear solver third_solver
       .
     End

  End
  .
End Sierra

See Solver Selection Guidelines for suggestions on how to select a type of linear solver for different kinds of problems Aria can solve.

4.7.3. 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 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 Preconditioners for details on the supported preconditioners and their available options.

4.7.4. 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.7.5. 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 Aria 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.7.6. Solver Selection Guidelines

Selecting an appropriate linear solver that is both robust and efficient for a particular Aria 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:

As a starting point preset solvers are provided for several common types of problems. These are recommended as a starting point for any analysis.

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

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, an equivalent solver block is written to the log file and can be used as a starting point for modifying additional solver or preconditioner settings. This solver block will contain the preset preconditioner along with default solver settings. In reusing the this block, users will often choose to increase the number of MAXIMUM ITERATIONS in hopes of lowering the solution time, but it is noted that increases beyond triple the default setting reap low reward in obtaining a solution.

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.
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 Preconditioners for setting the fill and overlap.
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 Aria, 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.

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. Aria 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. The default MUELU preconditioner is not recommended for general use. However, it is sometimes useful for steady state temperature simulations, Poisson problems that are nearly singular and difficult to solve with preset solvers.

4.7.7. 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 1e-6 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 Aria, 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 Aria 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

4.7.8. Block Preconditioning

Beta Capability

This features is currently under active development and research, and is not recommended for production use.

Aria supports block preconditioning through the Teko package. Currently the only way to use Teko is through the expert interface that uses an XML file specify the relevant Teko parameters. Please see https://trilinos.github.io/teko.html for more information on how to structure this XML. An example input file block is shown below.

BEGIN TPETRA EQUATION SOLVER  GMRES_DDILU
  BEGIN GMRES SOLVER
    BEGIN EXPERT TEKO PRECONDITIONER
      INVERSE METHOD = GS-Prec
      DEFINE MATRIX SUBBLOCK number=1 dof=Velocity
      DEFINE MATRIX SUBBLOCK number=2 dof=Pressure
      XML FILE = teko.xml
    END
    MAXIMUM ITERATIONS = 1000
    RESIDUAL SCALING = NONE
    CONVERGENCE TOLERANCE = AUTOMATIC MINIMUM = 1.0E-8
  END
END TPETRA EQUATION SOLVER

where INVERSE METHOD corresponds to the name of the Teko parameterlist found in the XML, and the DEFINE MATRIX SUBBLOCK commands define the dof-to-block mapping needed by the block preconditioner. An example of the Teko XML is shown as

<ParameterList>
<ParameterList name="GS-Prec">
  <Parameter name="Type" type="string" value="Block Gauss-Seidel"/>
  <Parameter name="Inverse Type 1" type="string" value="Diagonal 1"/>
  <Parameter name="Inverse Type 2" type="string" value="Diagonal 1"/>
</ParameterList>

<ParameterList name="Diagonal 1">
  <Parameter name="Type" type="string" value="Ifpack2"/>
  <ParameterList name="Ifpack2 Settings">
    <Parameter name="relaxation: type" type="string" value="Jacobi"/>
    <Parameter name="relaxation: sweeps" type="int" value="1"/>
  </ParameterList>
</ParameterList>
</ParameterList>

In the above block, “GS-Prec” is the main preconditioning method, with the two subblock inverses described by “Inverse Type 1” and “Inverse Type 2” corresponding to a Jacobi preconditioner in the “Diagonal 1” parameter list. Block preconditioning can also be used to accelerate incompressible Navier–Stokes flow through schemes inspired by the SIMPLE class of operating splitting techniques. An example Teko XML file for such a use case is shown as

<ParameterList name="SIMPLE">
  <Parameter name="Type" type="string" value="NS SIMPLE"/>
  <Parameter name="Inverse Velocity Type" type="string" value="ddilu"/>
  <Parameter name="Inverse Pressure Type" type="string" value="Amesos2"/>
  <Parameter name="Explicit Velocity Inverse Type" type="string" value="AbsRowSum"/>
</ParameterList>

<ParameterList name="ddilu">
 <Parameter name="Type" type="string" value="Ifpack2"/>
 <Parameter name="Prec Type" type="string" value="RILUK"/>
  <ParameterList name="Ifpack2 Settings">
    <Parameter name="fact: iluk level-of-fill" type="int" value="1"/>
    <Parameter name="schwarz: use reordering" type="bool" value="true"/>
  </ParameterList>
</ParameterList>

In the above case, a SIMPLEC method is used with the velocity inverse approximated using an RILUK preconditioner and the pressure (Schur complement approximation) inverse is solved directly using an Amesos2 direct solver. Please note that in practice, using MueLu to approximate the pressure inverse would demonstrate better scaling, but at the current state of development unknown problems exist using MueLu inside Teko leading to a lack of convergence.

To get additional diagnostics about the residual of the individual inverse approximations, Teko allows for the injection of a diagnostic inverse type. To do so, one can modify the above XML to inject a diagnostic block inside the pressure inverse as shown

<ParameterList name="MyPress">
  <Parameter name="Type" type="string" value="Diagnostic Inverse"/>
  <Parameter name="Inverse Factory" type="string" value="Amesos2"/>
  <Parameter name="Descriptive Label" type="string" value="Inverse Pressure"/>
  <Parameter name="Print Residual" type="bool" value="true"/>
</ParameterList>

<ParameterList name="SIMPLE">
  <Parameter name="Type" type="string" value="NS SIMPLE"/>
  <Parameter name="Inverse Velocity Type" type="string" value="ddilu"/>
  <Parameter name="Inverse Pressure Type" type="string" value="MyPress"/>
  <Parameter name="Explicit Velocity Inverse Type" type="string" value="AbsRowSum"/>
</ParameterList>

<ParameterList name="ddilu">
 <Parameter name="Type" type="string" value="Ifpack2"/>
 <Parameter name="Prec Type" type="string" value="RILUK"/>
  <ParameterList name="Ifpack2 Settings">
    <Parameter name="fact: iluk level-of-fill" type="int" value="1"/>
    <Parameter name="schwarz: use reordering" type="bool" value="true"/>
  </ParameterList>
</ParameterList>

where the pressure inverse residual will be printed to stdout with the “Inverse Pressure” label. Such a manner of injecting diagnostic inverses is helpful for determining which portion of a block preconditioning scheme is failing to converge.