equation_solver.hpp

Go to the documentation of this file.

// Copyright (c) Lawrence Livermore National Security, LLC and
// other Serac Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include <memory>
#include <optional>
#include <variant>
#include <utility>
 
#include "mpi.h"
#include "mfem.hpp"
 
#include "serac/infrastructure/input.hpp"
#include "serac/numerics/solver_config.hpp"
#include "serac/numerics/petsc_solvers.hpp"
 
namespace serac {
 
class EquationSolver {
 public:
  // _equationsolver_constructor_start
  EquationSolver(std::unique_ptr<mfem::NewtonSolver> nonlinear_solver, std::unique_ptr<mfem::Solver> linear_solver,
                 std::unique_ptr<mfem::Solver> preconditioner = nullptr);
  // _equationsolver_constructor_end
 
  // _build_equationsolver_start
  EquationSolver(NonlinearSolverOptions nonlinear_opts = {}, LinearSolverOptions lin_opts = {},
                 MPI_Comm comm = MPI_COMM_WORLD);
  // _build_equationsolver_end
 
  void setOperator(const mfem::Operator& op);
 
  void solve(mfem::Vector& x) const;
 
  mfem::NewtonSolver& nonlinearSolver() { return *nonlin_solver_; }
 
  const mfem::NewtonSolver& nonlinearSolver() const { return *nonlin_solver_; }
 
  mfem::Solver& linearSolver() { return *lin_solver_; }
 
  const mfem::Solver& linearSolver() const { return *lin_solver_; }
 
  mfem::Solver& preconditioner() { return *preconditioner_; }
 
  const mfem::Solver& preconditioner() const { return *preconditioner_; }
 
  static void defineInputFileSchema(axom::inlet::Container& container);
 
 private:
  std::unique_ptr<mfem::Solver> preconditioner_;
 
  std::unique_ptr<mfem::Solver> lin_solver_;
 
  std::unique_ptr<mfem::NewtonSolver> nonlin_solver_;
 
  bool nonlin_solver_set_solver_called_ = false;
};
 
class SuperLUSolver : public mfem::Solver {
 public:
  SuperLUSolver(int print_level, MPI_Comm comm) : superlu_solver_(comm)
  {
    superlu_solver_.SetColumnPermutation(mfem::superlu::PARMETIS);
    if (print_level == 0) {
      superlu_solver_.SetPrintStatistics(false);
    }
  }
 
  void Mult(const mfem::Vector& input, mfem::Vector& output) const;
 
  void SetOperator(const mfem::Operator& op);
 
 private:
  mutable std::unique_ptr<mfem::SuperLURowLocMatrix> superlu_mat_;
 
  mfem::SuperLUSolver superlu_solver_;
};
 
#ifdef MFEM_USE_STRUMPACK
class StrumpackSolver : public mfem::Solver {
 public:
  StrumpackSolver(int print_level, MPI_Comm comm) : strumpack_solver_(comm)
  {
    strumpack_solver_.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
    strumpack_solver_.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
 
    if (print_level == 1) {
      strumpack_solver_.SetPrintFactorStatistics(true);
      strumpack_solver_.SetPrintSolveStatistics(true);
    }
  }
 
  void Mult(const mfem::Vector& input, mfem::Vector& output) const;
 
  void SetOperator(const mfem::Operator& op);
 
 private:
  mutable std::unique_ptr<mfem::STRUMPACKRowLocMatrix> strumpack_mat_;
 
  mfem::STRUMPACKSolver strumpack_solver_;
};
 
#endif
 
std::unique_ptr<mfem::HypreParMatrix> buildMonolithicMatrix(const mfem::BlockOperator& block_operator);
 
std::unique_ptr<mfem::NewtonSolver> buildNonlinearSolver(const NonlinearSolverOptions& nonlinear_opts,
                                                         const LinearSolverOptions& linear_opts,
                                                         mfem::Solver& preconditioner, MPI_Comm comm = MPI_COMM_WORLD);
 
std::pair<std::unique_ptr<mfem::Solver>, std::unique_ptr<mfem::Solver>> buildLinearSolverAndPreconditioner(
    LinearSolverOptions linear_opts = {}, MPI_Comm comm = MPI_COMM_WORLD);
 
std::unique_ptr<mfem::Solver> buildPreconditioner(LinearSolverOptions linear_opts,
                                                  [[maybe_unused]] MPI_Comm comm = MPI_COMM_WORLD);
 
#ifdef MFEM_USE_AMGX
std::unique_ptr<mfem::AmgXSolver> buildAMGX(const AMGXOptions& options, const MPI_Comm comm);
#endif
 
}  // namespace serac
 
template <>
struct FromInlet<serac::LinearSolverOptions> {
  serac::LinearSolverOptions operator()(const axom::inlet::Container& base);
};
 
template <>
struct FromInlet<serac::NonlinearSolverOptions> {
  serac::NonlinearSolverOptions operator()(const axom::inlet::Container& base);
};
 
template <>
struct FromInlet<serac::EquationSolver> {
  serac::EquationSolver operator()(const axom::inlet::Container& base);
};