EquationSolver

Mathematical description

EquationSolver is an object used for solving nonlinear systems of equations of the form

\[F(\mathbf{X}) = \mathbf{0}\]

where \(\mathbf{X}\) is a parallel distributed vector, e.g. mfem::HypreParVector, and \(F\) is a square nonlinear operator, i.e. the dimension of the input vector equals the dimension of the output vector. These systems commonly arise from finite element discretizations such as the equations of heat transfer

\[\mathbf{Mu}_{n+1} + \Delta t (\mathbf{Ku}_{n+1} + f(\mathbf{u}_{n+1})) - \Delta t \mathbf{G} - \mathbf{Mu}_n = \mathbf{0}\]

where \(\mathbf{X} = \mathbf{u}_{n+1}\). These systems are commonly solved by Newton-type methods which have embedded linear solvers for a linearized approximation of the full nonlinear operator. A simple Newton algorithm for solving this system is given below:

Pick an initial guess X = X_0
Compute the residual r = F(X)
while (r is too large) {
  Compute the linearization (Jacobian) of F, J = dF/dX
  Solve J X_update = -r using a linear solver
  X = X + X_update
  Compute the updated residual r = F(X)
}

To perform this solve, we typically need to configure both the nonlinear solver algorithm and the linear solver algorithm. If we want to use an iterative solver for the linear part, we also may need to configure a preconditioner for the system.

Class design

EquationSolver provides an interface to the associated nonlinear and linear solver algorithms needed to solve these systems of equations. Note that while some nonlinear solvers do not depend on an embedded linear solver (e.g. L-BFGS), we require a linear solver to be specified as it is used to compute reasonable initial guesses and perform adjoint solves within Serac.

The key methods provided by this class are:

  1. void SetOperator(const mfem::Operator& op): This defines the nonlinear mfem::Operator representing the vector-valued nonlinear system of equations \(F\).

  2. void Solve(mfem::Vector& x): This solves the nonlinear system \(F(\mathbf{X}) = \mathbf{0}\) and stores the solution in-place in x.

  3. mfem::Solver& LinearSolver(): This returns the associated linear solver for adjoint solves and initial guess calculations.

Two workflows exist for defining the linear and nonlinear algorithms in a EquationSolver: predefined option structs for common use cases and fully custom mfem::Solvers for advanced users.

Common configurations via option structs

The easist way to build an equation solver is by providing structs containing parameters for common configurations of the linear and nonlinear solvers. They can also be passed directly to physics module constructors.

  /**
   * @brief Construct an equation solver object using nonlinear and linear solver option structs.
   *
   * @param nonlinear_opts The options to configure the nonlinear solution scheme
   * @param lin_opts The options to configure the underlying linear solution scheme to be used by the nonlinear solver
   * @param comm The MPI communicator for the supplied nonlinear operators and HypreParVectors
   */
  EquationSolver(NonlinearSolverOptions nonlinear_opts = {}, LinearSolverOptions lin_opts = {},
                 MPI_Comm comm = MPI_COMM_WORLD);

The nonlinear solver configuration options are provided by the NonlinearSolverOptions struct:

/// Nonlinear solution scheme parameters
struct NonlinearSolverOptions {
  /// Nonlinear solver selection
  NonlinearSolver nonlin_solver = NonlinearSolver::Newton;

  /// Relative tolerance
  double relative_tol = 1.0e-8;

  /// Absolute tolerance
  double absolute_tol = 1.0e-12;

  /// Maximum number of iterations
  int max_iterations = 20;

  /// Debug print level
  int print_level = 0;
};

The current possible nonlinear solution algorithms are:

/// Nonlinear solver method indicator
enum class NonlinearSolver
{
  Newton,                    /**< MFEM-native Newton-Raphson */
  LBFGS,                     /**< MFEM-native Limited memory BFGS */
  KINFullStep,               /**< KINSOL Full Newton (Sundials must be enabled) */
  KINBacktrackingLineSearch, /**< KINSOL Newton with Backtracking Line Search (Sundials must be enabled) */
  KINPicard                  /**< KINSOL Picard (Sundials must be enabled) */
};

The linear solver configuration options are provided by the LinearSolverOptions struct:

/// Parameters for an iterative linear solution scheme
struct LinearSolverOptions {
  /// Linear solver selection
  LinearSolver linear_solver = LinearSolver::GMRES;

  /// PreconditionerOptions selection
  Preconditioner preconditioner = Preconditioner::HypreJacobi;

  /// Relative tolerance
  double relative_tol = 1.0e-8;

  /// Absolute tolerance
  double absolute_tol = 1.0e-12;

  /// Maximum number of iterations
  int max_iterations = 300;

  /// Debugging print level for the linear solver
  int print_level = 0;

  /// Debugging print level for the preconditioner
  int preconditioner_print_level = 0;
};

The current possible linear solution algorithms are:

/// Linear solution method indicator
enum class LinearSolver
{
  CG,       /**< Conjugate gradient */
  GMRES,    /**< Generalized minimal residual method */
  SuperLU,  /**< SuperLU MPI-enabled direct nodal solver */
  Strumpack /**< Strumpack MPI-enabled direct frontal solver*/
};

The current possible preconditioners for iterative linear solvers are:

/// The type of preconditioner to be used
enum class Preconditioner
{
  HypreJacobi,      /**< Hypre-based Jacobi */
  HypreL1Jacobi,    /**< Hypre-based L1-scaled Jacobi */
  HypreGaussSeidel, /**< Hypre-based Gauss-Seidel */
  HypreAMG,         /**< Hypre's BoomerAMG algebraic multi-grid */
  HypreILU,         /**< Hypre's Incomplete LU */
  AMGX,             /**< NVIDIA's AMGX GPU-enabled algebraic multi-grid, GPU builds only */
  None              /**< No preconditioner used */
};

Custom configuration via pointers

If the nonlinear and linear solvers provided by the above options are not sufficient for your application, custom solvers can be written for both the nonlinear and linear solver objects. For this approach, the direct constructor for EquationSolver is used:

  /**
   * Constructs a new nonlinear equation solver
   * @param[in] nonlinear_solver A constructed nonlinear solver
   * @param[in] linear_solver A constructed linear solver to be called by the nonlinear algorithm and adjoint
   * equation solves
   * @param[in] preconditioner An optional constructed precondition to aid the linear solver
   */
  EquationSolver(std::unique_ptr<mfem::NewtonSolver> nonlinear_solver, std::unique_ptr<mfem::Solver> linear_solver,
                 std::unique_ptr<mfem::Solver> preconditioner = nullptr);

The nonlinear and linear solvers are required while the preconditioner is optional.

Although the nonlinear solver is expected to be of type mfem::NewtonSolver, it does not have to be a Newton-type method. This is simply the preferred MFEM container for nonlinear solvers. For example, the included L-BFGS solver is derived from this class. This class is preferred over type mfem::Solver as it has the appropriate methods for checking convergence as needed.

Use within physics modules

An example of configuring a SolidMechanics simulation module via options stucts is below:

  serac::LinearSolverOptions linear_options{.linear_solver = LinearSolver::SuperLU};

  serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver  = NonlinearSolver::Newton,
                                                  .relative_tol   = 1.0e-12,
                                                  .absolute_tol   = 1.0e-12,
                                                  .max_iterations = 5000,
                                                  .print_level    = 1};

  SolidMechanics<p, dim> solid_solver(nonlinear_options, linear_options, solid_mechanics::default_quasistatic_options,
                                      GeometricNonlinearities::Off, "solid_mechanics", mesh_tag);

Alternatively, you can build an EquationSolver using custom nonlinear and linear solvers if it is required by your application. An example of a parameterized SolidMechanics module that uses a custom EquationSolver is below:

  auto nonlinear_solver = std::make_unique<mfem::NewtonSolver>(pmesh.GetComm());
  nonlinear_solver->SetPrintLevel(1);
  nonlinear_solver->SetMaxIter(30);
  nonlinear_solver->SetAbsTol(1.0e-12);
  nonlinear_solver->SetRelTol(1.0e-10);

  auto linear_solver = std::make_unique<mfem::HypreGMRES>(pmesh.GetComm());
  linear_solver->SetPrintLevel(1);
  linear_solver->SetMaxIter(500);
  linear_solver->SetTol(1.0e-6);

  auto preconditioner = std::make_unique<mfem::HypreBoomerAMG>();
  linear_solver->SetPreconditioner(*preconditioner);

  auto equation_solver = std::make_unique<EquationSolver>(std::move(nonlinear_solver), std::move(linear_solver),
                                                          std::move(preconditioner));

  SolidMechanics<p, dim, Parameters<H1<1>, H1<1>>> solid_solver(
      std::move(equation_solver), solid_mechanics::default_quasistatic_options, GeometricNonlinearities::On,
      "parameterized_solid", mesh_tag, {"shear", "bulk"});

Warning

Each physics module must have its own EquationSolver. They cannot be reused between modules.