From b40e516aeabc535a1130c2051c201dd3183afa67 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 13:53:28 +0000 Subject: [PATCH 01/30] Added a ConstructJacobianSolver method to the ProblemOperator classes. Problem operators now have their own specific preconditioners and solvers. --- .../complex_maxwell_formulation.cpp | 8 +- .../complex_maxwell_formulation.hpp | 2 + src/formulations/Dual/dual_formulation.cpp | 30 ++++-- src/formulations/Dual/dual_formulation.hpp | 2 + src/formulations/HCurl/hcurl_formulation.cpp | 32 ++++-- src/formulations/HCurl/hcurl_formulation.hpp | 11 +++ .../Magnetostatic/statics_formulation.cpp | 32 ++++-- .../Magnetostatic/statics_formulation.hpp | 2 + src/problem_builders/problem_builder_base.cpp | 98 +------------------ src/problem_builders/problem_builder_base.hpp | 33 ------- .../problem_operator_base.cpp | 19 ++++ .../problem_operator_base.hpp | 3 + 12 files changed, 112 insertions(+), 160 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp index 238ca0c36..bdcecc983 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp @@ -27,7 +27,7 @@ ComplexMaxwellFormulation::ComplexMaxwellFormulation(std::string alpha_coef_name void ComplexMaxwellFormulation::ConstructJacobianSolver() { - ConstructJacobianSolverWithOptions(SolverType::SUPER_LU); + GetProblem()->GetOperator()->ConstructJacobianSolver(); } void @@ -226,4 +226,10 @@ ComplexMaxwellOperator::Solve(mfem::Vector & X) _problem._gridfunctions.GetRef(_trial_var_names.at(1)) = _u->imag(); } +void +ComplexMaxwellOperator::ConstructJacobianSolver() +{ + _jacobian_solver = std::make_unique(_problem._comm); +} + } // namespace hephaestus diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp index 5246fa18c..9ecce6d11 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp @@ -77,6 +77,8 @@ class ComplexMaxwellOperator : public ProblemOperator void Update() override; void Solve(mfem::Vector & X) override; + void ConstructJacobianSolver() override; + std::string _h_curl_var_complex_name, _h_curl_var_real_name, _h_curl_var_imag_name, _stiffness_coef_name, _mass_coef_name, _loss_coef_name; diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 56a21370c..abbd66668 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -61,15 +61,7 @@ DualFormulation::DualFormulation(std::string alpha_coef_name, void DualFormulation::ConstructJacobianSolver() { - auto precond = - std::make_unique(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - - ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000}); + GetProblem()->GetOperator()->ConstructJacobianSolver(); } void @@ -219,6 +211,26 @@ DualOperator::Update() _curl->Assemble(); } +void +DualOperator::ConstructJacobianSolver() +{ + auto precond = std::make_unique(GetEquationSystem()->_test_pfespaces.at(0)); + + precond->SetSingularProblem(); + precond->SetPrintLevel(-1); + + auto solver = std::make_unique(_problem._comm); + + solver->SetTol(1e-16); + solver->SetAbsTol(1e-16); + solver->SetMaxIter(1000); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); + + _jacobian_preconditioner = std::move(precond); + _jacobian_solver = std::move(solver); +} + void DualOperator::ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) { diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index 5a09f8923..11ab5bddc 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -61,6 +61,8 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator void Update() override; + void ConstructJacobianSolver() override; + void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override; mfem::ParFiniteElementSpace * _h_curl_fe_space{nullptr}; diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index bdb5acee0..2159993ba 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -62,22 +62,14 @@ HCurlFormulation::ConstructOperator() auto equation_system = std::make_unique(weak_form_params); - GetProblem()->SetOperator(std::make_unique( + GetProblem()->SetOperator(std::make_unique( *GetProblem(), std::move(equation_system))); } void HCurlFormulation::ConstructJacobianSolver() { - auto precond = - std::make_unique(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - - ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG); + GetProblem()->GetOperator()->ConstructJacobianSolver(); } void @@ -102,6 +94,26 @@ HCurlFormulation::RegisterGridFunctions() TimeDomainEquationSystemProblemBuilder::RegisterGridFunctions(); }; +void +HCurlProblemOperator::ConstructJacobianSolver() +{ + auto precond = std::make_unique(GetEquationSystem()->_test_pfespaces.at(0)); + + precond->SetSingularProblem(); + precond->SetPrintLevel(-1); + + auto solver = std::make_unique(_problem._comm); + + solver->SetTol(1e-16); + solver->SetAbsTol(1e-16); + solver->SetMaxIter(1000); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); + + _jacobian_preconditioner = std::move(precond); + _jacobian_solver = std::move(solver); +} + CurlCurlEquationSystem::CurlCurlEquationSystem(const hephaestus::InputParameters & params) : _h_curl_var_name(params.GetParam("HCurlVarName")), _alpha_coef_name(params.GetParam("AlphaCoefName")), diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index dc9e491a1..69a17602b 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -30,6 +30,17 @@ class HCurlFormulation : public TimeDomainEMFormulation const std::string _h_curl_var_name; }; +class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator +{ +public: + HCurlProblemOperator(hephaestus::Problem & problem, + std::unique_ptr equation_system) + : TimeDomainEquationSystemProblemOperator(problem, std::move(equation_system)) + { + } + void ConstructJacobianSolver() override; +}; + class CurlCurlEquationSystem : public TimeDependentEquationSystem { public: diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 25ddd21e5..cb61aa75c 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -40,16 +40,7 @@ StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string void StaticsFormulation::ConstructJacobianSolver() { - auto precond = std::make_unique( - GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace()); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - - ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES, - {._max_iteration = 100, ._k_dim = 10}); + GetProblem()->GetOperator()->ConstructJacobianSolver(); } void @@ -113,6 +104,27 @@ StaticsOperator::Init() ProblemOperator::Init(); } +void +StaticsOperator::ConstructJacobianSolver() +{ + auto precond = + std::make_unique(_problem._gridfunctions.Get(_h_curl_var_name)->ParFESpace()); + + precond->SetSingularProblem(); + precond->SetPrintLevel(-1); + + auto solver = std::make_unique(_problem._comm); + + solver->SetTol(1e-16); + solver->SetMaxIter(100); + solver->SetKDim(10); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); + + _jacobian_preconditioner = std::move(precond); + _jacobian_solver = std::move(solver); +} + /* This is the main method that solves for u. diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index b0a901695..5f20aa06c 100644 --- a/src/formulations/Magnetostatic/statics_formulation.hpp +++ b/src/formulations/Magnetostatic/statics_formulation.hpp @@ -40,6 +40,8 @@ class StaticsOperator : public ProblemOperator void Init() override; void Solve(mfem::Vector & X) override; + void ConstructJacobianSolver() override; + private: std::string _h_curl_var_name, _stiffness_coef_name; diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index d5e4ef8dd..8ecc27b20 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -218,103 +218,7 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr(); - precond->SetPrintLevel(GetGlobalPrintLevel()); - - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - - ConstructJacobianSolverWithOptions(SolverType::HYPRE_GMRES); -} - -void -ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams default_params) -{ - const auto & solver_options = GetProblem()->_solver_options; - - const auto tolerance = - solver_options.GetOptionalParam("Tolerance", default_params._tolerance); - const auto abs_tolerance = - solver_options.GetOptionalParam("AbsTolerance", default_params._abs_tolerance); - const auto max_iter = - solver_options.GetOptionalParam("MaxIter", default_params._max_iteration); - const auto print_level = - solver_options.GetOptionalParam("PrintLevel", default_params._print_level); - const auto k_dim = solver_options.GetOptionalParam("KDim", default_params._k_dim); - - auto preconditioner = GetProblem()->GetOperator()->JacobianPreconditioner(); - - switch (type) - { - case SolverType::HYPRE_PCG: - { - auto solver = std::make_unique(GetProblem()->_comm); - - solver->SetTol(tolerance); - solver->SetAbsTol(abs_tolerance); - solver->SetMaxIter(max_iter); - solver->SetPrintLevel(print_level); - - if (preconditioner) - solver->SetPreconditioner(*preconditioner); - - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - case SolverType::HYPRE_GMRES: - { - auto solver = std::make_unique(GetProblem()->_comm); - - solver->SetTol(tolerance); - solver->SetAbsTol(abs_tolerance); - solver->SetMaxIter(max_iter); - solver->SetKDim(k_dim); - solver->SetPrintLevel(print_level); - - if (preconditioner) - solver->SetPreconditioner(*preconditioner); - - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - case SolverType::HYPRE_FGMRES: - { - auto solver = std::make_unique(GetProblem()->_comm); - - solver->SetTol(tolerance); - solver->SetMaxIter(max_iter); - solver->SetKDim(k_dim); - solver->SetPrintLevel(print_level); - - if (preconditioner) - solver->SetPreconditioner(*preconditioner); - - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - case SolverType::HYPRE_AMG: - { - auto solver = std::make_unique(); - - solver->SetTol(tolerance); - solver->SetMaxIter(max_iter); - solver->SetPrintLevel(print_level); - - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - case SolverType::SUPER_LU: - { - auto solver = std::make_unique(GetProblem()->_comm); - - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - default: - { - MFEM_ABORT("Unsupported solver type specified."); - break; - } - } + GetProblem()->GetOperator()->ConstructJacobianSolver(); } void diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 39d85ac62..64a4a1979 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -115,39 +115,6 @@ class ProblemBuilder /// Protected constructor. Derived classes must call this constructor. ProblemBuilder(hephaestus::Problem * problem) : _problem{problem} {} - /// Supported Jacobian solver types. - enum class SolverType - { - HYPRE_PCG, - HYPRE_GMRES, - HYPRE_FGMRES, - HYPRE_AMG, - SUPER_LU - }; - - /// Structure containing default parameters which can be passed to @a ConstructJacobianSolverWithOptions. - /// These will be used if the user has not supplied their own values. - struct SolverParams - { - double _tolerance; - double _abs_tolerance; - - unsigned int _max_iteration; - - int _print_level; - int _k_dim; - }; - - /// Called in @a ConstructJacobianSolver. This will create a solver of the chosen type and use the user's input - /// parameters if they have been provided. - void ConstructJacobianSolverWithOptions(SolverType type, - SolverParams default_params = { - ._tolerance = 1e-16, - ._abs_tolerance = 1e-16, - ._max_iteration = 1000, - ._print_level = GetGlobalPrintLevel(), - ._k_dim = 10}); - /// Overridden in derived classes. [[nodiscard]] virtual hephaestus::Problem * GetProblem() const = 0; diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index 55aeeb746..b6c8e5331 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -8,6 +8,25 @@ ProblemOperatorBase::ProblemOperatorBase(hephaestus::Problem & problem) : _probl _block_vector = std::make_unique(); } +void +ProblemOperatorBase::ConstructJacobianSolver() +{ + auto precond = std::make_unique(); + precond->SetPrintLevel(GetGlobalPrintLevel()); + + auto solver = std::make_unique(_problem._comm); + + solver->SetTol(1e-16); + solver->SetAbsTol(1e-16); + solver->SetMaxIter(1000); + solver->SetKDim(10); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); + + _jacobian_preconditioner = std::move(precond); + _jacobian_solver = std::move(solver); +} + void ProblemOperatorBase::SetTrialVariables() { diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 323c373ab..8afc5f1d5 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -44,6 +44,9 @@ class ProblemOperatorBase return static_cast(_jacobian_preconditioner.get()); } + /// Override in derived classes to construct the Jacobian solver. + virtual void ConstructJacobianSolver(); + protected: /// Use of protected constructor to only allow construction by derived classes. /// All problem operator classes are built on-top of this class and it should not From 2d7c18a232e54cc94236508a305e501bc1ec40f6 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 14:00:34 +0000 Subject: [PATCH 02/30] Removed ConstructJacobianSolver method from problem operator classes; ConstructJacobianSolver method in problem operator classes is now called in the constructor of the base class. --- .../ComplexMaxwell/complex_maxwell_formulation.cpp | 6 ------ .../ComplexMaxwell/complex_maxwell_formulation.hpp | 2 -- src/formulations/Dual/dual_formulation.cpp | 6 ------ src/formulations/Dual/dual_formulation.hpp | 2 -- src/formulations/HCurl/hcurl_formulation.cpp | 6 ------ src/formulations/HCurl/hcurl_formulation.hpp | 2 -- src/formulations/Magnetostatic/statics_formulation.cpp | 6 ------ src/formulations/Magnetostatic/statics_formulation.hpp | 2 -- src/problem_builders/problem_builder_base.cpp | 7 ------- src/problem_builders/problem_builder_base.hpp | 1 - src/problem_operators/problem_operator_base.cpp | 8 ++++---- src/problem_operators/problem_operator_base.hpp | 7 ++++--- 12 files changed, 8 insertions(+), 47 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp index bdcecc983..3b3bc56bd 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp @@ -24,12 +24,6 @@ ComplexMaxwellFormulation::ComplexMaxwellFormulation(std::string alpha_coef_name { } -void -ComplexMaxwellFormulation::ConstructJacobianSolver() -{ - GetProblem()->GetOperator()->ConstructJacobianSolver(); -} - void ComplexMaxwellFormulation::ConstructOperator() { diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp index 9ecce6d11..7f8dba130 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp @@ -38,8 +38,6 @@ class ComplexMaxwellFormulation : public hephaestus::FrequencyDomainEMFormulatio ~ComplexMaxwellFormulation() override = default; - void ConstructJacobianSolver() override; - void ConstructOperator() override; void RegisterGridFunctions() override; diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index abbd66668..7bea7a5d5 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -58,12 +58,6 @@ DualFormulation::DualFormulation(std::string alpha_coef_name, { } -void -DualFormulation::ConstructJacobianSolver() -{ - GetProblem()->GetOperator()->ConstructJacobianSolver(); -} - void DualFormulation::ConstructOperator() { diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index 11ab5bddc..2db78eaef 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -16,8 +16,6 @@ class DualFormulation : public TimeDomainEMFormulation ~DualFormulation() override = default; - void ConstructJacobianSolver() override; - void ConstructOperator() override; void RegisterGridFunctions() override; diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 2159993ba..26eeca4d8 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -66,12 +66,6 @@ HCurlFormulation::ConstructOperator() *GetProblem(), std::move(equation_system))); } -void -HCurlFormulation::ConstructJacobianSolver() -{ - GetProblem()->GetOperator()->ConstructJacobianSolver(); -} - void HCurlFormulation::RegisterGridFunctions() { diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index 69a17602b..c7a11f782 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -18,8 +18,6 @@ class HCurlFormulation : public TimeDomainEMFormulation void ConstructOperator() override; - void ConstructJacobianSolver() override; - void RegisterGridFunctions() override; void RegisterCoefficients() override; diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index cb61aa75c..818ac5ae7 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -37,12 +37,6 @@ StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string { } -void -StaticsFormulation::ConstructJacobianSolver() -{ - GetProblem()->GetOperator()->ConstructJacobianSolver(); -} - void StaticsFormulation::ConstructOperator() { diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index 5f20aa06c..d336430c1 100644 --- a/src/formulations/Magnetostatic/statics_formulation.hpp +++ b/src/formulations/Magnetostatic/statics_formulation.hpp @@ -14,8 +14,6 @@ class StaticsFormulation : public SteadyStateEMFormulation ~StaticsFormulation() override = default; - void ConstructJacobianSolver() override; - void ConstructOperator() override; void RegisterGridFunctions() override; diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 8ecc27b20..4d0468103 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -215,12 +215,6 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr_sources.Register(source_name, std::move(source)); } -void -ProblemBuilder::ConstructJacobianSolver() -{ - GetProblem()->GetOperator()->ConstructJacobianSolver(); -} - void ProblemBuilder::ConstructNonlinearSolver() { @@ -280,7 +274,6 @@ ProblemBuilder::FinalizeProblem(bool build_operator) InitializeSources(); InitializeOperator(); - ConstructJacobianSolver(); ConstructNonlinearSolver(); ConstructTimestepper(); diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 64a4a1979..76e9035d9 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -93,7 +93,6 @@ class ProblemBuilder virtual void RegisterAuxSolvers() = 0; virtual void RegisterCoefficients() = 0; - virtual void ConstructJacobianSolver(); virtual void ConstructNonlinearSolver(); virtual void ConstructOperator() = 0; diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index b6c8e5331..bc7fc34b2 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -3,10 +3,7 @@ namespace hephaestus { -ProblemOperatorBase::ProblemOperatorBase(hephaestus::Problem & problem) : _problem{problem} -{ - _block_vector = std::make_unique(); -} +ProblemOperatorBase::ProblemOperatorBase(hephaestus::Problem & problem) : _problem{problem} {} void ProblemOperatorBase::ConstructJacobianSolver() @@ -100,6 +97,9 @@ ProblemOperatorBase::UpdateBlockVector(mfem::BlockVector & X) void ProblemOperatorBase::Init() { + _block_vector = std::make_unique(); + ConstructJacobianSolver(); + SetTrialVariables(); UpdateOffsets(); diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 8afc5f1d5..8ebbd164f 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -44,15 +44,16 @@ class ProblemOperatorBase return static_cast(_jacobian_preconditioner.get()); } - /// Override in derived classes to construct the Jacobian solver. - virtual void ConstructJacobianSolver(); - protected: /// Use of protected constructor to only allow construction by derived classes. /// All problem operator classes are built on-top of this class and it should not /// be possible to use directly. explicit ProblemOperatorBase(hephaestus::Problem & problem); + /// Override in derived classes to construct the Jacobian solver. Called in the + /// constructor. + virtual void ConstructJacobianSolver(); + /// Set trial variables names. Override in derived classes. virtual void SetTrialVariableNames() {} From 670fd60d5bc18b95e8b679afc1374359faa4e21b Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 14:02:58 +0000 Subject: [PATCH 03/30] Removes SetJacobianPreconditioner and SetJacobianSolver methods from problem builder classes. --- src/problem_builders/problem_builder_base.cpp | 12 ------------ src/problem_builders/problem_builder_base.hpp | 2 -- src/problem_operators/problem_operator_base.hpp | 12 ------------ 3 files changed, 26 deletions(-) diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 4d0468103..b8097127e 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -89,18 +89,6 @@ ProblemBuilder::SetSolverOptions(hephaestus::InputParameters & solver_options) GetProblem()->_solver_options = solver_options; } -void -ProblemBuilder::SetJacobianPreconditioner(std::unique_ptr preconditioner) -{ - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(preconditioner)); -} - -void -ProblemBuilder::SetJacobianSolver(std::unique_ptr jacobian_solver) -{ - GetProblem()->GetOperator()->SetJacobianSolver(std::move(jacobian_solver)); -} - void ProblemBuilder::SetCoefficients(hephaestus::Coefficients & coefficients) { diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 76e9035d9..c831cc048 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -73,8 +73,6 @@ class ProblemBuilder void SetSources(hephaestus::Sources & sources); void SetOutputs(hephaestus::Outputs & outputs); void SetSolverOptions(hephaestus::InputParameters & solver_options); - void SetJacobianPreconditioner(std::unique_ptr preconditioner); - void SetJacobianSolver(std::unique_ptr solver); void SetCoefficients(hephaestus::Coefficients & coefficients); void AddFESpace(std::string fespace_name, diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 8ebbd164f..92b55fe12 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,18 +19,6 @@ class ProblemOperatorBase /// Update the problem operator after a mesh change. virtual void Update(); - /// Set the Jacobian preconditioner. - void SetJacobianPreconditioner(std::unique_ptr preconditioner) - { - _jacobian_preconditioner = std::move(preconditioner); - } - - /// Set the Jacobian solver. - void SetJacobianSolver(std::unique_ptr solver) - { - _jacobian_solver = std::move(solver); - } - /// Set the nonlinear solver. void SetNonlinearSolver(std::unique_ptr nl_solver) { From a17308e12636cc62477a5bcc6cd0cb17652e5878 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 14:18:11 +0000 Subject: [PATCH 04/30] Removes ConstructTimestepper method from problem builder classes. Adds ConstructTimestepper to TimeDomainProblemOperator class; called in constructor. --- src/problem_builders/problem_builder_base.cpp | 1 - src/problem_builders/problem_builder_base.hpp | 1 - .../steady_state_problem_builder.hpp | 2 -- .../time_domain_problem_builder.cpp | 9 ------- .../time_domain_problem_builder.hpp | 2 -- .../time_domain_problem_operator.cpp | 25 ++++++++++++++++--- .../time_domain_problem_operator.hpp | 13 +++++----- 7 files changed, 27 insertions(+), 26 deletions(-) diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index b8097127e..af4302aa6 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -264,7 +264,6 @@ ProblemBuilder::FinalizeProblem(bool build_operator) ConstructNonlinearSolver(); - ConstructTimestepper(); InitializeOutputs(); } diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index c831cc048..3edec7ca3 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -94,7 +94,6 @@ class ProblemBuilder virtual void ConstructNonlinearSolver(); virtual void ConstructOperator() = 0; - virtual void ConstructTimestepper() = 0; virtual void InitializeSources(); diff --git a/src/problem_builders/steady_state_problem_builder.hpp b/src/problem_builders/steady_state_problem_builder.hpp index 994676aff..731caf314 100644 --- a/src/problem_builders/steady_state_problem_builder.hpp +++ b/src/problem_builders/steady_state_problem_builder.hpp @@ -56,8 +56,6 @@ class SteadyStateProblemBuilder : public ProblemBuilder void ConstructOperator() override; - void ConstructTimestepper() override {} - protected: // NB: constructor for derived classes. SteadyStateProblemBuilder(hephaestus::SteadyStateProblem * problem) : ProblemBuilder(problem) {} diff --git a/src/problem_builders/time_domain_problem_builder.cpp b/src/problem_builders/time_domain_problem_builder.cpp index 916857a8f..93fee3693 100644 --- a/src/problem_builders/time_domain_problem_builder.cpp +++ b/src/problem_builders/time_domain_problem_builder.cpp @@ -45,13 +45,4 @@ TimeDomainProblemBuilder::InitializeOperator() GetProblem()->GetOperator()->SetTime(0.0); } -void -TimeDomainProblemBuilder::ConstructTimestepper() -{ - auto ode_solver = std::make_unique(); - ode_solver->Init(*GetProblem()->GetOperator()); - - GetProblem()->GetOperator()->SetODESolver(std::move(ode_solver)); -} - } // namespace hephaestus diff --git a/src/problem_builders/time_domain_problem_builder.hpp b/src/problem_builders/time_domain_problem_builder.hpp index 231b20586..0134c339f 100644 --- a/src/problem_builders/time_domain_problem_builder.hpp +++ b/src/problem_builders/time_domain_problem_builder.hpp @@ -61,8 +61,6 @@ class TimeDomainProblemBuilder : public ProblemBuilder void InitializeOperator() override; - void ConstructTimestepper() override; - protected: /// NB: constructor called in derived classes. TimeDomainProblemBuilder(hephaestus::TimeDomainProblem * problem) : ProblemBuilder(problem) {} diff --git a/src/problem_operators/time_domain_problem_operator.cpp b/src/problem_operators/time_domain_problem_operator.cpp index a4ffa7211..2515b9e21 100644 --- a/src/problem_operators/time_domain_problem_operator.cpp +++ b/src/problem_operators/time_domain_problem_operator.cpp @@ -20,6 +20,26 @@ GetTimeDerivativeNames(std::vector gridfunction_names) return time_derivative_names; } +TimeDomainProblemOperator::TimeDomainProblemOperator(hephaestus::Problem & problem) + : ProblemOperatorBase(problem) +{ + ConstructTimestepper(); +} + +void +TimeDomainProblemOperator::ConstructTimestepper() +{ + _ode_solver = std::make_unique(); + _ode_solver->Init(*this); +} + +void +TimeDomainProblemOperator::Init() +{ + ProblemOperatorBase::Init(); + _ode_solver->Init(*this); +} + void TimeDomainProblemOperator::Update() { @@ -27,10 +47,7 @@ TimeDomainProblemOperator::Update() // The dimensions of the problem operator have now changed. We must call the // ODE solver's Init method to resize its internal vector. - if (_ode_solver) - { - _ode_solver->Init(*this); - } + _ode_solver->Init(*this); } } // namespace hephaestus \ No newline at end of file diff --git a/src/problem_operators/time_domain_problem_operator.hpp b/src/problem_operators/time_domain_problem_operator.hpp index d36f17365..e12c5f27f 100644 --- a/src/problem_operators/time_domain_problem_operator.hpp +++ b/src/problem_operators/time_domain_problem_operator.hpp @@ -16,29 +16,28 @@ std::vector GetTimeDerivativeNames(std::vector gridfun class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public ProblemOperatorBase { public: - TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {} + TimeDomainProblemOperator(hephaestus::Problem & problem); ~TimeDomainProblemOperator() override = default; void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override {} - /// Set the ODE solver. - void SetODESolver(std::unique_ptr ode_solver) - { - _ode_solver = std::move(ode_solver); - } - /// Wrapper around the ODE solver's Step method using the block vector. virtual void Step(mfem::real_t & t, mfem::real_t & dt) { _ode_solver->Step(*_block_vector, t, dt); } + void Init() override; + void Update() override; protected: int & Width() final { return mfem::TimeDependentOperator::width; } int & Height() final { return mfem::TimeDependentOperator::height; } + /// Construct the ODE solver. Called in constructor. Override in derived classes. + virtual void ConstructTimestepper(); + private: /// Store the ODE solver. std::unique_ptr _ode_solver{nullptr}; From 18a158d3178ef68eab8294d3263759664566b396 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 14:26:42 +0000 Subject: [PATCH 05/30] Removes ConstructNonlinearSolver from problem builder classes; adds ConstructNonlinearSolver to problem operator classes; called in constructor. --- src/problem_builders/problem_builder_base.cpp | 16 ---------------- src/problem_builders/problem_builder_base.hpp | 2 -- src/problem_operators/problem_operator_base.cpp | 12 ++++++++++++ src/problem_operators/problem_operator_base.hpp | 4 ++++ 4 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index af4302aa6..ea18616e9 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -203,19 +203,6 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr_sources.Register(source_name, std::move(source)); } -void -ProblemBuilder::ConstructNonlinearSolver() -{ - auto nl_solver = std::make_unique(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()->GetOperator()->SetNonlinearSolver(std::move(nl_solver)); -} - void ProblemBuilder::InitializeSources() { @@ -261,9 +248,6 @@ ProblemBuilder::FinalizeProblem(bool build_operator) InitializeAuxSolvers(); InitializeSources(); InitializeOperator(); - - ConstructNonlinearSolver(); - InitializeOutputs(); } diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 3edec7ca3..ef0651d5f 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -91,8 +91,6 @@ class ProblemBuilder virtual void RegisterAuxSolvers() = 0; virtual void RegisterCoefficients() = 0; - virtual void ConstructNonlinearSolver(); - virtual void ConstructOperator() = 0; virtual void InitializeSources(); diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index bc7fc34b2..4bf568ef7 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -24,6 +24,17 @@ ProblemOperatorBase::ConstructJacobianSolver() _jacobian_solver = std::move(solver); } +void +ProblemOperatorBase::ConstructNonlinearSolver() +{ + _nonlinear_solver = std::make_unique(_problem._comm); + + // Defaults to one iteration, without further nonlinear iterations + _nonlinear_solver->SetRelTol(0.0); + _nonlinear_solver->SetAbsTol(0.0); + _nonlinear_solver->SetMaxIter(1); +} + void ProblemOperatorBase::SetTrialVariables() { @@ -99,6 +110,7 @@ ProblemOperatorBase::Init() { _block_vector = std::make_unique(); ConstructJacobianSolver(); + ConstructNonlinearSolver(); SetTrialVariables(); diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 92b55fe12..98af8afcd 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -42,6 +42,10 @@ class ProblemOperatorBase /// constructor. virtual void ConstructJacobianSolver(); + /// Override in derived classes to construct the non-linear solver. Called in + /// the constructor. + virtual void ConstructNonlinearSolver(); + /// Set trial variables names. Override in derived classes. virtual void SetTrialVariableNames() {} From b3884765d914d9e70ebdf2e96e0a5c34072fd58c Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 14:29:32 +0000 Subject: [PATCH 06/30] Removed SetNonlinearSolver setter. --- src/problem_operators/problem_operator_base.hpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 98af8afcd..53104b14b 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,12 +19,6 @@ class ProblemOperatorBase /// Update the problem operator after a mesh change. virtual void Update(); - /// Set the nonlinear solver. - void SetNonlinearSolver(std::unique_ptr nl_solver) - { - _nonlinear_solver = std::move(nl_solver); - } - /// Accessor for Jacobian preconditioner. template TSolver * JacobianPreconditioner() const From 9c8a27a5cfe2b75f293a03034cec34807f2b2828 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 11:48:33 +0000 Subject: [PATCH 07/30] Moves ConstructTimestepper into TimeDomainProblemOperator::Init This will allow the virtual method to be overriden in derived classes properly. Virtual methods will not work properly in constructors since derived classes will not yet exist. --- src/problem_operators/time_domain_problem_operator.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/problem_operators/time_domain_problem_operator.cpp b/src/problem_operators/time_domain_problem_operator.cpp index 2515b9e21..78ee2b096 100644 --- a/src/problem_operators/time_domain_problem_operator.cpp +++ b/src/problem_operators/time_domain_problem_operator.cpp @@ -23,7 +23,6 @@ GetTimeDerivativeNames(std::vector gridfunction_names) TimeDomainProblemOperator::TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) { - ConstructTimestepper(); } void @@ -37,7 +36,7 @@ void TimeDomainProblemOperator::Init() { ProblemOperatorBase::Init(); - _ode_solver->Init(*this); + ConstructTimestepper(); } void From 9ce6f232bef8f6bf0c69a85a498ab0b32e769933 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 11:53:11 +0000 Subject: [PATCH 08/30] Adds SolverOptions structure and SetSolverOptions method. --- .../problem_operator_base.cpp | 21 +++++++++++++------ .../problem_operator_base.hpp | 15 +++++++++++++ 2 files changed, 30 insertions(+), 6 deletions(-) diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index 4bf568ef7..4b497150d 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -12,16 +12,25 @@ ProblemOperatorBase::ConstructJacobianSolver() precond->SetPrintLevel(GetGlobalPrintLevel()); auto solver = std::make_unique(_problem._comm); - - solver->SetTol(1e-16); - solver->SetAbsTol(1e-16); - solver->SetMaxIter(1000); - solver->SetKDim(10); - solver->SetPrintLevel(GetGlobalPrintLevel()); solver->SetPreconditioner(*precond); _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); + + SolverOptions default_options; + SetSolverOptions(default_options); +} + +void +ProblemOperatorBase::SetSolverOptions(SolverOptions & options) +{ + auto & solver = static_cast(*_jacobian_solver); + + solver.SetTol(options._tolerance); + solver.SetAbsTol(options._abs_tolerance); + solver.SetMaxIter(options._max_iteration); + solver.SetKDim(options._k_dim); + solver.SetPrintLevel(options._print_level); } void diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 53104b14b..0822c2ef2 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -26,6 +26,21 @@ class ProblemOperatorBase return static_cast(_jacobian_preconditioner.get()); } + /// A structure for setting solver options. Defaults have been set. + struct SolverOptions + { + double _tolerance{1e-16}; + double _abs_tolerance{1e-16}; + + unsigned int _max_iteration{1000}; + + int _print_level{GetGlobalPrintLevel()}; + int _k_dim{10}; + }; + + /// Sets the solver's options. Override in derived classes. + virtual void SetSolverOptions(SolverOptions & options); + protected: /// Use of protected constructor to only allow construction by derived classes. /// All problem operator classes are built on-top of this class and it should not From 6d30e486023cb48d6740788118ebd427153a5744 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 13:27:41 +0000 Subject: [PATCH 09/30] Makes ConstructJacobianSolver protected in derived classes. --- .../ComplexMaxwell/complex_maxwell_formulation.hpp | 5 +++-- src/formulations/Dual/dual_formulation.hpp | 4 ++-- src/formulations/HCurl/hcurl_formulation.hpp | 2 ++ src/formulations/Magnetostatic/statics_formulation.hpp | 1 + 4 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp index 7f8dba130..eab385247 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp @@ -75,8 +75,6 @@ class ComplexMaxwellOperator : public ProblemOperator void Update() override; void Solve(mfem::Vector & X) override; - void ConstructJacobianSolver() override; - std::string _h_curl_var_complex_name, _h_curl_var_real_name, _h_curl_var_imag_name, _stiffness_coef_name, _mass_coef_name, _loss_coef_name; @@ -88,6 +86,9 @@ class ComplexMaxwellOperator : public ProblemOperator mfem::Coefficient * _loss_coef{nullptr}; // omega sigma mfem::Array _ess_bdr_tdofs; + +protected: + void ConstructJacobianSolver() override; }; } // namespace hephaestus diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index 2db78eaef..eae9e00bc 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -59,8 +59,6 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator void Update() override; - void ConstructJacobianSolver() override; - void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override; mfem::ParFiniteElementSpace * _h_curl_fe_space{nullptr}; @@ -72,6 +70,8 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator mfem::ParGridFunction * _dv{nullptr}; // HDiv vector field protected: + void ConstructJacobianSolver() override; + int GetSolutionVectorSize() const override; std::unique_ptr _curl; diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index c7a11f782..e039166b9 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -36,6 +36,8 @@ class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator : TimeDomainEquationSystemProblemOperator(problem, std::move(equation_system)) { } + +protected: void ConstructJacobianSolver() override; }; diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index d336430c1..a84b31386 100644 --- a/src/formulations/Magnetostatic/statics_formulation.hpp +++ b/src/formulations/Magnetostatic/statics_formulation.hpp @@ -38,6 +38,7 @@ class StaticsOperator : public ProblemOperator void Init() override; void Solve(mfem::Vector & X) override; +protected: void ConstructJacobianSolver() override; private: From 84df05c599336daf5558160251a8825df3d45108 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 13:34:39 +0000 Subject: [PATCH 10/30] Implements SetSolverOptions for HCurlProblemOperator. --- src/formulations/HCurl/hcurl_formulation.cpp | 19 ++++++++++++++----- src/formulations/HCurl/hcurl_formulation.hpp | 2 ++ 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 26eeca4d8..75eba1e71 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -88,6 +88,17 @@ HCurlFormulation::RegisterGridFunctions() TimeDomainEquationSystemProblemBuilder::RegisterGridFunctions(); }; +void +HCurlProblemOperator::SetSolverOptions(SolverOptions & options) +{ + auto & solver = static_cast(*_jacobian_solver); + + solver.SetTol(options._tolerance); + solver.SetAbsTol(options._abs_tolerance); + solver.SetMaxIter(options._max_iteration); + solver.SetPrintLevel(options._print_level); +} + void HCurlProblemOperator::ConstructJacobianSolver() { @@ -97,15 +108,13 @@ HCurlProblemOperator::ConstructJacobianSolver() precond->SetPrintLevel(-1); auto solver = std::make_unique(_problem._comm); - - solver->SetTol(1e-16); - solver->SetAbsTol(1e-16); - solver->SetMaxIter(1000); - solver->SetPrintLevel(GetGlobalPrintLevel()); solver->SetPreconditioner(*precond); _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); + + SolverOptions default_options; + SetSolverOptions(default_options); } CurlCurlEquationSystem::CurlCurlEquationSystem(const hephaestus::InputParameters & params) diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index e039166b9..9bbbd936d 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -37,6 +37,8 @@ class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator { } + void SetSolverOptions(SolverOptions & options) override; + protected: void ConstructJacobianSolver() override; }; From 798739e1fa351627a027525aadc537a6ff15cd74 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 13:38:09 +0000 Subject: [PATCH 11/30] Implements SetSolverOptions for DualOperator. --- src/formulations/Dual/dual_formulation.cpp | 19 ++++++++++++++----- src/formulations/Dual/dual_formulation.hpp | 2 ++ 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 7bea7a5d5..980a6609a 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -214,15 +214,24 @@ DualOperator::ConstructJacobianSolver() precond->SetPrintLevel(-1); auto solver = std::make_unique(_problem._comm); - - solver->SetTol(1e-16); - solver->SetAbsTol(1e-16); - solver->SetMaxIter(1000); - solver->SetPrintLevel(GetGlobalPrintLevel()); solver->SetPreconditioner(*precond); _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); + + SolverOptions default_options; + SetSolverOptions(default_options); +} + +void +DualOperator::SetSolverOptions(SolverOptions & options) +{ + auto & solver = static_cast(*_jacobian_solver); + + solver.SetTol(options._tolerance); + solver.SetAbsTol(options._abs_tolerance); + solver.SetMaxIter(options._max_iteration); + solver.SetPrintLevel(options._print_level); } void diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index eae9e00bc..39089e629 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -61,6 +61,8 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override; + void SetSolverOptions(SolverOptions & options) override; + mfem::ParFiniteElementSpace * _h_curl_fe_space{nullptr}; mfem::ParFiniteElementSpace * _h_div_fe_space{nullptr}; From 2383fa3545e372e4f5ecd7312fc2652614051a3f Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 13:41:20 +0000 Subject: [PATCH 12/30] Implements SetSolverOptions for ComplexMaxwellOperator. --- src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp index eab385247..de5c44658 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp @@ -75,6 +75,8 @@ class ComplexMaxwellOperator : public ProblemOperator void Update() override; void Solve(mfem::Vector & X) override; + void SetSolverOptions(SolverOptions & options) override {} + std::string _h_curl_var_complex_name, _h_curl_var_real_name, _h_curl_var_imag_name, _stiffness_coef_name, _mass_coef_name, _loss_coef_name; From 14eb5d7d775bc78315e7f3a252eaa116501bddba Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 13:43:48 +0000 Subject: [PATCH 13/30] Implements SetSolverOptions for StaticsOperator. --- .../Magnetostatic/statics_formulation.cpp | 19 ++++++++++++++----- .../Magnetostatic/statics_formulation.hpp | 2 ++ 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 818ac5ae7..add7d1298 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -108,15 +108,24 @@ StaticsOperator::ConstructJacobianSolver() precond->SetPrintLevel(-1); auto solver = std::make_unique(_problem._comm); - - solver->SetTol(1e-16); - solver->SetMaxIter(100); - solver->SetKDim(10); - solver->SetPrintLevel(GetGlobalPrintLevel()); solver->SetPreconditioner(*precond); _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); + + SolverOptions default_options = {._max_iteration = 100}; + SetSolverOptions(default_options); +} + +void +StaticsOperator::SetSolverOptions(SolverOptions & options) +{ + auto & solver = static_cast(*_jacobian_solver); + + solver.SetTol(options._tolerance); + solver.SetMaxIter(options._max_iteration); + solver.SetKDim(options._k_dim); + solver.SetPrintLevel(options._print_level); } /* diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index a84b31386..14334c959 100644 --- a/src/formulations/Magnetostatic/statics_formulation.hpp +++ b/src/formulations/Magnetostatic/statics_formulation.hpp @@ -38,6 +38,8 @@ class StaticsOperator : public ProblemOperator void Init() override; void Solve(mfem::Vector & X) override; + void SetSolverOptions(SolverOptions & options) override; + protected: void ConstructJacobianSolver() override; From 29d3137c16b0f51869eab44fc3d7b7393b407484 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 13:55:42 +0000 Subject: [PATCH 14/30] SetSolverOptions argument now passed by value. --- src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp | 2 +- src/formulations/Dual/dual_formulation.cpp | 2 +- src/formulations/Dual/dual_formulation.hpp | 2 +- src/formulations/HCurl/hcurl_formulation.cpp | 2 +- src/formulations/HCurl/hcurl_formulation.hpp | 2 +- src/formulations/Magnetostatic/statics_formulation.cpp | 2 +- src/formulations/Magnetostatic/statics_formulation.hpp | 2 +- src/problem_operators/problem_operator_base.cpp | 2 +- src/problem_operators/problem_operator_base.hpp | 2 +- 9 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp index de5c44658..07239a267 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp @@ -75,7 +75,7 @@ class ComplexMaxwellOperator : public ProblemOperator void Update() override; void Solve(mfem::Vector & X) override; - void SetSolverOptions(SolverOptions & options) override {} + void SetSolverOptions(SolverOptions options) override {} std::string _h_curl_var_complex_name, _h_curl_var_real_name, _h_curl_var_imag_name, _stiffness_coef_name, _mass_coef_name, _loss_coef_name; diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 980a6609a..f30822cd4 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -224,7 +224,7 @@ DualOperator::ConstructJacobianSolver() } void -DualOperator::SetSolverOptions(SolverOptions & options) +DualOperator::SetSolverOptions(SolverOptions options) { auto & solver = static_cast(*_jacobian_solver); diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index 39089e629..bf79c6cae 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -61,7 +61,7 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override; - void SetSolverOptions(SolverOptions & options) override; + void SetSolverOptions(SolverOptions options) override; mfem::ParFiniteElementSpace * _h_curl_fe_space{nullptr}; mfem::ParFiniteElementSpace * _h_div_fe_space{nullptr}; diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 75eba1e71..b31c63718 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -89,7 +89,7 @@ HCurlFormulation::RegisterGridFunctions() }; void -HCurlProblemOperator::SetSolverOptions(SolverOptions & options) +HCurlProblemOperator::SetSolverOptions(SolverOptions options) { auto & solver = static_cast(*_jacobian_solver); diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index 9bbbd936d..20533acf3 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -37,7 +37,7 @@ class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator { } - void SetSolverOptions(SolverOptions & options) override; + void SetSolverOptions(SolverOptions options) override; protected: void ConstructJacobianSolver() override; diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index add7d1298..cf5adc569 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -118,7 +118,7 @@ StaticsOperator::ConstructJacobianSolver() } void -StaticsOperator::SetSolverOptions(SolverOptions & options) +StaticsOperator::SetSolverOptions(SolverOptions options) { auto & solver = static_cast(*_jacobian_solver); diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index 14334c959..8e6e5030d 100644 --- a/src/formulations/Magnetostatic/statics_formulation.hpp +++ b/src/formulations/Magnetostatic/statics_formulation.hpp @@ -38,7 +38,7 @@ class StaticsOperator : public ProblemOperator void Init() override; void Solve(mfem::Vector & X) override; - void SetSolverOptions(SolverOptions & options) override; + void SetSolverOptions(SolverOptions options) override; protected: void ConstructJacobianSolver() override; diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index 4b497150d..bd4879e31 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -22,7 +22,7 @@ ProblemOperatorBase::ConstructJacobianSolver() } void -ProblemOperatorBase::SetSolverOptions(SolverOptions & options) +ProblemOperatorBase::SetSolverOptions(SolverOptions options) { auto & solver = static_cast(*_jacobian_solver); diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 0822c2ef2..f67e6ea81 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -39,7 +39,7 @@ class ProblemOperatorBase }; /// Sets the solver's options. Override in derived classes. - virtual void SetSolverOptions(SolverOptions & options); + virtual void SetSolverOptions(SolverOptions options); protected: /// Use of protected constructor to only allow construction by derived classes. From cc09a407ffb9d1df1f22d6a65bf8496e556ce221 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 15:16:24 +0000 Subject: [PATCH 15/30] Updates examples to use new SetSolverOptions. Replaces clunky InputParameters setup. --- examples/closed_coil/Main.cpp | 10 ++++------ examples/complex_team7/Main.cpp | 7 ++----- examples/magnetostatic/Main.cpp | 8 +++----- examples/open_coil/Main.cpp | 10 ++++------ examples/team4/Main.cpp | 10 ++++------ examples/team7/Main.cpp | 8 +++----- 6 files changed, 20 insertions(+), 33 deletions(-) diff --git a/examples/closed_coil/Main.cpp b/examples/closed_coil/Main.cpp index d7118c0c1..11cf01ebd 100644 --- a/examples/closed_coil/Main.cpp +++ b/examples/closed_coil/Main.cpp @@ -132,15 +132,13 @@ main(int argc, char * argv[]) hephaestus::Outputs outputs = defineOutputs(); problem_builder->SetOutputs(outputs); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-13)); - solver_options.SetParam("AbsTolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)500); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions( + {._tolerance = 1.0e-13, ._abs_tolerance = 1.0e-16, ._max_iteration = 500}); + hephaestus::InputParameters exec_params; exec_params.SetParam("VisualisationSteps", int(1)); exec_params.SetParam("UseGLVis", true); diff --git a/examples/complex_team7/Main.cpp b/examples/complex_team7/Main.cpp index 5adf63c35..29402d00f 100644 --- a/examples/complex_team7/Main.cpp +++ b/examples/complex_team7/Main.cpp @@ -225,15 +225,12 @@ main(int argc, char * argv[]) linesamplerwriter->SetPriority(5); problem_builder->AddPostprocessor("LineSamplerWriter", linesamplerwriter); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("Problem", static_cast(problem.get())); diff --git a/examples/magnetostatic/Main.cpp b/examples/magnetostatic/Main.cpp index 7cef1fd23..53441e181 100644 --- a/examples/magnetostatic/Main.cpp +++ b/examples/magnetostatic/Main.cpp @@ -160,14 +160,12 @@ main(int argc, char * argv[]) hephaestus::Outputs outputs = defineOutputs(); problem_builder->SetOutputs(outputs); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("VisualisationSteps", int(1)); exec_params.SetParam("UseGLVis", true); diff --git a/examples/open_coil/Main.cpp b/examples/open_coil/Main.cpp index 1cc1966dc..85654d4d5 100644 --- a/examples/open_coil/Main.cpp +++ b/examples/open_coil/Main.cpp @@ -159,15 +159,13 @@ main(int argc, char * argv[]) hephaestus::Outputs outputs = defineOutputs(); problem_builder->SetOutputs(outputs); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-10)); - solver_options.SetParam("AbsTolerance", float(1.0e-10)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions( + {._tolerance = 1.0e-10, ._abs_tolerance = 1.0e-10, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("VisualisationSteps", int(1)); exec_params.SetParam("UseGLVis", true); diff --git a/examples/team4/Main.cpp b/examples/team4/Main.cpp index 7d6e35de9..5cc3259d6 100644 --- a/examples/team4/Main.cpp +++ b/examples/team4/Main.cpp @@ -153,15 +153,13 @@ main(int argc, char * argv[]) fluxmonitor->SetPriority(2); problem_builder->AddPostprocessor("FluxMonitor", fluxmonitor); - hephaestus::InputParameters solver_options; - solver_options.SetParam("AbsTolerance", float(1.0e-20)); - solver_options.SetParam("Tolerance", float(1.0e-20)); - solver_options.SetParam("MaxIter", (unsigned int)500); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions( + {._tolerance = 1.0e-20, ._abs_tolerance = 1.0e-20, ._max_iteration = 500}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.001)); exec_params.SetParam("StartTime", float(0.00)); diff --git a/examples/team7/Main.cpp b/examples/team7/Main.cpp index 4fa8d52f4..a19125d64 100644 --- a/examples/team7/Main.cpp +++ b/examples/team7/Main.cpp @@ -214,14 +214,12 @@ main(int argc, char * argv[]) linesamplerwriter->SetPriority(5); problem_builder->AddPostprocessor("LineSamplerWriter", linesamplerwriter); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.001)); exec_params.SetParam("StartTime", float(0.00)); From 7482759da563533da4b09a5055fa72937e9f3673 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 15:22:22 +0000 Subject: [PATCH 16/30] Updates integration tests to use SetSolverOptions. --- test/integration/test_aform_source.cpp | 10 ++-------- test/integration/test_avform_rod.cpp | 3 --- test/integration/test_avform_source.cpp | 3 --- test/integration/test_complex_aform_rod.cpp | 10 ++-------- test/integration/test_complex_ermes_mouse.cpp | 10 ++-------- test/integration/test_complex_iris_wg.cpp | 10 ++-------- test/integration/test_complex_team7.cpp | 10 ++-------- test/integration/test_ebform_coupled.cpp | 11 +++-------- test/integration/test_ebform_rod.cpp | 11 +++-------- test/integration/test_hform_rod.cpp | 11 +++-------- test/integration/test_hform_source.cpp | 10 ++-------- test/integration/test_mesh_updates.cpp | 10 ++-------- test/integration/test_team4_hform.cpp | 10 ++++------ 13 files changed, 27 insertions(+), 92 deletions(-) diff --git a/test/integration/test_aform_source.cpp b/test/integration/test_aform_source.cpp index b60093299..def039144 100644 --- a/test/integration/test_aform_source.cpp +++ b/test/integration/test_aform_source.cpp @@ -110,10 +110,6 @@ class TestAFormSource current_solver_options, false)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("Mesh", mfem::ParMesh(MPI_COMM_WORLD, mesh)); params.SetParam("BoundaryConditions", bc_map); @@ -121,7 +117,6 @@ class TestAFormSource params.SetParam("PostProcessors", postprocessors); params.SetParam("Outputs", outputs); params.SetParam("Sources", sources); - params.SetParam("SolverOptions", solver_options); return params; } @@ -156,8 +151,6 @@ TEST_CASE_METHOD(TestAFormSource, "TestAForm", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P2")); @@ -169,12 +162,13 @@ TEST_CASE_METHOD(TestAFormSource, "TestAForm", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.05)); exec_params.SetParam("StartTime", float(0.00)); diff --git a/test/integration/test_avform_rod.cpp b/test/integration/test_avform_rod.cpp index 84c09cb20..926c90bd6 100644 --- a/test/integration/test_avform_rod.cpp +++ b/test/integration/test_avform_rod.cpp @@ -118,8 +118,6 @@ TEST_CASE_METHOD(TestAVFormRod, "TestAVFormRod", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->SetBoundaryConditions(bc_map); @@ -128,7 +126,6 @@ TEST_CASE_METHOD(TestAVFormRod, "TestAVFormRod", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); diff --git a/test/integration/test_avform_source.cpp b/test/integration/test_avform_source.cpp index 2b2026812..2aedc08c2 100644 --- a/test/integration/test_avform_source.cpp +++ b/test/integration/test_avform_source.cpp @@ -171,8 +171,6 @@ TEST_CASE_METHOD(TestAVFormSource, "TestAVFormSource", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P2")); @@ -185,7 +183,6 @@ TEST_CASE_METHOD(TestAVFormSource, "TestAVFormSource", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); diff --git a/test/integration/test_complex_aform_rod.cpp b/test/integration/test_complex_aform_rod.cpp index 12280df4f..6d43ea30a 100644 --- a/test/integration/test_complex_aform_rod.cpp +++ b/test/integration/test_complex_aform_rod.cpp @@ -108,10 +108,6 @@ class TestComplexAFormRod -1, current_solver_options)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-9)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("UseGLVis", true); params.SetParam("Mesh", mfem::ParMesh(MPI_COMM_WORLD, mesh)); @@ -122,7 +118,6 @@ class TestComplexAFormRod params.SetParam("PostProcessors", postprocessors); params.SetParam("Sources", sources); params.SetParam("Outputs", outputs); - params.SetParam("SolverOptions", solver_options); return params; } @@ -149,8 +144,6 @@ TEST_CASE_METHOD(TestComplexAFormRod, "TestComplexAFormRod", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P1")); @@ -187,12 +180,13 @@ TEST_CASE_METHOD(TestComplexAFormRod, "TestComplexAFormRod", "[CheckRun]") problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-9, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("Problem", static_cast(problem.get())); diff --git a/test/integration/test_complex_ermes_mouse.cpp b/test/integration/test_complex_ermes_mouse.cpp index ddac16f11..ffdafb2a5 100644 --- a/test/integration/test_complex_ermes_mouse.cpp +++ b/test/integration/test_complex_ermes_mouse.cpp @@ -101,10 +101,6 @@ class TestComplexERMESMouse hephaestus::AuxSolvers preprocessors; hephaestus::Sources sources; - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("UseGLVis", true); @@ -117,7 +113,6 @@ class TestComplexERMESMouse params.SetParam("PostProcessors", postprocessors); params.SetParam("Outputs", outputs); params.SetParam("Sources", sources); - params.SetParam("SolverOptions", solver_options); return params; } @@ -144,8 +139,6 @@ TEST_CASE_METHOD(TestComplexERMESMouse, "TestComplexERMESMouse", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->SetBoundaryConditions(bc_map); @@ -154,12 +147,13 @@ TEST_CASE_METHOD(TestComplexERMESMouse, "TestComplexERMESMouse", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("Problem", static_cast(problem.get())); diff --git a/test/integration/test_complex_iris_wg.cpp b/test/integration/test_complex_iris_wg.cpp index e2611cd7c..7b06130ba 100644 --- a/test/integration/test_complex_iris_wg.cpp +++ b/test/integration/test_complex_iris_wg.cpp @@ -99,10 +99,6 @@ class TestComplexIrisWaveguide hephaestus::AuxSolvers preprocessors; hephaestus::Sources sources; - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("UseGLVis", true); @@ -115,7 +111,6 @@ class TestComplexIrisWaveguide params.SetParam("PostProcessors", postprocessors); params.SetParam("Outputs", outputs); params.SetParam("Sources", sources); - params.SetParam("SolverOptions", solver_options); return params; } @@ -142,8 +137,6 @@ TEST_CASE_METHOD(TestComplexIrisWaveguide, "TestComplexIrisWaveguide", "[CheckRu auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->SetBoundaryConditions(bc_map); @@ -172,12 +165,13 @@ TEST_CASE_METHOD(TestComplexIrisWaveguide, "TestComplexIrisWaveguide", "[CheckRu problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("Problem", static_cast(problem.get())); diff --git a/test/integration/test_complex_team7.cpp b/test/integration/test_complex_team7.cpp index 923390be9..4d57f6f08 100644 --- a/test/integration/test_complex_team7.cpp +++ b/test/integration/test_complex_team7.cpp @@ -135,10 +135,6 @@ class TestComplexTeam7 std::make_shared( "source", "source", "HCurl", "H1", "electric_potential", current_solver_options)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("UseGLVis", true); @@ -147,7 +143,6 @@ class TestComplexTeam7 params.SetParam("Coefficients", coefficients); params.SetParam("Outputs", outputs); params.SetParam("Sources", sources); - params.SetParam("SolverOptions", solver_options); hephaestus::logger.info("Created params "); return params; } @@ -166,8 +161,6 @@ TEST_CASE_METHOD(TestComplexTeam7, "TestComplexTeam7", "[CheckRun]") auto coefficients(params.GetParam("Coefficients")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); auto problem_builder = std::make_unique("magnetic_reluctivity", @@ -192,7 +185,6 @@ TEST_CASE_METHOD(TestComplexTeam7, "TestComplexTeam7", "[CheckRun]") problem_builder->SetCoefficients(coefficients); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); // Call LineSampler to save values std::string gridfunction_name("magnetic_flux_density_real"); @@ -220,6 +212,8 @@ TEST_CASE_METHOD(TestComplexTeam7, "TestComplexTeam7", "[CheckRun]") auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("Problem", static_cast(problem.get())); diff --git a/test/integration/test_ebform_coupled.cpp b/test/integration/test_ebform_coupled.cpp index a1ee79e86..9178a28c7 100644 --- a/test/integration/test_ebform_coupled.cpp +++ b/test/integration/test_ebform_coupled.cpp @@ -125,10 +125,6 @@ class TestEBFormCoupled -1, current_solver_options)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-9)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("Mesh", mfem::ParMesh(MPI_COMM_WORLD, mesh)); params.SetParam("BoundaryConditions", bc_map); @@ -139,7 +135,6 @@ class TestEBFormCoupled params.SetParam("PostProcessors", postprocessors); params.SetParam("Sources", sources); params.SetParam("Outputs", outputs); - params.SetParam("SolverOptions", solver_options); return params; } @@ -161,8 +156,6 @@ TEST_CASE_METHOD(TestEBFormCoupled, "TestEBFormCoupled", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); std::shared_ptr pmesh = std::make_shared(params.GetParam("Mesh")); @@ -185,11 +178,13 @@ TEST_CASE_METHOD(TestEBFormCoupled, "TestEBFormCoupled", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-9, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.5)); exec_params.SetParam("StartTime", float(0.00)); diff --git a/test/integration/test_ebform_rod.cpp b/test/integration/test_ebform_rod.cpp index d9a7bed54..f8d026c07 100644 --- a/test/integration/test_ebform_rod.cpp +++ b/test/integration/test_ebform_rod.cpp @@ -103,10 +103,6 @@ class TestEBFormRod -1, current_solver_options)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-9)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("Mesh", mfem::ParMesh(MPI_COMM_WORLD, mesh)); params.SetParam("BoundaryConditions", bc_map); @@ -116,7 +112,6 @@ class TestEBFormRod params.SetParam("PostProcessors", postprocessors); params.SetParam("Sources", sources); params.SetParam("Outputs", outputs); - params.SetParam("SolverOptions", solver_options); return params; } @@ -138,8 +133,6 @@ TEST_CASE_METHOD(TestEBFormRod, "TestEBFormRod", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); std::shared_ptr pmesh = std::make_shared(params.GetParam("Mesh")); @@ -153,11 +146,13 @@ TEST_CASE_METHOD(TestEBFormRod, "TestEBFormRod", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-9, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.5)); exec_params.SetParam("StartTime", float(0.00)); diff --git a/test/integration/test_hform_rod.cpp b/test/integration/test_hform_rod.cpp index 74e163773..8444df3e9 100644 --- a/test/integration/test_hform_rod.cpp +++ b/test/integration/test_hform_rod.cpp @@ -88,10 +88,6 @@ class TestHFormRod outputs.Register("ParaViewDataCollection", std::make_shared("HFormParaView")); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::GridFunctions gridfunctions; hephaestus::AuxSolvers preprocessors; hephaestus::AuxSolvers postprocessors; @@ -118,7 +114,6 @@ class TestHFormRod params.SetParam("PostProcessors", postprocessors); params.SetParam("Outputs", outputs); params.SetParam("Sources", sources); - params.SetParam("SolverOptions", solver_options); return params; } @@ -141,8 +136,6 @@ TEST_CASE_METHOD(TestHFormRod, "TestHFormRod", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->AddFESpace(std::string("H1"), std::string("H1_3D_P2")); @@ -152,11 +145,13 @@ TEST_CASE_METHOD(TestHFormRod, "TestHFormRod", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.5)); exec_params.SetParam("StartTime", float(0.00)); diff --git a/test/integration/test_hform_source.cpp b/test/integration/test_hform_source.cpp index dbf26f531..604036bc9 100644 --- a/test/integration/test_hform_source.cpp +++ b/test/integration/test_hform_source.cpp @@ -116,10 +116,6 @@ class TestHFormSource current_solver_options, false)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("TimeStep", float(0.05)); params.SetParam("StartTime", float(0.00)); @@ -134,7 +130,6 @@ class TestHFormSource params.SetParam("PostProcessors", postprocessors); params.SetParam("Sources", sources); params.SetParam("Outputs", outputs); - params.SetParam("SolverOptions", solver_options); return params; } @@ -171,8 +166,6 @@ TEST_CASE_METHOD(TestHFormSource, "TestHFormSource", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P2")); @@ -184,12 +177,13 @@ TEST_CASE_METHOD(TestHFormSource, "TestHFormSource", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.05)); exec_params.SetParam("StartTime", float(0.00)); diff --git a/test/integration/test_mesh_updates.cpp b/test/integration/test_mesh_updates.cpp index 277488a19..d15982c98 100644 --- a/test/integration/test_mesh_updates.cpp +++ b/test/integration/test_mesh_updates.cpp @@ -112,10 +112,6 @@ class TestMeshUpdates current_solver_options, false)); - hephaestus::InputParameters solver_options; - solver_options.SetParam("Tolerance", float(1.0e-16)); - solver_options.SetParam("MaxIter", (unsigned int)1000); - hephaestus::InputParameters params; params.SetParam("Mesh", mfem::ParMesh(MPI_COMM_WORLD, mesh)); params.SetParam("BoundaryConditions", bc_map); @@ -123,7 +119,6 @@ class TestMeshUpdates params.SetParam("PostProcessors", postprocessors); params.SetParam("Outputs", outputs); params.SetParam("Sources", sources); - params.SetParam("SolverOptions", solver_options); return params; } @@ -146,8 +141,6 @@ TEST_CASE_METHOD(TestMeshUpdates, "TestMeshUpdates", "[CheckRun]") auto postprocessors(params.GetParam("PostProcessors")); auto sources(params.GetParam("Sources")); auto outputs(params.GetParam("Outputs")); - auto solver_options(params.GetOptionalParam( - "SolverOptions", hephaestus::InputParameters())); problem_builder->SetMesh(pmesh); problem_builder->AddFESpace(std::string("HCurl"), std::string("ND_3D_P2")); @@ -158,12 +151,13 @@ TEST_CASE_METHOD(TestMeshUpdates, "TestMeshUpdates", "[CheckRun]") problem_builder->SetPostprocessors(postprocessors); problem_builder->SetSources(sources); problem_builder->SetOutputs(outputs); - problem_builder->SetSolverOptions(solver_options); problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + problem->GetOperator()->SetSolverOptions({._tolerance = 1.0e-16, ._max_iteration = 1000}); + hephaestus::InputParameters exec_params; // Single timestep. diff --git a/test/integration/test_team4_hform.cpp b/test/integration/test_team4_hform.cpp index 8c226d203..63813bc30 100644 --- a/test/integration/test_team4_hform.cpp +++ b/test/integration/test_team4_hform.cpp @@ -147,15 +147,13 @@ TEST_CASE_METHOD(TestTEAM4HForm, "TestTEAM4HForm", "[CheckRun]") fluxmonitor->SetPriority(2); problem_builder->AddPostprocessor("FluxMonitor", fluxmonitor); - hephaestus::InputParameters solver_options; - solver_options.SetParam("AbsTolerance", float(1.0e-20)); - solver_options.SetParam("Tolerance", float(1.0e-20)); - solver_options.SetParam("MaxIter", (unsigned int)500); - problem_builder->SetSolverOptions(solver_options); - problem_builder->FinalizeProblem(); auto problem = problem_builder->ReturnProblem(); + + problem->GetOperator()->SetSolverOptions( + {._tolerance = 1.0e-20, ._abs_tolerance = 1.0e-20, ._max_iteration = 500}); + hephaestus::InputParameters exec_params; exec_params.SetParam("TimeStep", float(0.001)); exec_params.SetParam("StartTime", float(0.00)); From b9f3ae1abb8977cf5ca8c194900bcfb61832c517 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 15:30:26 +0000 Subject: [PATCH 17/30] Removes solver input parameters from Problem and SetSolverOptions from ProblemBuilder. --- src/problem_builders/problem_builder_base.cpp | 7 ------- src/problem_builders/problem_builder_base.hpp | 2 -- 2 files changed, 9 deletions(-) diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index ea18616e9..a125a5e24 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -82,13 +82,6 @@ ProblemBuilder::SetOutputs(hephaestus::Outputs & outputs) GetProblem()->_outputs = outputs; } -void -ProblemBuilder::SetSolverOptions(hephaestus::InputParameters & solver_options) -{ - logger.info("Setting Solver Options"); - GetProblem()->_solver_options = solver_options; -} - void ProblemBuilder::SetCoefficients(hephaestus::Coefficients & coefficients) { diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index ef0651d5f..79423f8ce 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -28,7 +28,6 @@ class Problem hephaestus::AuxSolvers _postprocessors; hephaestus::Sources _sources; hephaestus::Outputs _outputs; - hephaestus::InputParameters _solver_options; hephaestus::FECollections _fecs; hephaestus::FESpaces _fespaces; @@ -72,7 +71,6 @@ class ProblemBuilder void SetPostprocessors(hephaestus::AuxSolvers & postprocessors); void SetSources(hephaestus::Sources & sources); void SetOutputs(hephaestus::Outputs & outputs); - void SetSolverOptions(hephaestus::InputParameters & solver_options); void SetCoefficients(hephaestus::Coefficients & coefficients); void AddFESpace(std::string fespace_name, From fe35b4244f510e68a1e37404a7c4484bb2b6ebc2 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 16:08:08 +0000 Subject: [PATCH 18/30] Adds custom Update implementation for HCurlProblemOperator --- src/formulations/HCurl/hcurl_formulation.cpp | 26 +++++++++++++++++ src/formulations/HCurl/hcurl_formulation.hpp | 3 ++ ...omain_equation_system_problem_operator.cpp | 28 +------------------ 3 files changed, 30 insertions(+), 27 deletions(-) diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index b31c63718..3dec5ea03 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -117,6 +117,32 @@ HCurlProblemOperator::ConstructJacobianSolver() SetSolverOptions(default_options); } +void +HCurlProblemOperator::Update() +{ + // 0. Call superclass' update method. + TimeDomainEquationSystemProblemOperator::Update(); + + // 1. Rebuild jacobian preconditioner. + auto precond = std::make_unique(GetEquationSystem()->_test_pfespaces.at(0)); + + precond->SetSingularProblem(); + precond->SetPrintLevel(-1); + + // 2. Set new preconditioner. + auto & solver = static_cast(*_jacobian_solver); + solver.SetPreconditioner(*precond); + + // 3. Rebuild jacobian matrix and set. + GetEquationSystem()->BuildJacobian(_true_x, _true_rhs); + + auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); + solver.SetOperator(*matrix); + + // 4. Set preconditioner. + _jacobian_preconditioner = std::move(precond); +} + CurlCurlEquationSystem::CurlCurlEquationSystem(const hephaestus::InputParameters & params) : _h_curl_var_name(params.GetParam("HCurlVarName")), _alpha_coef_name(params.GetParam("AlphaCoefName")), diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index 20533acf3..c216644e2 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -39,6 +39,9 @@ class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator void SetSolverOptions(SolverOptions options) override; + /// Provide a custom update method for chosen solvers. + void Update() override; + protected: void ConstructJacobianSolver() override; }; diff --git a/src/problem_operators/time_domain_equation_system_problem_operator.cpp b/src/problem_operators/time_domain_equation_system_problem_operator.cpp index ccc210c42..d39872839 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -43,33 +43,7 @@ TimeDomainEquationSystemProblemOperator::Update() TimeDomainProblemOperator::Update(); - // TODO: - we need to update the size of the jacobian_solver here after the parent class' Update - // method is called which ensures that we've updated the _true_x, _true_rhs. This is a hacky first - // attempt to avoid a segfault with test_mesh_updates. This is very bad code and basically - // copies methods from the AFormulation problem builder used. We need to generalize this. - // Note that this approach will only work for this specific problem! Static pointer casts - // have been used to ensure that this will fail if used on a different problem. - { - GetEquationSystem()->BuildJacobian(_true_x, _true_rhs); - - // Rebuild the jacobian preconditioner. - auto first_pfespace = GetEquationSystem()->_test_pfespaces.at(0); - - auto precond = std::make_unique(first_pfespace); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - _jacobian_preconditioner = std::move(precond); - - // Set new preconditioner. - static_cast(_jacobian_solver.get()) - ->SetPreconditioner(*static_cast(_jacobian_preconditioner.get())); - - // Set Jacobian matrix. - auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); - _jacobian_solver->SetOperator(*matrix); - } + // NB: derived classes should update Jacobian solver and preconditioner here. } void From 1b582547e26ab3c1eab8d51b166b05b3fbd7a572 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 16:18:25 +0000 Subject: [PATCH 19/30] Adds DualOperator update implementation. --- src/formulations/Dual/dual_formulation.cpp | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index f30822cd4..394ef3c9d 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -199,8 +199,29 @@ DualOperator::Init() void DualOperator::Update() { + // 0. Call superclass' update method. TimeDomainEquationSystemProblemOperator::Update(); + // 1. Rebuild jacobian preconditioner. + auto precond = std::make_unique(GetEquationSystem()->_test_pfespaces.at(0)); + + precond->SetSingularProblem(); + precond->SetPrintLevel(-1); + + // 2. Set new preconditioner. + auto & solver = static_cast(*_jacobian_solver); + solver.SetPreconditioner(*precond); + + // 3. Rebuild jacobian matrix and set. + GetEquationSystem()->BuildJacobian(_true_x, _true_rhs); + + auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); + solver.SetOperator(*matrix); + + // 4. Set preconditioner. + _jacobian_preconditioner = std::move(precond); + + // 5. Update curl. _curl->Update(); _curl->Assemble(); } From c801b5642dac507deb5a5d42f5cf83443ce064e8 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 16:22:00 +0000 Subject: [PATCH 20/30] Corrected comments in ProblemOperatorBase. --- src/problem_operators/problem_operator_base.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index f67e6ea81..a4996a83a 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -48,11 +48,11 @@ class ProblemOperatorBase explicit ProblemOperatorBase(hephaestus::Problem & problem); /// Override in derived classes to construct the Jacobian solver. Called in the - /// constructor. + /// Init method. virtual void ConstructJacobianSolver(); /// Override in derived classes to construct the non-linear solver. Called in - /// the constructor. + /// the Init method. virtual void ConstructNonlinearSolver(); /// Set trial variables names. Override in derived classes. From 73f9c965ab85cf6f68647d680b904d3dfd1a305c Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 17:41:31 +0000 Subject: [PATCH 21/30] Switched order in ComplexMaxwellOperator::Update to be consistent --- src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp index 3b3bc56bd..6d60953e5 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp @@ -162,8 +162,8 @@ ComplexMaxwellOperator::Init() void ComplexMaxwellOperator::Update() { - _u->Update(); ProblemOperator::Update(); + _u->Update(); } void From f574312958d2c6a393ac3a5c144d084f20bf24eb Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 17:47:03 +0000 Subject: [PATCH 22/30] Switches order in TimeDomainEquationSystemProblemOperator::Update to be consistent. --- .../time_domain_equation_system_problem_operator.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/problem_operators/time_domain_equation_system_problem_operator.cpp b/src/problem_operators/time_domain_equation_system_problem_operator.cpp index d39872839..86588af50 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -39,11 +39,9 @@ TimeDomainEquationSystemProblemOperator::Init() void TimeDomainEquationSystemProblemOperator::Update() { - GetEquationSystem()->Update(_problem._bc_map, _problem._sources); - TimeDomainProblemOperator::Update(); - // NB: derived classes should update Jacobian solver and preconditioner here. + GetEquationSystem()->Update(_problem._bc_map, _problem._sources); } void From 673d1ebc462b6036c30287b07f6b0ce5039c8f9b Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 17:48:45 +0000 Subject: [PATCH 23/30] Adds solver options member variable to store parameters; now call ConstructJacobianSolver in base class for updates. This is faster and much simpler than having to provide a custom update for each solver. --- src/formulations/Dual/dual_formulation.cpp | 25 ++------------- src/formulations/HCurl/hcurl_formulation.cpp | 31 ++----------------- src/formulations/HCurl/hcurl_formulation.hpp | 3 -- .../Magnetostatic/statics_formulation.cpp | 5 +-- .../problem_operator_base.cpp | 10 ++++-- .../problem_operator_base.hpp | 3 ++ 6 files changed, 20 insertions(+), 57 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 394ef3c9d..46f57341a 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -202,26 +202,6 @@ DualOperator::Update() // 0. Call superclass' update method. TimeDomainEquationSystemProblemOperator::Update(); - // 1. Rebuild jacobian preconditioner. - auto precond = std::make_unique(GetEquationSystem()->_test_pfespaces.at(0)); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - // 2. Set new preconditioner. - auto & solver = static_cast(*_jacobian_solver); - solver.SetPreconditioner(*precond); - - // 3. Rebuild jacobian matrix and set. - GetEquationSystem()->BuildJacobian(_true_x, _true_rhs); - - auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); - solver.SetOperator(*matrix); - - // 4. Set preconditioner. - _jacobian_preconditioner = std::move(precond); - - // 5. Update curl. _curl->Update(); _curl->Assemble(); } @@ -240,8 +220,7 @@ DualOperator::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - SolverOptions default_options; - SetSolverOptions(default_options); + SetSolverOptions(_solver_options); } void @@ -253,6 +232,8 @@ DualOperator::SetSolverOptions(SolverOptions options) solver.SetAbsTol(options._abs_tolerance); solver.SetMaxIter(options._max_iteration); solver.SetPrintLevel(options._print_level); + + _solver_options = options; } void diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 3dec5ea03..f08457bd7 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -97,6 +97,8 @@ HCurlProblemOperator::SetSolverOptions(SolverOptions options) solver.SetAbsTol(options._abs_tolerance); solver.SetMaxIter(options._max_iteration); solver.SetPrintLevel(options._print_level); + + _solver_options = options; } void @@ -113,34 +115,7 @@ HCurlProblemOperator::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - SolverOptions default_options; - SetSolverOptions(default_options); -} - -void -HCurlProblemOperator::Update() -{ - // 0. Call superclass' update method. - TimeDomainEquationSystemProblemOperator::Update(); - - // 1. Rebuild jacobian preconditioner. - auto precond = std::make_unique(GetEquationSystem()->_test_pfespaces.at(0)); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - // 2. Set new preconditioner. - auto & solver = static_cast(*_jacobian_solver); - solver.SetPreconditioner(*precond); - - // 3. Rebuild jacobian matrix and set. - GetEquationSystem()->BuildJacobian(_true_x, _true_rhs); - - auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); - solver.SetOperator(*matrix); - - // 4. Set preconditioner. - _jacobian_preconditioner = std::move(precond); + SetSolverOptions(_solver_options); } CurlCurlEquationSystem::CurlCurlEquationSystem(const hephaestus::InputParameters & params) diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index c216644e2..20533acf3 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -39,9 +39,6 @@ class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator void SetSolverOptions(SolverOptions options) override; - /// Provide a custom update method for chosen solvers. - void Update() override; - protected: void ConstructJacobianSolver() override; }; diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index cf5adc569..156243c3e 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -113,8 +113,7 @@ StaticsOperator::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - SolverOptions default_options = {._max_iteration = 100}; - SetSolverOptions(default_options); + SetSolverOptions(_solver_options); } void @@ -126,6 +125,8 @@ StaticsOperator::SetSolverOptions(SolverOptions options) solver.SetMaxIter(options._max_iteration); solver.SetKDim(options._k_dim); solver.SetPrintLevel(options._print_level); + + _solver_options = options; } /* diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index bd4879e31..bace1bd33 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -17,8 +17,8 @@ ProblemOperatorBase::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - SolverOptions default_options; - SetSolverOptions(default_options); + // Set default options or user-options if available. + SetSolverOptions(_solver_options); } void @@ -31,6 +31,9 @@ ProblemOperatorBase::SetSolverOptions(SolverOptions options) solver.SetMaxIter(options._max_iteration); solver.SetKDim(options._k_dim); solver.SetPrintLevel(options._print_level); + + // Store the options for future. + _solver_options = options; } void @@ -136,6 +139,9 @@ ProblemOperatorBase::Update() UpdateOffsets(); UpdateBlockVector(*_block_vector); + + // Update the Jacobian solver. + ConstructJacobianSolver(); } } \ No newline at end of file diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index a4996a83a..47b51e7fb 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -102,6 +102,9 @@ class ProblemOperatorBase /// Store the block vector. std::unique_ptr _block_vector{nullptr}; + /// The current solver options. + SolverOptions _solver_options; + private: /// Update the block vectors and offsets after a mesh change. void UpdateOffsets(); From 3d9071c20ac2ce4c1c3d12d44f2c5e067617724b Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 17:52:43 +0000 Subject: [PATCH 24/30] Removes jacobian preconditioner template accessor from ProblemOperatorBase. --- src/problem_operators/problem_operator_base.hpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 47b51e7fb..8ed005b1e 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,13 +19,6 @@ class ProblemOperatorBase /// Update the problem operator after a mesh change. virtual void Update(); - /// Accessor for Jacobian preconditioner. - template - TSolver * JacobianPreconditioner() const - { - return static_cast(_jacobian_preconditioner.get()); - } - /// A structure for setting solver options. Defaults have been set. struct SolverOptions { From 5046623534a36aa599c335bb3893b8785acae8a6 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 17:53:34 +0000 Subject: [PATCH 25/30] Removes Jacobian operator handle from EquationSystem. --- src/equation_systems/equation_system.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/equation_systems/equation_system.hpp b/src/equation_systems/equation_system.hpp index 613bbe706..d1278b242 100644 --- a/src/equation_systems/equation_system.hpp +++ b/src/equation_systems/equation_system.hpp @@ -97,9 +97,6 @@ class EquationSystem : public mfem::Operator virtual void RecoverFEMSolution(mfem::BlockVector & trueX, hephaestus::GridFunctions & gridfunctions); - /// Returns a reference to the jacobian operator handle. - [[nodiscard]] mfem::OperatorHandle & JacobianOperatorHandle() { return _jacobian; } - std::vector> _ess_tdof_lists; protected: From f57eb82c28c50169123ddd8cc042b058d4a05425 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Thu, 23 May 2024 18:08:31 +0000 Subject: [PATCH 26/30] Removes ConstructOperator from problem classes. Implemented directly in problem builders. --- src/problem_builders/problem_builder_base.hpp | 3 --- ...ady_state_equation_system_problem_builder.cpp | 16 ++++++++++++++++ ...ady_state_equation_system_problem_builder.hpp | 11 ++--------- .../steady_state_problem_builder.cpp | 3 ++- .../steady_state_problem_builder.hpp | 6 ------ ...me_domain_equation_system_problem_builder.cpp | 15 +++++++++++++++ ...me_domain_equation_system_problem_builder.hpp | 11 ++--------- .../time_domain_problem_builder.cpp | 3 ++- .../time_domain_problem_builder.hpp | 6 ------ 9 files changed, 39 insertions(+), 35 deletions(-) create mode 100644 src/problem_builders/steady_state_equation_system_problem_builder.cpp create mode 100644 src/problem_builders/time_domain_equation_system_problem_builder.cpp diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 79423f8ce..1778d84e6 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -40,9 +40,6 @@ class Problem /// Returns a pointer to the operator. See derived classes. [[nodiscard]] virtual hephaestus::ProblemOperatorBase * GetOperator() const = 0; - /// Virtual method to construct the operator. Call for default problems. - virtual void ConstructOperator() = 0; - /// Call to update on mesh change. virtual void Update(); }; diff --git a/src/problem_builders/steady_state_equation_system_problem_builder.cpp b/src/problem_builders/steady_state_equation_system_problem_builder.cpp new file mode 100644 index 000000000..0c0f94557 --- /dev/null +++ b/src/problem_builders/steady_state_equation_system_problem_builder.cpp @@ -0,0 +1,16 @@ +#include "steady_state_equation_system_problem_builder.hpp" + +namespace hephaestus +{ + +void +SteadyStateEquationSystemProblemBuilder::ConstructOperator() +{ + auto equation_system = std::make_unique(); + auto problem_operator = std::make_unique( + *GetProblem(), std::move(equation_system)); + + GetProblem()->SetOperator(std::move(problem_operator)); +} + +} // namespace hephaestus diff --git a/src/problem_builders/steady_state_equation_system_problem_builder.hpp b/src/problem_builders/steady_state_equation_system_problem_builder.hpp index 32c4c57cb..37d1484ea 100644 --- a/src/problem_builders/steady_state_equation_system_problem_builder.hpp +++ b/src/problem_builders/steady_state_equation_system_problem_builder.hpp @@ -24,15 +24,6 @@ class SteadyStateEquationSystemProblem : public SteadyStateProblem, public Equat SteadyStateProblem::SetOperator(std::move(problem_operator)); } - void ConstructOperator() override - { - auto equation_system = std::make_unique(); - auto problem_operator = std::make_unique( - *this, std::move(equation_system)); - - SetOperator(std::move(problem_operator)); - } - [[nodiscard]] hephaestus::EquationSystem * GetEquationSystem() const override { return GetOperator()->GetEquationSystem(); @@ -54,6 +45,8 @@ class SteadyStateEquationSystemProblemBuilder : public SteadyStateProblemBuilder auto ReturnProblem() { return ProblemBuilder::ReturnProblem(); } + void ConstructOperator() override; + protected: [[nodiscard]] hephaestus::SteadyStateEquationSystemProblem * GetProblem() const override { diff --git a/src/problem_builders/steady_state_problem_builder.cpp b/src/problem_builders/steady_state_problem_builder.cpp index a9fbe8ecd..60eca97bd 100644 --- a/src/problem_builders/steady_state_problem_builder.cpp +++ b/src/problem_builders/steady_state_problem_builder.cpp @@ -6,7 +6,8 @@ namespace hephaestus void SteadyStateProblemBuilder::ConstructOperator() { - GetProblem()->ConstructOperator(); + auto problem_operator = std::make_unique(*GetProblem()); + GetProblem()->SetOperator(std::move(problem_operator)); } } // namespace hephaestus diff --git a/src/problem_builders/steady_state_problem_builder.hpp b/src/problem_builders/steady_state_problem_builder.hpp index 731caf314..9765fe975 100644 --- a/src/problem_builders/steady_state_problem_builder.hpp +++ b/src/problem_builders/steady_state_problem_builder.hpp @@ -27,12 +27,6 @@ class SteadyStateProblem : public Problem _problem_operator = std::move(problem_operator); } - void ConstructOperator() override - { - _problem_operator.reset(); - _problem_operator = std::make_unique(*this); - } - private: std::unique_ptr _problem_operator{nullptr}; }; diff --git a/src/problem_builders/time_domain_equation_system_problem_builder.cpp b/src/problem_builders/time_domain_equation_system_problem_builder.cpp new file mode 100644 index 000000000..b26680b0c --- /dev/null +++ b/src/problem_builders/time_domain_equation_system_problem_builder.cpp @@ -0,0 +1,15 @@ +#include "time_domain_equation_system_problem_builder.hpp" + +namespace hephaestus +{ +void +TimeDomainEquationSystemProblemBuilder::ConstructOperator() +{ + auto equation_system = std::make_unique(); + auto problem_operator = std::make_unique( + *GetProblem(), std::move(equation_system)); + + GetProblem()->SetOperator(std::move(problem_operator)); +} + +} // namespace hephaestus diff --git a/src/problem_builders/time_domain_equation_system_problem_builder.hpp b/src/problem_builders/time_domain_equation_system_problem_builder.hpp index a96620911..dc436199a 100644 --- a/src/problem_builders/time_domain_equation_system_problem_builder.hpp +++ b/src/problem_builders/time_domain_equation_system_problem_builder.hpp @@ -25,15 +25,6 @@ class TimeDomainEquationSystemProblem : public TimeDomainProblem, public Equatio TimeDomainProblem::SetOperator(std::move(problem_operator)); } - void ConstructOperator() override - { - auto equation_system = std::make_unique(); - auto problem_operator = std::make_unique( - *this, std::move(equation_system)); - - SetOperator(std::move(problem_operator)); - } - [[nodiscard]] TimeDependentEquationSystem * GetEquationSystem() const override { return GetOperator()->GetEquationSystem(); @@ -55,6 +46,8 @@ class TimeDomainEquationSystemProblemBuilder : public TimeDomainProblemBuilder, auto ReturnProblem() { return ProblemBuilder::ReturnProblem(); } + void ConstructOperator() override; + protected: [[nodiscard]] hephaestus::TimeDomainEquationSystemProblem * GetProblem() const override { diff --git a/src/problem_builders/time_domain_problem_builder.cpp b/src/problem_builders/time_domain_problem_builder.cpp index 93fee3693..8ebe4e105 100644 --- a/src/problem_builders/time_domain_problem_builder.cpp +++ b/src/problem_builders/time_domain_problem_builder.cpp @@ -35,7 +35,8 @@ TimeDomainProblemBuilder::RegisterGridFunctions() void TimeDomainProblemBuilder::ConstructOperator() { - GetProblem()->ConstructOperator(); + auto problem_operator = std::make_unique(*GetProblem()); + GetProblem()->SetOperator(std::move(problem_operator)); } void diff --git a/src/problem_builders/time_domain_problem_builder.hpp b/src/problem_builders/time_domain_problem_builder.hpp index 0134c339f..d21cd422c 100644 --- a/src/problem_builders/time_domain_problem_builder.hpp +++ b/src/problem_builders/time_domain_problem_builder.hpp @@ -25,12 +25,6 @@ class TimeDomainProblem : public Problem _problem_operator = std::move(problem_operator); } - void ConstructOperator() override - { - _problem_operator.reset(); - _problem_operator = std::make_unique(*this); - } - private: std::unique_ptr _problem_operator{nullptr}; }; From e4c038e46fccca1b0c14292b469c82cc26858c8d Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Fri, 24 May 2024 10:04:34 +0000 Subject: [PATCH 27/30] Removed TODO comment from EquationSystemProblemOperator. --- src/problem_operators/equation_system_problem_operator.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/problem_operators/equation_system_problem_operator.cpp b/src/problem_operators/equation_system_problem_operator.cpp index 0a2a4793c..7d7d944f8 100644 --- a/src/problem_operators/equation_system_problem_operator.cpp +++ b/src/problem_operators/equation_system_problem_operator.cpp @@ -25,8 +25,6 @@ EquationSystemProblemOperator::Update() { GetEquationSystem()->Update(_problem._bc_map, _problem._sources); - // TODO: - rebuild jacobian, jacobian_solver, etc... - ProblemOperator::Update(); } From 2f972989c2f8ee76e89ecc91147bc57b0011effc Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Fri, 24 May 2024 10:59:53 +0000 Subject: [PATCH 28/30] Adds DefaultSolverOptions virtual method Allows ProblemOperatorBase derived classes to set custom default options. --- src/problem_operators/problem_operator_base.cpp | 13 ++++++++++++- src/problem_operators/problem_operator_base.hpp | 16 ++++++++++------ 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index bace1bd33..211cdcb89 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -5,6 +5,16 @@ namespace hephaestus ProblemOperatorBase::ProblemOperatorBase(hephaestus::Problem & problem) : _problem{problem} {} +ProblemOperatorBase::SolverOptions +ProblemOperatorBase::DefaultSolverOptions() const +{ + return {._tolerance = 1e-16, + ._abs_tolerance = 1e-16, + ._max_iteration = 1000, + ._print_level = GetGlobalPrintLevel(), + ._k_dim = 10}; +} + void ProblemOperatorBase::ConstructJacobianSolver() { @@ -121,7 +131,8 @@ void ProblemOperatorBase::Init() { _block_vector = std::make_unique(); - ConstructJacobianSolver(); + _solver_options = DefaultSolverOptions(); + ConstructNonlinearSolver(); SetTrialVariables(); diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 8ed005b1e..85485a3bd 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,16 +19,16 @@ class ProblemOperatorBase /// Update the problem operator after a mesh change. virtual void Update(); - /// A structure for setting solver options. Defaults have been set. + /// A structure for setting solver options. struct SolverOptions { - double _tolerance{1e-16}; - double _abs_tolerance{1e-16}; + double _tolerance; + double _abs_tolerance; - unsigned int _max_iteration{1000}; + unsigned int _max_iteration; - int _print_level{GetGlobalPrintLevel()}; - int _k_dim{10}; + int _print_level; + int _k_dim; }; /// Sets the solver's options. Override in derived classes. @@ -40,6 +40,10 @@ class ProblemOperatorBase /// be possible to use directly. explicit ProblemOperatorBase(hephaestus::Problem & problem); + /// Override in derived classes to set default solver options. These will be + /// applied following solver construction. + virtual SolverOptions DefaultSolverOptions() const; + /// Override in derived classes to construct the Jacobian solver. Called in the /// Init method. virtual void ConstructJacobianSolver(); From 01a2ceb5199339a9f4e290000ab7739a9c545b1d Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Fri, 24 May 2024 11:03:09 +0000 Subject: [PATCH 29/30] Adds ApplySolverOptions; SetSolverOptions no longer virtual and wraps around ApplySolverOptions; adds ConstructJacobianSolverAndApplyOptions with wraps around ConstructJacobianSolver and ensures that solver options are applied. --- .../complex_maxwell_formulation.hpp | 4 +-- src/formulations/Dual/dual_formulation.cpp | 14 +++----- src/formulations/Dual/dual_formulation.hpp | 4 +-- src/formulations/HCurl/hcurl_formulation.cpp | 14 +++----- src/formulations/HCurl/hcurl_formulation.hpp | 4 +-- .../Magnetostatic/statics_formulation.cpp | 14 +++----- .../Magnetostatic/statics_formulation.hpp | 4 +-- .../problem_operator_base.cpp | 32 ++++++++++++------- .../problem_operator_base.hpp | 11 +++++-- 9 files changed, 53 insertions(+), 48 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp index 07239a267..710b31108 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.hpp @@ -75,8 +75,6 @@ class ComplexMaxwellOperator : public ProblemOperator void Update() override; void Solve(mfem::Vector & X) override; - void SetSolverOptions(SolverOptions options) override {} - std::string _h_curl_var_complex_name, _h_curl_var_real_name, _h_curl_var_imag_name, _stiffness_coef_name, _mass_coef_name, _loss_coef_name; @@ -90,6 +88,8 @@ class ComplexMaxwellOperator : public ProblemOperator mfem::Array _ess_bdr_tdofs; protected: + void ApplySolverOptions() override {} + void ConstructJacobianSolver() override; }; diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 46f57341a..fa1bb4bea 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -219,21 +219,17 @@ DualOperator::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - - SetSolverOptions(_solver_options); } void -DualOperator::SetSolverOptions(SolverOptions options) +DualOperator::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(options._tolerance); - solver.SetAbsTol(options._abs_tolerance); - solver.SetMaxIter(options._max_iteration); - solver.SetPrintLevel(options._print_level); - - _solver_options = options; + solver.SetTol(_solver_options._tolerance); + solver.SetAbsTol(_solver_options._abs_tolerance); + solver.SetMaxIter(_solver_options._max_iteration); + solver.SetPrintLevel(_solver_options._print_level); } void diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index bf79c6cae..51988024a 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -61,8 +61,6 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override; - void SetSolverOptions(SolverOptions options) override; - mfem::ParFiniteElementSpace * _h_curl_fe_space{nullptr}; mfem::ParFiniteElementSpace * _h_div_fe_space{nullptr}; @@ -74,6 +72,8 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator protected: void ConstructJacobianSolver() override; + void ApplySolverOptions() override; + int GetSolutionVectorSize() const override; std::unique_ptr _curl; diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index f08457bd7..c9239b1c9 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -89,16 +89,14 @@ HCurlFormulation::RegisterGridFunctions() }; void -HCurlProblemOperator::SetSolverOptions(SolverOptions options) +HCurlProblemOperator::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(options._tolerance); - solver.SetAbsTol(options._abs_tolerance); - solver.SetMaxIter(options._max_iteration); - solver.SetPrintLevel(options._print_level); - - _solver_options = options; + solver.SetTol(_solver_options._tolerance); + solver.SetAbsTol(_solver_options._abs_tolerance); + solver.SetMaxIter(_solver_options._max_iteration); + solver.SetPrintLevel(_solver_options._print_level); } void @@ -114,8 +112,6 @@ HCurlProblemOperator::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - - SetSolverOptions(_solver_options); } CurlCurlEquationSystem::CurlCurlEquationSystem(const hephaestus::InputParameters & params) diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index 20533acf3..25c86b69d 100644 --- a/src/formulations/HCurl/hcurl_formulation.hpp +++ b/src/formulations/HCurl/hcurl_formulation.hpp @@ -37,9 +37,9 @@ class HCurlProblemOperator : public TimeDomainEquationSystemProblemOperator { } - void SetSolverOptions(SolverOptions options) override; - protected: + void ApplySolverOptions() override; + void ConstructJacobianSolver() override; }; diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 156243c3e..bf8e76f7f 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -112,21 +112,17 @@ StaticsOperator::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); - - SetSolverOptions(_solver_options); } void -StaticsOperator::SetSolverOptions(SolverOptions options) +StaticsOperator::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(options._tolerance); - solver.SetMaxIter(options._max_iteration); - solver.SetKDim(options._k_dim); - solver.SetPrintLevel(options._print_level); - - _solver_options = options; + solver.SetTol(_solver_options._tolerance); + solver.SetMaxIter(_solver_options._max_iteration); + solver.SetKDim(_solver_options._k_dim); + solver.SetPrintLevel(_solver_options._print_level); } /* diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index 8e6e5030d..2f3f78eab 100644 --- a/src/formulations/Magnetostatic/statics_formulation.hpp +++ b/src/formulations/Magnetostatic/statics_formulation.hpp @@ -38,9 +38,9 @@ class StaticsOperator : public ProblemOperator void Init() override; void Solve(mfem::Vector & X) override; - void SetSolverOptions(SolverOptions options) override; - protected: + void ApplySolverOptions() override; + void ConstructJacobianSolver() override; private: diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index 211cdcb89..a8b720110 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -26,24 +26,33 @@ ProblemOperatorBase::ConstructJacobianSolver() _jacobian_preconditioner = std::move(precond); _jacobian_solver = std::move(solver); +} - // Set default options or user-options if available. - SetSolverOptions(_solver_options); +void +ProblemOperatorBase::ConstructJacobianSolverAndApplyOptions() +{ + ConstructJacobianSolver(); + ApplySolverOptions(); } void ProblemOperatorBase::SetSolverOptions(SolverOptions options) { - auto & solver = static_cast(*_jacobian_solver); + // Store the options for future. + _solver_options = std::move(options); + ApplySolverOptions(); +} - solver.SetTol(options._tolerance); - solver.SetAbsTol(options._abs_tolerance); - solver.SetMaxIter(options._max_iteration); - solver.SetKDim(options._k_dim); - solver.SetPrintLevel(options._print_level); +void +ProblemOperatorBase::ApplySolverOptions() +{ + auto & solver = static_cast(*_jacobian_solver); - // Store the options for future. - _solver_options = options; + solver.SetTol(_solver_options._tolerance); + solver.SetAbsTol(_solver_options._abs_tolerance); + solver.SetMaxIter(_solver_options._max_iteration); + solver.SetKDim(_solver_options._k_dim); + solver.SetPrintLevel(_solver_options._print_level); } void @@ -133,6 +142,7 @@ ProblemOperatorBase::Init() _block_vector = std::make_unique(); _solver_options = DefaultSolverOptions(); + ConstructJacobianSolverAndApplyOptions(); ConstructNonlinearSolver(); SetTrialVariables(); @@ -152,7 +162,7 @@ ProblemOperatorBase::Update() UpdateBlockVector(*_block_vector); // Update the Jacobian solver. - ConstructJacobianSolver(); + ConstructJacobianSolverAndApplyOptions(); } } \ No newline at end of file diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 85485a3bd..04cc6e21a 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -31,8 +31,8 @@ class ProblemOperatorBase int _k_dim; }; - /// Sets the solver's options. Override in derived classes. - virtual void SetSolverOptions(SolverOptions options); + /// Sets the solver's options. Then calls ApplySolverOptions. + void SetSolverOptions(SolverOptions options); protected: /// Use of protected constructor to only allow construction by derived classes. @@ -44,6 +44,9 @@ class ProblemOperatorBase /// applied following solver construction. virtual SolverOptions DefaultSolverOptions() const; + /// Applies the current solver options. Override in derived classes. + virtual void ApplySolverOptions(); + /// Override in derived classes to construct the Jacobian solver. Called in the /// Init method. virtual void ConstructJacobianSolver(); @@ -103,6 +106,10 @@ class ProblemOperatorBase SolverOptions _solver_options; private: + /// Calls ConstructJacobianSolver followed by ApplySolverOptions. This ensures + /// that ApplySolverOptions is always called! + void ConstructJacobianSolverAndApplyOptions(); + /// Update the block vectors and offsets after a mesh change. void UpdateOffsets(); From 636cdfd34a289bd2623a8ae4c624510df79d7352 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Fri, 24 May 2024 11:16:26 +0000 Subject: [PATCH 30/30] SolverOptions are now private; adds GetSolverOptions accessor. --- src/formulations/Dual/dual_formulation.cpp | 8 ++++---- src/formulations/HCurl/hcurl_formulation.cpp | 8 ++++---- src/formulations/Magnetostatic/statics_formulation.cpp | 8 ++++---- src/problem_operators/problem_operator_base.cpp | 10 +++++----- src/problem_operators/problem_operator_base.hpp | 9 ++++++--- 5 files changed, 23 insertions(+), 20 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index fa1bb4bea..b9c6cfd70 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -226,10 +226,10 @@ DualOperator::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(_solver_options._tolerance); - solver.SetAbsTol(_solver_options._abs_tolerance); - solver.SetMaxIter(_solver_options._max_iteration); - solver.SetPrintLevel(_solver_options._print_level); + solver.SetTol(GetSolverOptions()._tolerance); + solver.SetAbsTol(GetSolverOptions()._abs_tolerance); + solver.SetMaxIter(GetSolverOptions()._max_iteration); + solver.SetPrintLevel(GetSolverOptions()._print_level); } void diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index c9239b1c9..dc0442a0f 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -93,10 +93,10 @@ HCurlProblemOperator::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(_solver_options._tolerance); - solver.SetAbsTol(_solver_options._abs_tolerance); - solver.SetMaxIter(_solver_options._max_iteration); - solver.SetPrintLevel(_solver_options._print_level); + solver.SetTol(GetSolverOptions()._tolerance); + solver.SetAbsTol(GetSolverOptions()._abs_tolerance); + solver.SetMaxIter(GetSolverOptions()._max_iteration); + solver.SetPrintLevel(GetSolverOptions()._print_level); } void diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index bf8e76f7f..4335b3389 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -119,10 +119,10 @@ StaticsOperator::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(_solver_options._tolerance); - solver.SetMaxIter(_solver_options._max_iteration); - solver.SetKDim(_solver_options._k_dim); - solver.SetPrintLevel(_solver_options._print_level); + solver.SetTol(GetSolverOptions()._tolerance); + solver.SetMaxIter(GetSolverOptions()._max_iteration); + solver.SetKDim(GetSolverOptions()._k_dim); + solver.SetPrintLevel(GetSolverOptions()._print_level); } /* diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index a8b720110..965b837db 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -48,11 +48,11 @@ ProblemOperatorBase::ApplySolverOptions() { auto & solver = static_cast(*_jacobian_solver); - solver.SetTol(_solver_options._tolerance); - solver.SetAbsTol(_solver_options._abs_tolerance); - solver.SetMaxIter(_solver_options._max_iteration); - solver.SetKDim(_solver_options._k_dim); - solver.SetPrintLevel(_solver_options._print_level); + solver.SetTol(GetSolverOptions()._tolerance); + solver.SetAbsTol(GetSolverOptions()._abs_tolerance); + solver.SetMaxIter(GetSolverOptions()._max_iteration); + solver.SetKDim(GetSolverOptions()._k_dim); + solver.SetPrintLevel(GetSolverOptions()._print_level); } void diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 04cc6e21a..ae7ce5c13 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -74,6 +74,9 @@ class ProblemOperatorBase /// Returns the number of trial variables. int GetTrialVariablesSize() const; + /// Solver options accessor. + const SolverOptions & GetSolverOptions() const { return _solver_options; } + // Reference to the current problem. hephaestus::Problem & _problem; @@ -102,9 +105,6 @@ class ProblemOperatorBase /// Store the block vector. std::unique_ptr _block_vector{nullptr}; - /// The current solver options. - SolverOptions _solver_options; - private: /// Calls ConstructJacobianSolver followed by ApplySolverOptions. This ensures /// that ApplySolverOptions is always called! @@ -115,5 +115,8 @@ class ProblemOperatorBase /// Update a block vector. Should be called after the offsets have been updated. void UpdateBlockVector(mfem::BlockVector & X); + + /// The current solver options. + SolverOptions _solver_options; }; } \ No newline at end of file