Move solvers into ProblemOperator classes

src/executioners/steady_executioner.cpp

@@ -13,7 +13,7 @@ SteadyExecutioner::Solve() const
 {
   // Advance time step.
   _problem->_preprocessors.Solve();
-  _problem->GetOperator()->Solve(*(_problem->_f));
+  _problem->GetOperator()->Solve();
   _problem->_postprocessors.Solve();
 
   // Output data

src/executioners/transient_executioner.cpp

@@ -27,7 +27,7 @@ TransientExecutioner::Step(double dt, int it) const
 
   // Advance time step.
   _problem->_preprocessors.Solve(_t);
-  _problem->_ode_solver->Step(*(_problem->_f), _t, dt);
+  _problem->GetOperator()->Step(_t, dt);
   _problem->_postprocessors.Solve(_t);
 
   // Output data

src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp

@@ -218,8 +218,8 @@ ComplexMaxwellOperator::Solve(mfem::Vector & X)
 
   auto * jac_z = jac.As<mfem::ComplexHypreParMatrix>()->GetSystemMatrix();
 
-  _problem._jacobian_solver->SetOperator(*jac_z);
-  _problem._jacobian_solver->Mult(rhs, u);
+  _jacobian_solver->SetOperator(*jac_z);
+  _jacobian_solver->Mult(rhs, u);
   sqlf.RecoverFEMSolution(u, lf, *_u);
 
   _problem._gridfunctions.GetRef(_trial_var_names.at(0)) = _u->real();

src/formulations/Dual/dual_formulation.cpp

@@ -59,20 +59,16 @@ DualFormulation::DualFormulation(std::string alpha_coef_name,
 }
 
 void
-DualFormulation::ConstructJacobianPreconditioner()
+DualFormulation::ConstructJacobianSolver()
 {
   auto precond =
-      std::make_shared<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
+      std::make_unique<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
 
   precond->SetSingularProblem();
   precond->SetPrintLevel(-1);
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-DualFormulation::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000});
 }
 
@@ -238,12 +234,12 @@ DualOperator::ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vecto
   logger.info("{} ImplicitSolve: {} seconds", typeid(this).name(), sw);
 }
 
-void
-DualOperator::UpdateOffsets()
+int
+DualOperator::GetSolutionVectorSize() const
 {
   // Blocks for solution vector are smaller than the operator size for DualOperator,
   // as curl is stored separately. Block operator only has the HCurl TrueVSize;
-  return UpdateOffsetsWithSize(_trial_variables.size() - 1);
+  return static_cast<int>(GetTrialVariablesSize() - 1);
 }
 
 } // namespace hephaestus

src/formulations/Dual/dual_formulation.hpp

@@ -16,8 +16,6 @@ class DualFormulation : public TimeDomainEMFormulation
 
   ~DualFormulation() override = default;
 
-  void ConstructJacobianPreconditioner() override;
-
   void ConstructJacobianSolver() override;
 
   void ConstructOperator() override;
@@ -74,7 +72,7 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator
   mfem::ParGridFunction * _dv{nullptr}; // HDiv vector field
 
 protected:
-  void UpdateOffsets() override;
+  int GetSolutionVectorSize() const override;
 
   std::unique_ptr<mfem::ParDiscreteLinearOperator> _curl;
 };

src/formulations/HCurl/hcurl_formulation.cpp

@@ -67,20 +67,16 @@ HCurlFormulation::ConstructOperator()
 }
 
 void
-HCurlFormulation::ConstructJacobianPreconditioner()
+HCurlFormulation::ConstructJacobianSolver()
 {
   auto precond =
-      std::make_shared<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
+      std::make_unique<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
 
   precond->SetSingularProblem();
   precond->SetPrintLevel(-1);
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-HCurlFormulation::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG);
 }
 

src/formulations/HCurl/hcurl_formulation.hpp

@@ -18,8 +18,6 @@ class HCurlFormulation : public TimeDomainEMFormulation
 
   void ConstructOperator() override;
 
-  void ConstructJacobianPreconditioner() override;
-
   void ConstructJacobianSolver() override;
 
   void RegisterGridFunctions() override;

src/formulations/Magnetostatic/statics_formulation.cpp

@@ -38,20 +38,16 @@ StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string
 }
 
 void
-StaticsFormulation::ConstructJacobianPreconditioner()
+StaticsFormulation::ConstructJacobianSolver()
 {
-  std::shared_ptr<mfem::HypreAMS> precond{std::make_shared<mfem::HypreAMS>(
-      GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace())};
+  auto precond = std::make_unique<mfem::HypreAMS>(
+      GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace());
 
   precond->SetSingularProblem();
   precond->SetPrintLevel(-1);
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-StaticsFormulation::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES,
                                      {._max_iteration = 100, ._k_dim = 10});
 }
@@ -152,8 +148,8 @@ StaticsOperator::Solve(mfem::Vector & X)
 
   // Define and apply a parallel FGMRES solver for AX=B with the AMS
   // preconditioner from hypre.
-  _problem._jacobian_solver->SetOperator(curl_mu_inv_curl);
-  _problem._jacobian_solver->Mult(rhs_tdofs, sol_tdofs);
+  _jacobian_solver->SetOperator(curl_mu_inv_curl);
+  _jacobian_solver->Mult(rhs_tdofs, sol_tdofs);
   blf.RecoverFEMSolution(sol_tdofs, lf, gf);
 
   logger.info("{} Solve: {} seconds", typeid(this).name(), sw);

src/formulations/Magnetostatic/statics_formulation.hpp

@@ -14,8 +14,6 @@ class StaticsFormulation : public SteadyStateEMFormulation
 
   ~StaticsFormulation() override = default;
 
-  void ConstructJacobianPreconditioner() override;
-
   void ConstructJacobianSolver() override;
 
   void ConstructOperator() override;

src/problem_builders/problem_builder_base.cpp

@@ -3,13 +3,6 @@
 namespace hephaestus
 {
 
-Problem::~Problem()
-{
-  // Ensure that all owned memory is properly freed!
-  _f.reset();
-  _ode_solver.reset();
-}
-
 void
 Problem::Update()
 {
@@ -97,15 +90,15 @@ ProblemBuilder::SetSolverOptions(hephaestus::InputParameters & solver_options)
 }
 
 void
-ProblemBuilder::SetJacobianPreconditioner(std::shared_ptr<mfem::Solver> preconditioner)
+ProblemBuilder::SetJacobianPreconditioner(std::unique_ptr<mfem::Solver> preconditioner)
 {
-  GetProblem()->_jacobian_preconditioner = preconditioner;
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(preconditioner));
 }
 
 void
-ProblemBuilder::SetJacobianSolver(std::shared_ptr<mfem::Solver> jacobian_solver)
+ProblemBuilder::SetJacobianSolver(std::unique_ptr<mfem::Solver> jacobian_solver)
 {
-  GetProblem()->_jacobian_solver = jacobian_solver;
+  GetProblem()->GetOperator()->SetJacobianSolver(std::move(jacobian_solver));
 }
 
 void
@@ -223,26 +216,16 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr<hephaestus::S
 }
 
 void
-ProblemBuilder::ConstructJacobianPreconditioner()
+ProblemBuilder::ConstructJacobianSolver()
 {
-  auto precond = std::make_shared<mfem::HypreBoomerAMG>();
+  auto precond = std::make_unique<mfem::HypreBoomerAMG>();
   precond->SetPrintLevel(GetGlobalPrintLevel());
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-ProblemBuilder::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_GMRES);
 }
 
-void
-ProblemBuilder::ConstructBlockVector()
-{
-  GetProblem()->_f = std::make_unique<mfem::BlockVector>();
-}
-
 void
 ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams default_params)
 {
@@ -258,14 +241,13 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       solver_options.GetOptionalParam<int>("PrintLevel", default_params._print_level);
   const auto k_dim = solver_options.GetOptionalParam<unsigned int>("KDim", default_params._k_dim);
 
-  auto preconditioner =
-      std::dynamic_pointer_cast<mfem::HypreSolver>(GetProblem()->_jacobian_preconditioner);
+  auto preconditioner = GetProblem()->GetOperator()->JacobianPreconditioner<mfem::HypreSolver>();
 
   switch (type)
   {
     case SolverType::HYPRE_PCG:
     {
-      auto solver = std::make_shared<mfem::HyprePCG>(GetProblem()->_comm);
+      auto solver = std::make_unique<mfem::HyprePCG>(GetProblem()->_comm);
 
       solver->SetTol(tolerance);
       solver->SetAbsTol(abs_tolerance);
@@ -275,12 +257,12 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       if (preconditioner)
         solver->SetPreconditioner(*preconditioner);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::HYPRE_GMRES:
     {
-      auto solver = std::make_shared<mfem::HypreGMRES>(GetProblem()->_comm);
+      auto solver = std::make_unique<mfem::HypreGMRES>(GetProblem()->_comm);
 
       solver->SetTol(tolerance);
       solver->SetAbsTol(abs_tolerance);
@@ -291,12 +273,12 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       if (preconditioner)
         solver->SetPreconditioner(*preconditioner);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::HYPRE_FGMRES:
     {
-      auto solver = std::make_shared<mfem::HypreFGMRES>(GetProblem()->_comm);
+      auto solver = std::make_unique<mfem::HypreFGMRES>(GetProblem()->_comm);
 
       solver->SetTol(tolerance);
       solver->SetMaxIter(max_iter);
@@ -306,25 +288,25 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       if (preconditioner)
         solver->SetPreconditioner(*preconditioner);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::HYPRE_AMG:
     {
-      auto solver = std::make_shared<mfem::HypreBoomerAMG>();
+      auto solver = std::make_unique<mfem::HypreBoomerAMG>();
 
       solver->SetTol(tolerance);
       solver->SetMaxIter(max_iter);
       solver->SetPrintLevel(print_level);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::SUPER_LU:
     {
-      auto solver = std::make_shared<hephaestus::SuperLUSolver>(GetProblem()->_comm);
+      auto solver = std::make_unique<hephaestus::SuperLUSolver>(GetProblem()->_comm);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     default:
@@ -338,14 +320,14 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
 void
 ProblemBuilder::ConstructNonlinearSolver()
 {
-  auto nl_solver = std::make_shared<mfem::NewtonSolver>(GetProblem()->_comm);
+  auto nl_solver = std::make_unique<mfem::NewtonSolver>(GetProblem()->_comm);
 
   // Defaults to one iteration, without further nonlinear iterations
   nl_solver->SetRelTol(0.0);
   nl_solver->SetAbsTol(0.0);
   nl_solver->SetMaxIter(1);
 
-  GetProblem()->_nonlinear_solver = nl_solver;
+  GetProblem()->GetOperator()->SetNonlinearSolver(std::move(nl_solver));
 }
 
 void
@@ -390,13 +372,10 @@ ProblemBuilder::FinalizeProblem(bool build_operator)
     ConstructOperator();
   }
 
-  ConstructBlockVector();
-
   InitializeAuxSolvers();
   InitializeSources();
   InitializeOperator();
 
-  ConstructJacobianPreconditioner();
   ConstructJacobianSolver();
   ConstructNonlinearSolver();
 

src/problem_builders/problem_builder_base.hpp

@@ -12,14 +12,14 @@ namespace hephaestus
 {
 
 /// Forwards declaration.
-class ProblemOperatorInterface;
+class ProblemOperatorBase;
 
 /// Base problem class.
 class Problem
 {
 public:
   Problem() = default;
-  virtual ~Problem();
+  virtual ~Problem() = default;
 
   std::shared_ptr<mfem::ParMesh> _pmesh{nullptr};
   hephaestus::BCMap _bc_map;
@@ -30,13 +30,6 @@ class Problem
   hephaestus::Outputs _outputs;
   hephaestus::InputParameters _solver_options;
 
-  std::unique_ptr<mfem::ODESolver> _ode_solver{nullptr};
-  std::unique_ptr<mfem::BlockVector> _f{nullptr};
-
-  std::shared_ptr<mfem::Solver> _jacobian_preconditioner{nullptr};
-  std::shared_ptr<mfem::Solver> _jacobian_solver{nullptr};
-  std::shared_ptr<mfem::NewtonSolver> _nonlinear_solver{nullptr};
-
   hephaestus::FECollections _fecs;
   hephaestus::FESpaces _fespaces;
   hephaestus::GridFunctions _gridfunctions;
@@ -46,7 +39,7 @@ class Problem
   int _num_procs;
 
   /// Returns a pointer to the operator. See derived classes.
-  [[nodiscard]] virtual hephaestus::ProblemOperatorInterface * GetOperator() const = 0;
+  [[nodiscard]] virtual hephaestus::ProblemOperatorBase * GetOperator() const = 0;
 
   /// Virtual method to construct the operator. Call for default problems.
   virtual void ConstructOperator() = 0;
@@ -80,8 +73,8 @@ class ProblemBuilder
   void SetSources(hephaestus::Sources & sources);
   void SetOutputs(hephaestus::Outputs & outputs);
   void SetSolverOptions(hephaestus::InputParameters & solver_options);
-  void SetJacobianPreconditioner(std::shared_ptr<mfem::Solver> preconditioner);
-  void SetJacobianSolver(std::shared_ptr<mfem::Solver> solver);
+  void SetJacobianPreconditioner(std::unique_ptr<mfem::Solver> preconditioner);
+  void SetJacobianSolver(std::unique_ptr<mfem::Solver> solver);
   void SetCoefficients(hephaestus::Coefficients & coefficients);
 
   void AddFESpace(std::string fespace_name,
@@ -100,10 +93,8 @@ class ProblemBuilder
   virtual void RegisterAuxSolvers() = 0;
   virtual void RegisterCoefficients() = 0;
 
-  virtual void ConstructJacobianPreconditioner();
   virtual void ConstructJacobianSolver();
   virtual void ConstructNonlinearSolver();
-  virtual void ConstructBlockVector();
 
   virtual void ConstructOperator() = 0;
   virtual void ConstructTimestepper() = 0;