From 2552c43f630ccee514925b0878e1a78e3b0b1b63 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 01/20] Adds ProblemOperatorBase class; ProblemOperatorInterface now only contains pure virtual methods. --- src/problem_operators/problem_operator.hpp | 6 +-- ...nterface.cpp => problem_operator_base.cpp} | 16 ++++---- .../problem_operator_base.hpp | 40 +++++++++++++++++++ .../problem_operator_interface.hpp | 30 ++++---------- .../time_domain_problem_operator.hpp | 7 ++-- 5 files changed, 62 insertions(+), 37 deletions(-) rename src/problem_operators/{problem_operator_interface.cpp => problem_operator_base.cpp} (80%) create mode 100644 src/problem_operators/problem_operator_base.hpp diff --git a/src/problem_operators/problem_operator.hpp b/src/problem_operators/problem_operator.hpp index 1aba00690..1923acbb8 100644 --- a/src/problem_operators/problem_operator.hpp +++ b/src/problem_operators/problem_operator.hpp @@ -2,15 +2,15 @@ #include "../common/pfem_extras.hpp" #include "hephaestus_solvers.hpp" #include "problem_builder_base.hpp" -#include "problem_operator_interface.hpp" +#include "problem_operator_base.hpp" namespace hephaestus { /// Steady-state problem operator with no equation system. -class ProblemOperator : public mfem::Operator, public ProblemOperatorInterface +class ProblemOperator : public mfem::Operator, public ProblemOperatorBase { public: - ProblemOperator(hephaestus::Problem & problem) : ProblemOperatorInterface(problem) {} + ProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {} ~ProblemOperator() override = default; virtual void Solve(mfem::Vector & X) {} diff --git a/src/problem_operators/problem_operator_interface.cpp b/src/problem_operators/problem_operator_base.cpp similarity index 80% rename from src/problem_operators/problem_operator_interface.cpp rename to src/problem_operators/problem_operator_base.cpp index 5fd4d291b..dbffc730d 100644 --- a/src/problem_operators/problem_operator_interface.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -1,10 +1,10 @@ -#include "problem_operator_interface.hpp" +#include "problem_operator_base.hpp" namespace hephaestus { void -ProblemOperatorInterface::SetTrialVariables() +ProblemOperatorBase::SetTrialVariables() { SetTrialVariableNames(); @@ -12,13 +12,13 @@ ProblemOperatorInterface::SetTrialVariables() } void -ProblemOperatorInterface::UpdateOffsets() +ProblemOperatorBase::UpdateOffsets() { UpdateOffsetsWithSize(_trial_variables.size()); } void -ProblemOperatorInterface::UpdateOffsetsWithSize(const size_t block_vector_size) +ProblemOperatorBase::UpdateOffsetsWithSize(const size_t block_vector_size) { if (block_vector_size > _trial_variables.size()) { @@ -54,7 +54,7 @@ ProblemOperatorInterface::UpdateOffsetsWithSize(const size_t block_vector_size) } void -ProblemOperatorInterface::UpdateBlockVector(mfem::BlockVector & X) +ProblemOperatorBase::UpdateBlockVector(mfem::BlockVector & X) { X.Update(_true_offsets); @@ -68,7 +68,7 @@ ProblemOperatorInterface::UpdateBlockVector(mfem::BlockVector & X) } void -ProblemOperatorInterface::Init() +ProblemOperatorBase::Init() { SetTrialVariables(); @@ -77,9 +77,9 @@ ProblemOperatorInterface::Init() } void -ProblemOperatorInterface::Update() +ProblemOperatorBase::Update() { - logger.debug("Update called for ProblemOperatorInterface."); + logger.debug("Update called for ProblemOperatorBase."); // Recalculate the offsets from gridfunction trial variables. UpdateOffsets(); diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp new file mode 100644 index 000000000..6fa7a9b1a --- /dev/null +++ b/src/problem_operators/problem_operator_base.hpp @@ -0,0 +1,40 @@ +#pragma once +#include "problem_operator_interface.hpp" + +namespace hephaestus +{ +class ProblemOperatorBase : public ProblemOperatorInterface +{ +public: + ProblemOperatorBase(hephaestus::Problem & problem) : _problem(problem) {} + ~ProblemOperatorBase() override = default; + + void UpdateOperatorWidthAndHeight() override {} + + void Init() override; + + void Update() override; + + mfem::Array _true_offsets, _block_true_offsets; + + mfem::BlockVector _true_x, _true_rhs; + mfem::OperatorHandle _equation_system_operator; + +protected: + void SetTrialVariableNames() override {} + void SetTrialVariables() override; + void UpdateOffsets() override; + + void UpdateBlockVector(mfem::BlockVector & X) override; + + void UpdateOffsetsWithSize(size_t soln_vector_size) override; + + // Reference to the current problem. + hephaestus::Problem & _problem; + + // Vector of names of state gridfunctions used in formulation, ordered by appearance in block + // vector during solve. + std::vector _trial_var_names; + std::vector _trial_variables; +}; +} \ No newline at end of file diff --git a/src/problem_operators/problem_operator_interface.hpp b/src/problem_operators/problem_operator_interface.hpp index af91e7759..16958904b 100644 --- a/src/problem_operators/problem_operator_interface.hpp +++ b/src/problem_operators/problem_operator_interface.hpp @@ -7,35 +7,21 @@ namespace hephaestus class ProblemOperatorInterface { public: - ProblemOperatorInterface(hephaestus::Problem & problem) : _problem(problem) {} virtual ~ProblemOperatorInterface() = default; - virtual void UpdateOperatorWidthAndHeight() {} + virtual void UpdateOperatorWidthAndHeight() = 0; - virtual void Init(); + virtual void Init() = 0; - virtual void Update(); - - mfem::Array _true_offsets, _block_true_offsets; - - mfem::BlockVector _true_x, _true_rhs; - mfem::OperatorHandle _equation_system_operator; + virtual void Update() = 0; protected: - virtual void SetTrialVariableNames() {} - virtual void SetTrialVariables(); - virtual void UpdateOffsets(); - - virtual void UpdateBlockVector(mfem::BlockVector & X); - - void UpdateOffsetsWithSize(size_t soln_vector_size); + virtual void SetTrialVariableNames() = 0; + virtual void SetTrialVariables() = 0; + virtual void UpdateOffsets() = 0; - // Reference to the current problem. - hephaestus::Problem & _problem; + virtual void UpdateBlockVector(mfem::BlockVector & X) = 0; - // Vector of names of state gridfunctions used in formulation, ordered by appearance in block - // vector during solve. - std::vector _trial_var_names; - std::vector _trial_variables; + virtual void UpdateOffsetsWithSize(size_t soln_vector_size) = 0; }; } \ 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 095bd4fc6..992d4fc9c 100644 --- a/src/problem_operators/time_domain_problem_operator.hpp +++ b/src/problem_operators/time_domain_problem_operator.hpp @@ -2,7 +2,7 @@ #include "../common/pfem_extras.hpp" #include "hephaestus_solvers.hpp" #include "problem_builder_base.hpp" -#include "problem_operator_interface.hpp" +#include "problem_operator_base.hpp" namespace hephaestus { @@ -13,11 +13,10 @@ std::vector GetTimeDerivativeNames(std::vector gridfun /// Problem operator for time-dependent problems with no equation system. The user will need to subclass this since the solve is not /// implemented. -class TimeDomainProblemOperator : public mfem::TimeDependentOperator, - public ProblemOperatorInterface +class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public ProblemOperatorBase { public: - TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorInterface(problem) {} + TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {} ~TimeDomainProblemOperator() override = default; void UpdateOperatorWidthAndHeight() final; From 1b9319124b039ac7884365084bd3c0cd34ac01b0 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 02/20] Removes unused equation system operator handle from ProblemOperatorBase. --- src/problem_operators/problem_operator_base.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 6fa7a9b1a..e82addf7f 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -18,7 +18,6 @@ class ProblemOperatorBase : public ProblemOperatorInterface mfem::Array _true_offsets, _block_true_offsets; mfem::BlockVector _true_x, _true_rhs; - mfem::OperatorHandle _equation_system_operator; protected: void SetTrialVariableNames() override {} From 803284fd2d55a6c4f5f94050c8ae68d73bcbdb75 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 03/20] Adds Width and Height getter/setter to ProblemOperatorInterface; Removes UpdateOperatorWidthAndHeight method. --- src/problem_operators/problem_operator.cpp | 12 ------------ src/problem_operators/problem_operator.hpp | 4 +++- src/problem_operators/problem_operator_base.cpp | 2 +- src/problem_operators/problem_operator_base.hpp | 2 -- src/problem_operators/problem_operator_interface.hpp | 8 ++++++-- .../time_domain_problem_operator.cpp | 6 ------ .../time_domain_problem_operator.hpp | 6 ++++-- 7 files changed, 14 insertions(+), 26 deletions(-) delete mode 100644 src/problem_operators/problem_operator.cpp diff --git a/src/problem_operators/problem_operator.cpp b/src/problem_operators/problem_operator.cpp deleted file mode 100644 index 5e9efa441..000000000 --- a/src/problem_operators/problem_operator.cpp +++ /dev/null @@ -1,12 +0,0 @@ -#include "problem_operator.hpp" - -namespace hephaestus -{ - -void -ProblemOperator::UpdateOperatorWidthAndHeight() -{ - width = height = _true_offsets.Last(); -} - -} // namespace hephaestus \ No newline at end of file diff --git a/src/problem_operators/problem_operator.hpp b/src/problem_operators/problem_operator.hpp index 1923acbb8..ace77b42e 100644 --- a/src/problem_operators/problem_operator.hpp +++ b/src/problem_operators/problem_operator.hpp @@ -16,7 +16,9 @@ class ProblemOperator : public mfem::Operator, public ProblemOperatorBase virtual void Solve(mfem::Vector & X) {} void Mult(const mfem::Vector & x, mfem::Vector & y) const override {} - void UpdateOperatorWidthAndHeight() final; +protected: + int & Width() final { return mfem::Operator::width; } + int & Height() final { return mfem::Operator::height; } }; } // namespace hephaestus \ No newline at end of file diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index dbffc730d..2def6e993 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -50,7 +50,7 @@ ProblemOperatorBase::UpdateOffsetsWithSize(const size_t block_vector_size) _true_rhs.Update(_block_true_offsets); // Set width and height member variables. - UpdateOperatorWidthAndHeight(); + Width() = Height() = _true_offsets.Last(); } void diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index e82addf7f..356505ab9 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -9,8 +9,6 @@ class ProblemOperatorBase : public ProblemOperatorInterface ProblemOperatorBase(hephaestus::Problem & problem) : _problem(problem) {} ~ProblemOperatorBase() override = default; - void UpdateOperatorWidthAndHeight() override {} - void Init() override; void Update() override; diff --git a/src/problem_operators/problem_operator_interface.hpp b/src/problem_operators/problem_operator_interface.hpp index 16958904b..d5f59ebb4 100644 --- a/src/problem_operators/problem_operator_interface.hpp +++ b/src/problem_operators/problem_operator_interface.hpp @@ -9,8 +9,6 @@ class ProblemOperatorInterface public: virtual ~ProblemOperatorInterface() = default; - virtual void UpdateOperatorWidthAndHeight() = 0; - virtual void Init() = 0; virtual void Update() = 0; @@ -23,5 +21,11 @@ class ProblemOperatorInterface virtual void UpdateBlockVector(mfem::BlockVector & X) = 0; virtual void UpdateOffsetsWithSize(size_t soln_vector_size) = 0; + + /// Returns a reference to the operator's width. + virtual int & Width() = 0; + + /// Returns a reference to the operator's height. + virtual int & Height() = 0; }; } \ No newline at end of file diff --git a/src/problem_operators/time_domain_problem_operator.cpp b/src/problem_operators/time_domain_problem_operator.cpp index ecde0de8c..0c02bcb36 100644 --- a/src/problem_operators/time_domain_problem_operator.cpp +++ b/src/problem_operators/time_domain_problem_operator.cpp @@ -20,10 +20,4 @@ GetTimeDerivativeNames(std::vector gridfunction_names) return time_derivative_names; } -void -TimeDomainProblemOperator::UpdateOperatorWidthAndHeight() -{ - width = height = _true_offsets.Last(); -} - } // 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 992d4fc9c..78804d313 100644 --- a/src/problem_operators/time_domain_problem_operator.hpp +++ b/src/problem_operators/time_domain_problem_operator.hpp @@ -19,9 +19,11 @@ class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public Pro TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {} ~TimeDomainProblemOperator() override = default; - void UpdateOperatorWidthAndHeight() final; - void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override {} + +protected: + int & Width() final { return mfem::TimeDependentOperator::width; } + int & Height() final { return mfem::TimeDependentOperator::height; } }; } // namespace hephaestus \ No newline at end of file From 21e5c7528beaa7a963f6805ab6f4154d8ddd8747 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 04/20] Removes ProblemOperatorInterface; now using base class ProblemOperatorBase. --- src/problem_builders/problem_builder_base.hpp | 4 +-- .../equation_system_problem_operator.hpp | 2 +- .../problem_operator_base.hpp | 29 +++++++++++------ .../problem_operator_interface.hpp | 31 ------------------- ...omain_equation_system_problem_operator.hpp | 2 +- 5 files changed, 23 insertions(+), 45 deletions(-) delete mode 100644 src/problem_operators/problem_operator_interface.hpp diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 55d881def..aec0484e7 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -12,7 +12,7 @@ namespace hephaestus { /// Forwards declaration. -class ProblemOperatorInterface; +class ProblemOperatorBase; /// Base problem class. class Problem @@ -46,7 +46,7 @@ class Problem int _num_procs; /// Returns a pointer to the operator. See derived classes. - [[nodiscard]] virtual hephaestus::ProblemOperatorInterface * GetOperator() const = 0; + [[nodiscard]] virtual hephaestus::ProblemOperatorBase * GetOperator() const = 0; /// Virtual method to construct the operator. Call for default problems. virtual void ConstructOperator() = 0; diff --git a/src/problem_operators/equation_system_problem_operator.hpp b/src/problem_operators/equation_system_problem_operator.hpp index 5e4ec50e6..5f583b4e8 100644 --- a/src/problem_operators/equation_system_problem_operator.hpp +++ b/src/problem_operators/equation_system_problem_operator.hpp @@ -1,6 +1,6 @@ #pragma once #include "problem_operator.hpp" -#include "problem_operator_interface.hpp" +#include "problem_operator_base.hpp" #include "equation_system_interface.hpp" namespace hephaestus diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 356505ab9..7c2ee81c0 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -1,30 +1,39 @@ #pragma once -#include "problem_operator_interface.hpp" +#include "mfem.hpp" +#include "problem_builder_base.hpp" +#include +#include namespace hephaestus { -class ProblemOperatorBase : public ProblemOperatorInterface +class ProblemOperatorBase { public: ProblemOperatorBase(hephaestus::Problem & problem) : _problem(problem) {} - ~ProblemOperatorBase() override = default; + virtual ~ProblemOperatorBase() = default; - void Init() override; + virtual void Init(); - void Update() override; + virtual void Update(); mfem::Array _true_offsets, _block_true_offsets; mfem::BlockVector _true_x, _true_rhs; protected: - void SetTrialVariableNames() override {} - void SetTrialVariables() override; - void UpdateOffsets() override; + virtual void SetTrialVariableNames() {} + virtual void SetTrialVariables(); + virtual void UpdateOffsets(); - void UpdateBlockVector(mfem::BlockVector & X) override; + virtual void UpdateBlockVector(mfem::BlockVector & X); - void UpdateOffsetsWithSize(size_t soln_vector_size) override; + virtual void UpdateOffsetsWithSize(size_t soln_vector_size); + + /// Returns a reference to the operator's width. + virtual int & Width() = 0; + + /// Returns a reference to the operator's height. + virtual int & Height() = 0; // Reference to the current problem. hephaestus::Problem & _problem; diff --git a/src/problem_operators/problem_operator_interface.hpp b/src/problem_operators/problem_operator_interface.hpp deleted file mode 100644 index d5f59ebb4..000000000 --- a/src/problem_operators/problem_operator_interface.hpp +++ /dev/null @@ -1,31 +0,0 @@ -#pragma once -#include "problem_builder_base.hpp" - -namespace hephaestus -{ -/// Interface inherited by ProblemOperator and TimeDomainProblemOperator. Removes duplicated code in both classes. -class ProblemOperatorInterface -{ -public: - virtual ~ProblemOperatorInterface() = default; - - virtual void Init() = 0; - - virtual void Update() = 0; - -protected: - virtual void SetTrialVariableNames() = 0; - virtual void SetTrialVariables() = 0; - virtual void UpdateOffsets() = 0; - - virtual void UpdateBlockVector(mfem::BlockVector & X) = 0; - - virtual void UpdateOffsetsWithSize(size_t soln_vector_size) = 0; - - /// Returns a reference to the operator's width. - virtual int & Width() = 0; - - /// Returns a reference to the operator's height. - virtual int & Height() = 0; -}; -} \ No newline at end of file diff --git a/src/problem_operators/time_domain_equation_system_problem_operator.hpp b/src/problem_operators/time_domain_equation_system_problem_operator.hpp index 3aa518753..5e83e6cbd 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.hpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.hpp @@ -1,7 +1,7 @@ #pragma once #include "../common/pfem_extras.hpp" #include "time_domain_problem_operator.hpp" -#include "problem_operator_interface.hpp" +#include "problem_operator_base.hpp" #include "equation_system_interface.hpp" namespace hephaestus From 3562141911f7ba95626e96b3b25f1e439e3c810c Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 05/20] Added comments to ProblemOperatorBase methods. --- .../problem_operator_base.hpp | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 7c2ee81c0..3e3bc17e9 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -10,23 +10,35 @@ class ProblemOperatorBase { public: ProblemOperatorBase(hephaestus::Problem & problem) : _problem(problem) {} + virtual ~ProblemOperatorBase() = default; + /// Initialize the problem operator. virtual void Init(); + /// Update the problem operator after a mesh change. virtual void Update(); + /// Arrays of offsets. mfem::Array _true_offsets, _block_true_offsets; + /// Block vectors. mfem::BlockVector _true_x, _true_rhs; protected: + /// Set trial variables names. Override in derived classes. virtual void SetTrialVariableNames() {} + + /// Set gridfunction of trial variables from the trial variable names. virtual void SetTrialVariables(); + + /// Update the block vectors and offsets after a mesh change. virtual void UpdateOffsets(); + /// Update a block vector. Should be called after the offsets have been updated. virtual void UpdateBlockVector(mfem::BlockVector & X); + /// Called by UpdateOffsets. virtual void UpdateOffsetsWithSize(size_t soln_vector_size); /// Returns a reference to the operator's width. @@ -38,9 +50,11 @@ class ProblemOperatorBase // Reference to the current problem. hephaestus::Problem & _problem; - // Vector of names of state gridfunctions used in formulation, ordered by appearance in block - // vector during solve. + /// Vector of names of state gridfunctions used in formulation, ordered by + /// appearance in block vector during solve. std::vector _trial_var_names; + + /// Vector of trial variables set by SetTrialVariables. std::vector _trial_variables; }; } \ No newline at end of file From 19f10a4c967bb46f9b7e1c93155447a0ce0edd19 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 06/20] Deleted default ProblemOperatorBase constructor and made constructor protected to prevent direct construction of the base class. --- src/problem_operators/problem_operator_base.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 3e3bc17e9..22f644bd9 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -9,7 +9,7 @@ namespace hephaestus class ProblemOperatorBase { public: - ProblemOperatorBase(hephaestus::Problem & problem) : _problem(problem) {} + ProblemOperatorBase() = delete; virtual ~ProblemOperatorBase() = default; @@ -26,6 +26,11 @@ class ProblemOperatorBase mfem::BlockVector _true_x, _true_rhs; 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) : _problem(problem) {} + /// Set trial variables names. Override in derived classes. virtual void SetTrialVariableNames() {} From 0bffc6992ca3b235046fdba92f24c28f014fc3a0 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 07/20] Made UpdateOffsetsWithSize non-virtual. --- .../problem_operator_base.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 22f644bd9..a5d241cf4 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(); - /// Arrays of offsets. - mfem::Array _true_offsets, _block_true_offsets; - - /// Block vectors. - mfem::BlockVector _true_x, _true_rhs; - 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 @@ -43,15 +37,15 @@ class ProblemOperatorBase /// Update a block vector. Should be called after the offsets have been updated. virtual void UpdateBlockVector(mfem::BlockVector & X); - /// Called by UpdateOffsets. - virtual void UpdateOffsetsWithSize(size_t soln_vector_size); - /// Returns a reference to the operator's width. virtual int & Width() = 0; /// Returns a reference to the operator's height. virtual int & Height() = 0; + /// Called by UpdateOffsets. + void UpdateOffsetsWithSize(size_t soln_vector_size); + // Reference to the current problem. hephaestus::Problem & _problem; @@ -61,5 +55,11 @@ class ProblemOperatorBase /// Vector of trial variables set by SetTrialVariables. std::vector _trial_variables; + + /// Arrays of offsets. + mfem::Array _true_offsets, _block_true_offsets; + + /// Block vectors. + mfem::BlockVector _true_x, _true_rhs; }; } \ No newline at end of file From 8d57a1846594e6f26f563297d6bd5193207bd6f9 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 08/20] Removed UpdateOffsetsWithSize; added GetSolutionVectorSize and GetTrialVariablesSize methods. --- src/formulations/Dual/dual_formulation.cpp | 6 ++-- src/formulations/Dual/dual_formulation.hpp | 2 +- .../problem_operator_base.cpp | 30 +++++++++++-------- .../problem_operator_base.hpp | 8 +++-- 4 files changed, 28 insertions(+), 18 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index e85002533..6b8ddb936 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -238,12 +238,12 @@ DualOperator::ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vecto logger.info("{} ImplicitSolve: {} seconds", typeid(this).name(), sw); } -void -DualOperator::UpdateOffsets() +int +DualOperator::GetSolutionVectorSize() const { // Blocks for solution vector are smaller than the operator size for DualOperator, // as curl is stored separately. Block operator only has the HCurl TrueVSize; - return UpdateOffsetsWithSize(_trial_variables.size() - 1); + return static_cast(GetTrialVariablesSize() - 1); } } // namespace hephaestus diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index 58d55005d..78dd23a08 100644 --- a/src/formulations/Dual/dual_formulation.hpp +++ b/src/formulations/Dual/dual_formulation.hpp @@ -74,7 +74,7 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator mfem::ParGridFunction * _dv{nullptr}; // HDiv vector field protected: - void UpdateOffsets() override; + int GetSolutionVectorSize() const override; std::unique_ptr _curl; }; diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index 2def6e993..ba3887bd7 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -11,31 +11,37 @@ ProblemOperatorBase::SetTrialVariables() _trial_variables = _problem._gridfunctions.Get(_trial_var_names); } -void -ProblemOperatorBase::UpdateOffsets() +int +ProblemOperatorBase::GetSolutionVectorSize() const +{ + return static_cast(GetTrialVariablesSize()); +} + +int +ProblemOperatorBase::GetTrialVariablesSize() const { - UpdateOffsetsWithSize(_trial_variables.size()); + return static_cast(_trial_variables.size()); } void -ProblemOperatorBase::UpdateOffsetsWithSize(const size_t block_vector_size) +ProblemOperatorBase::UpdateOffsets() { - if (block_vector_size > _trial_variables.size()) + if (GetSolutionVectorSize() > GetTrialVariablesSize()) { - MFEM_ABORT("Solution vector size (" << block_vector_size + MFEM_ABORT("Solution vector size (" << GetSolutionVectorSize() << ") cannot exceed the number of trial variables (" - << _trial_variables.size() << ")."); + << GetTrialVariablesSize() << ")."); } - _block_true_offsets.SetSize(block_vector_size + 1); - _true_offsets.SetSize(_trial_variables.size() + 1); + _block_true_offsets.SetSize(GetSolutionVectorSize() + 1); + _true_offsets.SetSize(GetTrialVariablesSize() + 1); _block_true_offsets[0] = _true_offsets[0] = 0; - for (size_t i = 0; i < _trial_variables.size(); ++i) + for (int i = 0; i < GetTrialVariablesSize(); ++i) { mfem::ParFiniteElementSpace * fespace = _trial_variables.at(i)->ParFESpace(); - if (i < block_vector_size) + if (i < GetSolutionVectorSize()) _block_true_offsets[i + 1] = fespace->GetTrueVSize(); _true_offsets[i + 1] = fespace->GetVSize(); @@ -58,7 +64,7 @@ ProblemOperatorBase::UpdateBlockVector(mfem::BlockVector & X) { X.Update(_true_offsets); - for (size_t i = 0; i < _trial_variables.size(); ++i) + for (int i = 0; i < GetTrialVariablesSize(); ++i) { mfem::ParGridFunction * trial_var = _trial_variables.at(i); diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index a5d241cf4..3a2433f54 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -43,8 +43,12 @@ class ProblemOperatorBase /// Returns a reference to the operator's height. virtual int & Height() = 0; - /// Called by UpdateOffsets. - void UpdateOffsetsWithSize(size_t soln_vector_size); + /// Returns the solution vector size. By default this will be the same as the + /// number of trial variables. + virtual int GetSolutionVectorSize() const; + + /// Returns the number of trial variables. + int GetTrialVariablesSize() const; // Reference to the current problem. hephaestus::Problem & _problem; From c76733b073d650af754a8d3b2a37750f78d8fdef Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 14:10:19 +0000 Subject: [PATCH 09/20] UpdateOffsets and UpdateBlockVector are now private methods; no longer virtual --- src/problem_operators/problem_operator_base.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 3a2433f54..159fac372 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -31,12 +31,6 @@ class ProblemOperatorBase /// Set gridfunction of trial variables from the trial variable names. virtual void SetTrialVariables(); - /// Update the block vectors and offsets after a mesh change. - virtual void UpdateOffsets(); - - /// Update a block vector. Should be called after the offsets have been updated. - virtual void UpdateBlockVector(mfem::BlockVector & X); - /// Returns a reference to the operator's width. virtual int & Width() = 0; @@ -65,5 +59,12 @@ class ProblemOperatorBase /// Block vectors. mfem::BlockVector _true_x, _true_rhs; + +private: + /// Update the block vectors and offsets after a mesh change. + void UpdateOffsets(); + + /// Update a block vector. Should be called after the offsets have been updated. + void UpdateBlockVector(mfem::BlockVector & X); }; } \ No newline at end of file From b89084020a2f794838305f8ab2c924c36475d791 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 15:08:16 +0000 Subject: [PATCH 10/20] Jacobian solver is now stored inside ProblemOperatorBase. --- .../ComplexMaxwell/complex_maxwell_formulation.cpp | 4 ++-- .../Magnetostatic/statics_formulation.cpp | 4 ++-- src/problem_builders/problem_builder_base.cpp | 12 ++++++------ src/problem_builders/problem_builder_base.hpp | 1 - src/problem_operators/problem_operator_base.hpp | 6 ++++++ .../time_domain_equation_system_problem_operator.cpp | 6 +++--- 6 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp index d4bcbdd0e..238ca0c36 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp @@ -218,8 +218,8 @@ ComplexMaxwellOperator::Solve(mfem::Vector & X) auto * jac_z = jac.As()->GetSystemMatrix(); - _problem._jacobian_solver->SetOperator(*jac_z); - _problem._jacobian_solver->Mult(rhs, u); + _jacobian_solver->SetOperator(*jac_z); + _jacobian_solver->Mult(rhs, u); sqlf.RecoverFEMSolution(u, lf, *_u); _problem._gridfunctions.GetRef(_trial_var_names.at(0)) = _u->real(); diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 05e0f0d08..64049b52c 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -152,8 +152,8 @@ StaticsOperator::Solve(mfem::Vector & X) // Define and apply a parallel FGMRES solver for AX=B with the AMS // preconditioner from hypre. - _problem._jacobian_solver->SetOperator(curl_mu_inv_curl); - _problem._jacobian_solver->Mult(rhs_tdofs, sol_tdofs); + _jacobian_solver->SetOperator(curl_mu_inv_curl); + _jacobian_solver->Mult(rhs_tdofs, sol_tdofs); blf.RecoverFEMSolution(sol_tdofs, lf, gf); logger.info("{} Solve: {} seconds", typeid(this).name(), sw); diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 49c40900b..c0d9e861c 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -105,7 +105,7 @@ ProblemBuilder::SetJacobianPreconditioner(std::shared_ptr precondi void ProblemBuilder::SetJacobianSolver(std::shared_ptr jacobian_solver) { - GetProblem()->_jacobian_solver = jacobian_solver; + GetProblem()->GetOperator()->JacobianSolver() = jacobian_solver; } void @@ -275,7 +275,7 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams if (preconditioner) solver->SetPreconditioner(*preconditioner); - GetProblem()->_jacobian_solver = solver; + GetProblem()->GetOperator()->JacobianSolver() = solver; break; } case SolverType::HYPRE_GMRES: @@ -291,7 +291,7 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams if (preconditioner) solver->SetPreconditioner(*preconditioner); - GetProblem()->_jacobian_solver = solver; + GetProblem()->GetOperator()->JacobianSolver() = solver; break; } case SolverType::HYPRE_FGMRES: @@ -306,7 +306,7 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams if (preconditioner) solver->SetPreconditioner(*preconditioner); - GetProblem()->_jacobian_solver = solver; + GetProblem()->GetOperator()->JacobianSolver() = solver; break; } case SolverType::HYPRE_AMG: @@ -317,14 +317,14 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams solver->SetMaxIter(max_iter); solver->SetPrintLevel(print_level); - GetProblem()->_jacobian_solver = solver; + GetProblem()->GetOperator()->JacobianSolver() = solver; break; } case SolverType::SUPER_LU: { auto solver = std::make_shared(GetProblem()->_comm); - GetProblem()->_jacobian_solver = solver; + GetProblem()->GetOperator()->JacobianSolver() = solver; break; } default: diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index aec0484e7..ecc70b4d0 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -34,7 +34,6 @@ class Problem std::unique_ptr _f{nullptr}; std::shared_ptr _jacobian_preconditioner{nullptr}; - std::shared_ptr _jacobian_solver{nullptr}; std::shared_ptr _nonlinear_solver{nullptr}; hephaestus::FECollections _fecs; diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 159fac372..d0a7779f8 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,6 +19,9 @@ class ProblemOperatorBase /// Update the problem operator after a mesh change. virtual void Update(); + /// Returns a reference to the Jacobian solver which allows it to be set. NB: bad code! + std::shared_ptr & JacobianSolver() { return _jacobian_solver; } + 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 @@ -60,6 +63,9 @@ class ProblemOperatorBase /// Block vectors. mfem::BlockVector _true_x, _true_rhs; + /// Store the Jacobian solver. + std::shared_ptr _jacobian_solver{nullptr}; + private: /// Update the block vectors and offsets after a mesh change. void UpdateOffsets(); 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 94b6383f9..467b3f3c4 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -63,13 +63,13 @@ TimeDomainEquationSystemProblemOperator::Update() _problem._jacobian_preconditioner = precond; // Set new preconditioner. - std::static_pointer_cast(_problem._jacobian_solver) + std::static_pointer_cast(_jacobian_solver) ->SetPreconditioner( *std::static_pointer_cast(_problem._jacobian_preconditioner)); // Set Jacobian matrix. auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); - _problem._jacobian_solver->SetOperator(*matrix); + _jacobian_solver->SetOperator(*matrix); } } @@ -89,7 +89,7 @@ TimeDomainEquationSystemProblemOperator::ImplicitSolve(const double dt, _problem._coefficients.SetTime(GetTime()); BuildEquationSystemOperator(dt); - _problem._nonlinear_solver->SetSolver(*_problem._jacobian_solver); + _problem._nonlinear_solver->SetSolver(*_jacobian_solver); _problem._nonlinear_solver->SetOperator(*GetEquationSystem()); _problem._nonlinear_solver->Mult(_true_rhs, _true_x); From 9c6e1fbbec6bd8ab6e5c4fb56475ec7fddead968 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 15:16:58 +0000 Subject: [PATCH 11/20] Moves the Jacobian preconditioner into ProblemOperatorBase. --- src/formulations/Dual/dual_formulation.cpp | 2 +- src/formulations/HCurl/hcurl_formulation.cpp | 2 +- src/formulations/Magnetostatic/statics_formulation.cpp | 2 +- src/problem_builders/problem_builder_base.cpp | 8 ++++---- src/problem_builders/problem_builder_base.hpp | 1 - src/problem_operators/problem_operator_base.hpp | 6 ++++++ .../time_domain_equation_system_problem_operator.cpp | 5 ++--- 7 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 6b8ddb936..4f9868afa 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -67,7 +67,7 @@ DualFormulation::ConstructJacobianPreconditioner() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->_jacobian_preconditioner = precond; + GetProblem()->GetOperator()->JacobianPreconditioner() = precond; } void diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 435e932f7..2be335051 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -75,7 +75,7 @@ HCurlFormulation::ConstructJacobianPreconditioner() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->_jacobian_preconditioner = precond; + GetProblem()->GetOperator()->JacobianPreconditioner() = precond; } void diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 64049b52c..8232d9458 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -46,7 +46,7 @@ StaticsFormulation::ConstructJacobianPreconditioner() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->_jacobian_preconditioner = precond; + GetProblem()->GetOperator()->JacobianPreconditioner() = precond; } void diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index c0d9e861c..a5d076250 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -99,7 +99,7 @@ ProblemBuilder::SetSolverOptions(hephaestus::InputParameters & solver_options) void ProblemBuilder::SetJacobianPreconditioner(std::shared_ptr preconditioner) { - GetProblem()->_jacobian_preconditioner = preconditioner; + GetProblem()->GetOperator()->JacobianPreconditioner() = preconditioner; } void @@ -228,7 +228,7 @@ ProblemBuilder::ConstructJacobianPreconditioner() auto precond = std::make_shared(); precond->SetPrintLevel(GetGlobalPrintLevel()); - GetProblem()->_jacobian_preconditioner = precond; + GetProblem()->GetOperator()->JacobianPreconditioner() = precond; } void @@ -258,8 +258,8 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams solver_options.GetOptionalParam("PrintLevel", default_params._print_level); const auto k_dim = solver_options.GetOptionalParam("KDim", default_params._k_dim); - auto preconditioner = - std::dynamic_pointer_cast(GetProblem()->_jacobian_preconditioner); + auto preconditioner = std::dynamic_pointer_cast( + GetProblem()->GetOperator()->JacobianPreconditioner()); switch (type) { diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index ecc70b4d0..db8a8d8e3 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -33,7 +33,6 @@ class Problem std::unique_ptr _ode_solver{nullptr}; std::unique_ptr _f{nullptr}; - std::shared_ptr _jacobian_preconditioner{nullptr}; std::shared_ptr _nonlinear_solver{nullptr}; hephaestus::FECollections _fecs; diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index d0a7779f8..26e448fe7 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -22,6 +22,9 @@ class ProblemOperatorBase /// Returns a reference to the Jacobian solver which allows it to be set. NB: bad code! std::shared_ptr & JacobianSolver() { return _jacobian_solver; } + /// Returns a reference to the Jacobian preconditioner which allows it to be set. + std::shared_ptr & JacobianPreconditioner() { return _jacobian_preconditioner; } + 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 @@ -63,6 +66,9 @@ class ProblemOperatorBase /// Block vectors. mfem::BlockVector _true_x, _true_rhs; + /// Store the Jacobian preconditioner. + std::shared_ptr _jacobian_preconditioner{nullptr}; + /// Store the Jacobian solver. std::shared_ptr _jacobian_solver{nullptr}; 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 467b3f3c4..86ad86c93 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -60,12 +60,11 @@ TimeDomainEquationSystemProblemOperator::Update() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - _problem._jacobian_preconditioner = precond; + _jacobian_preconditioner = precond; // Set new preconditioner. std::static_pointer_cast(_jacobian_solver) - ->SetPreconditioner( - *std::static_pointer_cast(_problem._jacobian_preconditioner)); + ->SetPreconditioner(*std::static_pointer_cast(_jacobian_preconditioner)); // Set Jacobian matrix. auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); From db4e8c62990a793371c848402159ddd6de76621f Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 15:21:41 +0000 Subject: [PATCH 12/20] Moved the nonlinear solver into ProblemOperatorBase. --- src/problem_builders/problem_builder_base.cpp | 2 +- src/problem_builders/problem_builder_base.hpp | 2 -- src/problem_operators/problem_operator_base.hpp | 6 ++++++ .../time_domain_equation_system_problem_operator.cpp | 6 +++--- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index a5d076250..079e787fc 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -345,7 +345,7 @@ ProblemBuilder::ConstructNonlinearSolver() nl_solver->SetAbsTol(0.0); nl_solver->SetMaxIter(1); - GetProblem()->_nonlinear_solver = nl_solver; + GetProblem()->GetOperator()->NonlinearSolver() = nl_solver; } void diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index db8a8d8e3..fece81b63 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -33,8 +33,6 @@ class Problem std::unique_ptr _ode_solver{nullptr}; std::unique_ptr _f{nullptr}; - std::shared_ptr _nonlinear_solver{nullptr}; - hephaestus::FECollections _fecs; hephaestus::FESpaces _fespaces; hephaestus::GridFunctions _gridfunctions; diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 26e448fe7..16b54c3e4 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -25,6 +25,9 @@ class ProblemOperatorBase /// Returns a reference to the Jacobian preconditioner which allows it to be set. std::shared_ptr & JacobianPreconditioner() { return _jacobian_preconditioner; } + /// Returns a reference to the non-linear solver which allows it to be set. + std::shared_ptr & NonlinearSolver() { return _nonlinear_solver; } + 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 @@ -72,6 +75,9 @@ class ProblemOperatorBase /// Store the Jacobian solver. std::shared_ptr _jacobian_solver{nullptr}; + /// Store the non-linear solver. + std::shared_ptr _nonlinear_solver{nullptr}; + private: /// Update the block vectors and offsets after a mesh change. void UpdateOffsets(); 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 86ad86c93..655ed0d1e 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -88,9 +88,9 @@ TimeDomainEquationSystemProblemOperator::ImplicitSolve(const double dt, _problem._coefficients.SetTime(GetTime()); BuildEquationSystemOperator(dt); - _problem._nonlinear_solver->SetSolver(*_jacobian_solver); - _problem._nonlinear_solver->SetOperator(*GetEquationSystem()); - _problem._nonlinear_solver->Mult(_true_rhs, _true_x); + _nonlinear_solver->SetSolver(*_jacobian_solver); + _nonlinear_solver->SetOperator(*GetEquationSystem()); + _nonlinear_solver->Mult(_true_rhs, _true_x); GetEquationSystem()->RecoverFEMSolution(_true_x, _problem._gridfunctions); } From 661adaa76368cbefbd2e469ce61baae963127e6c Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 15:45:02 +0000 Subject: [PATCH 13/20] Switched to unique pointers for jacobian solvers; Added proper setters. --- src/formulations/Dual/dual_formulation.cpp | 4 +- src/formulations/HCurl/hcurl_formulation.cpp | 4 +- .../Magnetostatic/statics_formulation.cpp | 6 +-- src/problem_builders/problem_builder_base.cpp | 40 +++++++++---------- src/problem_builders/problem_builder_base.hpp | 4 +- .../problem_operator_base.hpp | 28 ++++++++----- ...omain_equation_system_problem_operator.cpp | 8 ++-- 7 files changed, 52 insertions(+), 42 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 4f9868afa..785634773 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -62,12 +62,12 @@ void DualFormulation::ConstructJacobianPreconditioner() { auto precond = - std::make_shared(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); + std::make_unique(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->GetOperator()->JacobianPreconditioner() = precond; + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); } void diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 2be335051..4c4110256 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -70,12 +70,12 @@ void HCurlFormulation::ConstructJacobianPreconditioner() { auto precond = - std::make_shared(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); + std::make_unique(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->GetOperator()->JacobianPreconditioner() = precond; + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); } void diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 8232d9458..3d06d1fe8 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -40,13 +40,13 @@ StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string void StaticsFormulation::ConstructJacobianPreconditioner() { - std::shared_ptr precond{std::make_shared( - GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace())}; + auto precond = std::make_unique( + GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace()); precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->GetOperator()->JacobianPreconditioner() = precond; + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); } void diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 079e787fc..ba3a0eca3 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -97,15 +97,15 @@ ProblemBuilder::SetSolverOptions(hephaestus::InputParameters & solver_options) } void -ProblemBuilder::SetJacobianPreconditioner(std::shared_ptr preconditioner) +ProblemBuilder::SetJacobianPreconditioner(std::unique_ptr preconditioner) { - GetProblem()->GetOperator()->JacobianPreconditioner() = preconditioner; + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(preconditioner)); } void -ProblemBuilder::SetJacobianSolver(std::shared_ptr jacobian_solver) +ProblemBuilder::SetJacobianSolver(std::unique_ptr jacobian_solver) { - GetProblem()->GetOperator()->JacobianSolver() = jacobian_solver; + GetProblem()->GetOperator()->SetJacobianSolver(std::move(jacobian_solver)); } void @@ -225,10 +225,10 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr(); + auto precond = std::make_unique(); precond->SetPrintLevel(GetGlobalPrintLevel()); - GetProblem()->GetOperator()->JacobianPreconditioner() = precond; + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); } void @@ -258,14 +258,14 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams solver_options.GetOptionalParam("PrintLevel", default_params._print_level); const auto k_dim = solver_options.GetOptionalParam("KDim", default_params._k_dim); - auto preconditioner = std::dynamic_pointer_cast( - GetProblem()->GetOperator()->JacobianPreconditioner()); + auto preconditioner = + dynamic_cast(GetProblem()->GetOperator()->JacobianPreconditioner()); switch (type) { case SolverType::HYPRE_PCG: { - auto solver = std::make_shared(GetProblem()->_comm); + auto solver = std::make_unique(GetProblem()->_comm); solver->SetTol(tolerance); solver->SetAbsTol(abs_tolerance); @@ -275,12 +275,12 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams if (preconditioner) solver->SetPreconditioner(*preconditioner); - GetProblem()->GetOperator()->JacobianSolver() = solver; + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); break; } case SolverType::HYPRE_GMRES: { - auto solver = std::make_shared(GetProblem()->_comm); + auto solver = std::make_unique(GetProblem()->_comm); solver->SetTol(tolerance); solver->SetAbsTol(abs_tolerance); @@ -291,12 +291,12 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams if (preconditioner) solver->SetPreconditioner(*preconditioner); - GetProblem()->GetOperator()->JacobianSolver() = solver; + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); break; } case SolverType::HYPRE_FGMRES: { - auto solver = std::make_shared(GetProblem()->_comm); + auto solver = std::make_unique(GetProblem()->_comm); solver->SetTol(tolerance); solver->SetMaxIter(max_iter); @@ -306,25 +306,25 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams if (preconditioner) solver->SetPreconditioner(*preconditioner); - GetProblem()->GetOperator()->JacobianSolver() = solver; + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); break; } case SolverType::HYPRE_AMG: { - auto solver = std::make_shared(); + auto solver = std::make_unique(); solver->SetTol(tolerance); solver->SetMaxIter(max_iter); solver->SetPrintLevel(print_level); - GetProblem()->GetOperator()->JacobianSolver() = solver; + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); break; } case SolverType::SUPER_LU: { - auto solver = std::make_shared(GetProblem()->_comm); + auto solver = std::make_unique(GetProblem()->_comm); - GetProblem()->GetOperator()->JacobianSolver() = solver; + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); break; } default: @@ -338,14 +338,14 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams void ProblemBuilder::ConstructNonlinearSolver() { - auto nl_solver = std::make_shared(GetProblem()->_comm); + 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()->NonlinearSolver() = nl_solver; + GetProblem()->GetOperator()->SetNonlinearSolver(std::move(nl_solver)); } void diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index fece81b63..7aa800551 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -76,8 +76,8 @@ class ProblemBuilder void SetSources(hephaestus::Sources & sources); void SetOutputs(hephaestus::Outputs & outputs); void SetSolverOptions(hephaestus::InputParameters & solver_options); - void SetJacobianPreconditioner(std::shared_ptr preconditioner); - void SetJacobianSolver(std::shared_ptr solver); + 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 16b54c3e4..5f5269a96 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,14 +19,24 @@ class ProblemOperatorBase /// Update the problem operator after a mesh change. virtual void Update(); - /// Returns a reference to the Jacobian solver which allows it to be set. NB: bad code! - std::shared_ptr & JacobianSolver() { return _jacobian_solver; } + void SetJacobianPreconditioner(std::unique_ptr preconditioner) + { + _jacobian_preconditioner = std::move(preconditioner); + } - /// Returns a reference to the Jacobian preconditioner which allows it to be set. - std::shared_ptr & JacobianPreconditioner() { return _jacobian_preconditioner; } + void SetJacobianSolver(std::unique_ptr solver) + { + _jacobian_solver = std::move(solver); + } - /// Returns a reference to the non-linear solver which allows it to be set. - std::shared_ptr & NonlinearSolver() { return _nonlinear_solver; } + void SetNonlinearSolver(std::unique_ptr nl_solver) + { + _nonlinear_solver = std::move(nl_solver); + } + + inline mfem::Solver * JacobianPreconditioner() const { return _jacobian_preconditioner.get(); } + inline mfem::Solver * JacobianSolver() const { return _jacobian_solver.get(); } + inline mfem::NewtonSolver * NonlinearSolver() const { return _nonlinear_solver.get(); } protected: /// Use of protected constructor to only allow construction by derived classes. @@ -70,13 +80,13 @@ class ProblemOperatorBase mfem::BlockVector _true_x, _true_rhs; /// Store the Jacobian preconditioner. - std::shared_ptr _jacobian_preconditioner{nullptr}; + std::unique_ptr _jacobian_preconditioner{nullptr}; /// Store the Jacobian solver. - std::shared_ptr _jacobian_solver{nullptr}; + std::unique_ptr _jacobian_solver{nullptr}; /// Store the non-linear solver. - std::shared_ptr _nonlinear_solver{nullptr}; + std::unique_ptr _nonlinear_solver{nullptr}; private: /// Update the block vectors and offsets after a mesh change. 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 655ed0d1e..f8426d7a5 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -55,16 +55,16 @@ TimeDomainEquationSystemProblemOperator::Update() // Rebuild the jacobian preconditioner. auto first_pfespace = GetEquationSystem()->_test_pfespaces.at(0); - auto precond = std::make_shared(first_pfespace); + auto precond = std::make_unique(first_pfespace); precond->SetSingularProblem(); precond->SetPrintLevel(-1); - _jacobian_preconditioner = precond; + _jacobian_preconditioner = std::move(precond); // Set new preconditioner. - std::static_pointer_cast(_jacobian_solver) - ->SetPreconditioner(*std::static_pointer_cast(_jacobian_preconditioner)); + static_cast(JacobianSolver()) + ->SetPreconditioner(*static_cast(JacobianPreconditioner())); // Set Jacobian matrix. auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); From af947c9d5fd38a7376ab3ff3ac3ca0cde2ac5f43 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Mon, 20 May 2024 16:12:26 +0000 Subject: [PATCH 14/20] Added accessor for JacobianPreconditioner; added documentation. --- src/problem_builders/problem_builder_base.cpp | 3 +-- src/problem_operators/problem_operator_base.hpp | 12 +++++++++--- .../time_domain_equation_system_problem_operator.cpp | 4 ++-- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index ba3a0eca3..17bb873d3 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -258,8 +258,7 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams solver_options.GetOptionalParam("PrintLevel", default_params._print_level); const auto k_dim = solver_options.GetOptionalParam("KDim", default_params._k_dim); - auto preconditioner = - dynamic_cast(GetProblem()->GetOperator()->JacobianPreconditioner()); + auto preconditioner = GetProblem()->GetOperator()->JacobianPreconditioner(); switch (type) { diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp index 5f5269a96..df13a6b7f 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -19,24 +19,30 @@ 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) { _nonlinear_solver = std::move(nl_solver); } - inline mfem::Solver * JacobianPreconditioner() const { return _jacobian_preconditioner.get(); } - inline mfem::Solver * JacobianSolver() const { return _jacobian_solver.get(); } - inline mfem::NewtonSolver * NonlinearSolver() const { return _nonlinear_solver.get(); } + /// Accessor for Jacobian preconditioner. + template + TSolver * JacobianPreconditioner() const + { + return static_cast(_jacobian_preconditioner.get()); + } protected: /// Use of protected constructor to only allow construction by derived classes. 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 f8426d7a5..ccc210c42 100644 --- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp +++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp @@ -63,8 +63,8 @@ TimeDomainEquationSystemProblemOperator::Update() _jacobian_preconditioner = std::move(precond); // Set new preconditioner. - static_cast(JacobianSolver()) - ->SetPreconditioner(*static_cast(JacobianPreconditioner())); + static_cast(_jacobian_solver.get()) + ->SetPreconditioner(*static_cast(_jacobian_preconditioner.get())); // Set Jacobian matrix. auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As(); From c1e60ea75c5660535988bdffcc11977d3f71dba6 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 09:46:00 +0000 Subject: [PATCH 15/20] Moved ODE solver into TimeDomainProblemOperator. --- src/executioners/transient_executioner.cpp | 2 +- src/problem_builders/problem_builder_base.cpp | 1 - src/problem_builders/problem_builder_base.hpp | 1 - .../time_domain_problem_builder.cpp | 17 ++++------------- .../time_domain_problem_builder.hpp | 2 -- .../time_domain_problem_operator.cpp | 13 +++++++++++++ .../time_domain_problem_operator.hpp | 15 +++++++++++++++ 7 files changed, 33 insertions(+), 18 deletions(-) diff --git a/src/executioners/transient_executioner.cpp b/src/executioners/transient_executioner.cpp index ddd7fa067..7ce7a99de 100644 --- a/src/executioners/transient_executioner.cpp +++ b/src/executioners/transient_executioner.cpp @@ -27,7 +27,7 @@ TransientExecutioner::Step(double dt, int it) const // Advance time step. _problem->_preprocessors.Solve(_t); - _problem->_ode_solver->Step(*(_problem->_f), _t, dt); + _problem->GetOperator()->ODESolver()->Step(*(_problem->_f), _t, dt); _problem->_postprocessors.Solve(_t); // Output data diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 17bb873d3..3a159c42c 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -7,7 +7,6 @@ Problem::~Problem() { // Ensure that all owned memory is properly freed! _f.reset(); - _ode_solver.reset(); } void diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 7aa800551..587dcc57d 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -30,7 +30,6 @@ class Problem hephaestus::Outputs _outputs; hephaestus::InputParameters _solver_options; - std::unique_ptr _ode_solver{nullptr}; std::unique_ptr _f{nullptr}; hephaestus::FECollections _fecs; diff --git a/src/problem_builders/time_domain_problem_builder.cpp b/src/problem_builders/time_domain_problem_builder.cpp index b1a62099e..916857a8f 100644 --- a/src/problem_builders/time_domain_problem_builder.cpp +++ b/src/problem_builders/time_domain_problem_builder.cpp @@ -3,17 +3,6 @@ namespace hephaestus { -void -TimeDomainProblem::Update() -{ - // 1. Call superclass. Updates the offsets and block vector. - Problem::Update(); - - // 2. The dimensions of the problem operator have now changed. We must call - // the ode solver's Init method to resize its internal vector. - _ode_solver->Init(*GetOperator()); -} - std::vector TimeDomainProblemBuilder::RegisterTimeDerivatives(std::vector gridfunction_names, hephaestus::GridFunctions & gridfunctions) @@ -59,8 +48,10 @@ TimeDomainProblemBuilder::InitializeOperator() void TimeDomainProblemBuilder::ConstructTimestepper() { - GetProblem()->_ode_solver = std::make_unique(); - GetProblem()->_ode_solver->Init(*(GetProblem()->GetOperator())); + 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 50c5e3caf..231b20586 100644 --- a/src/problem_builders/time_domain_problem_builder.hpp +++ b/src/problem_builders/time_domain_problem_builder.hpp @@ -31,8 +31,6 @@ class TimeDomainProblem : public Problem _problem_operator = std::make_unique(*this); } - void Update() override; - private: std::unique_ptr _problem_operator{nullptr}; }; diff --git a/src/problem_operators/time_domain_problem_operator.cpp b/src/problem_operators/time_domain_problem_operator.cpp index 0c02bcb36..a4ffa7211 100644 --- a/src/problem_operators/time_domain_problem_operator.cpp +++ b/src/problem_operators/time_domain_problem_operator.cpp @@ -20,4 +20,17 @@ GetTimeDerivativeNames(std::vector gridfunction_names) return time_derivative_names; } +void +TimeDomainProblemOperator::Update() +{ + ProblemOperatorBase::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); + } +} + } // 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 78804d313..d1ecbc4f8 100644 --- a/src/problem_operators/time_domain_problem_operator.hpp +++ b/src/problem_operators/time_domain_problem_operator.hpp @@ -21,9 +21,24 @@ class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public Pro 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); + } + + /// Returns a pointer to the ODE solver. + mfem::ODESolver * ODESolver() const { return _ode_solver.get(); } + + void Update() override; + protected: int & Width() final { return mfem::TimeDependentOperator::width; } int & Height() final { return mfem::TimeDependentOperator::height; } + +private: + /// Store the ODE solver. + std::unique_ptr _ode_solver{nullptr}; }; } // namespace hephaestus \ No newline at end of file From 63d0c6a9758890bd442c21ef1852b3046779c9d2 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 10:02:45 +0000 Subject: [PATCH 16/20] Block vector is now stored in ProblemOperatorBase; Added wrapper Step and Solve methods to problem operators which use stored ode solver. --- src/executioners/steady_executioner.cpp | 2 +- src/executioners/transient_executioner.cpp | 2 +- src/problem_builders/problem_builder_base.cpp | 14 -------------- src/problem_builders/problem_builder_base.hpp | 5 +---- src/problem_operators/problem_operator.hpp | 2 ++ src/problem_operators/problem_operator_base.cpp | 9 +++++++-- src/problem_operators/problem_operator_base.hpp | 5 ++++- .../time_domain_problem_operator.hpp | 6 ++++++ 8 files changed, 22 insertions(+), 23 deletions(-) diff --git a/src/executioners/steady_executioner.cpp b/src/executioners/steady_executioner.cpp index b5994d430..f8ac5d293 100644 --- a/src/executioners/steady_executioner.cpp +++ b/src/executioners/steady_executioner.cpp @@ -13,7 +13,7 @@ SteadyExecutioner::Solve() const { // Advance time step. _problem->_preprocessors.Solve(); - _problem->GetOperator()->Solve(*(_problem->_f)); + _problem->GetOperator()->Solve(); _problem->_postprocessors.Solve(); // Output data diff --git a/src/executioners/transient_executioner.cpp b/src/executioners/transient_executioner.cpp index 7ce7a99de..a947395e6 100644 --- a/src/executioners/transient_executioner.cpp +++ b/src/executioners/transient_executioner.cpp @@ -27,7 +27,7 @@ TransientExecutioner::Step(double dt, int it) const // Advance time step. _problem->_preprocessors.Solve(_t); - _problem->GetOperator()->ODESolver()->Step(*(_problem->_f), _t, dt); + _problem->GetOperator()->Step(_t, dt); _problem->_postprocessors.Solve(_t); // Output data diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 3a159c42c..50d794ccb 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -3,12 +3,6 @@ namespace hephaestus { -Problem::~Problem() -{ - // Ensure that all owned memory is properly freed! - _f.reset(); -} - void Problem::Update() { @@ -236,12 +230,6 @@ ProblemBuilder::ConstructJacobianSolver() ConstructJacobianSolverWithOptions(SolverType::HYPRE_GMRES); } -void -ProblemBuilder::ConstructBlockVector() -{ - GetProblem()->_f = std::make_unique(); -} - void ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams default_params) { @@ -388,8 +376,6 @@ ProblemBuilder::FinalizeProblem(bool build_operator) ConstructOperator(); } - ConstructBlockVector(); - InitializeAuxSolvers(); InitializeSources(); InitializeOperator(); diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 587dcc57d..cf0a50326 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -19,7 +19,7 @@ class Problem { public: Problem() = default; - virtual ~Problem(); + virtual ~Problem() = default; std::shared_ptr _pmesh{nullptr}; hephaestus::BCMap _bc_map; @@ -30,8 +30,6 @@ class Problem hephaestus::Outputs _outputs; hephaestus::InputParameters _solver_options; - std::unique_ptr _f{nullptr}; - hephaestus::FECollections _fecs; hephaestus::FESpaces _fespaces; hephaestus::GridFunctions _gridfunctions; @@ -98,7 +96,6 @@ class ProblemBuilder virtual void ConstructJacobianPreconditioner(); virtual void ConstructJacobianSolver(); virtual void ConstructNonlinearSolver(); - virtual void ConstructBlockVector(); virtual void ConstructOperator() = 0; virtual void ConstructTimestepper() = 0; diff --git a/src/problem_operators/problem_operator.hpp b/src/problem_operators/problem_operator.hpp index ace77b42e..12d151272 100644 --- a/src/problem_operators/problem_operator.hpp +++ b/src/problem_operators/problem_operator.hpp @@ -13,6 +13,8 @@ class ProblemOperator : public mfem::Operator, public ProblemOperatorBase ProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {} ~ProblemOperator() override = default; + virtual void Solve() { Solve(*_block_vector); } + virtual void Solve(mfem::Vector & X) {} void Mult(const mfem::Vector & x, mfem::Vector & y) const override {} diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp index ba3887bd7..55aeeb746 100644 --- a/src/problem_operators/problem_operator_base.cpp +++ b/src/problem_operators/problem_operator_base.cpp @@ -3,6 +3,11 @@ namespace hephaestus { +ProblemOperatorBase::ProblemOperatorBase(hephaestus::Problem & problem) : _problem{problem} +{ + _block_vector = std::make_unique(); +} + void ProblemOperatorBase::SetTrialVariables() { @@ -79,7 +84,7 @@ ProblemOperatorBase::Init() SetTrialVariables(); UpdateOffsets(); - UpdateBlockVector(*_problem._f); + UpdateBlockVector(*_block_vector); } void @@ -90,7 +95,7 @@ ProblemOperatorBase::Update() // Recalculate the offsets from gridfunction trial variables. UpdateOffsets(); - UpdateBlockVector(*_problem._f); + UpdateBlockVector(*_block_vector); } } \ 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 df13a6b7f..323c373ab 100644 --- a/src/problem_operators/problem_operator_base.hpp +++ b/src/problem_operators/problem_operator_base.hpp @@ -48,7 +48,7 @@ class ProblemOperatorBase /// 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) : _problem(problem) {} + explicit ProblemOperatorBase(hephaestus::Problem & problem); /// Set trial variables names. Override in derived classes. virtual void SetTrialVariableNames() {} @@ -94,6 +94,9 @@ class ProblemOperatorBase /// Store the non-linear solver. std::unique_ptr _nonlinear_solver{nullptr}; + /// Store the block vector. + std::unique_ptr _block_vector{nullptr}; + private: /// Update the block vectors and offsets after a mesh change. void UpdateOffsets(); diff --git a/src/problem_operators/time_domain_problem_operator.hpp b/src/problem_operators/time_domain_problem_operator.hpp index d1ecbc4f8..4f2da3482 100644 --- a/src/problem_operators/time_domain_problem_operator.hpp +++ b/src/problem_operators/time_domain_problem_operator.hpp @@ -27,6 +27,12 @@ class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public Pro _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); + } + /// Returns a pointer to the ODE solver. mfem::ODESolver * ODESolver() const { return _ode_solver.get(); } From 99903156f307af2f7ae08a6d5c3563df2011f034 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 10:04:22 +0000 Subject: [PATCH 17/20] Removed ODE solver accessor from TimeDomainProblemOperator. --- src/problem_operators/time_domain_problem_operator.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/problem_operators/time_domain_problem_operator.hpp b/src/problem_operators/time_domain_problem_operator.hpp index 4f2da3482..d36f17365 100644 --- a/src/problem_operators/time_domain_problem_operator.hpp +++ b/src/problem_operators/time_domain_problem_operator.hpp @@ -33,9 +33,6 @@ class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public Pro _ode_solver->Step(*_block_vector, t, dt); } - /// Returns a pointer to the ODE solver. - mfem::ODESolver * ODESolver() const { return _ode_solver.get(); } - void Update() override; protected: From 58814f85c2a739171ee365b2b8515dd6b7f782a3 Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 10:16:30 +0000 Subject: [PATCH 18/20] Removes ConstructJacobianPreconditioner method; moves preconditioner setup into ConstructJacobianSolver. --- 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 - 8 files changed, 4 insertions(+), 28 deletions(-) diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 785634773..56a21370c 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -59,7 +59,7 @@ DualFormulation::DualFormulation(std::string alpha_coef_name, } void -DualFormulation::ConstructJacobianPreconditioner() +DualFormulation::ConstructJacobianSolver() { auto precond = std::make_unique(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); @@ -68,11 +68,7 @@ DualFormulation::ConstructJacobianPreconditioner() precond->SetPrintLevel(-1); GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); -} -void -DualFormulation::ConstructJacobianSolver() -{ ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000}); } diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp index 78dd23a08..5a09f8923 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 ConstructJacobianPreconditioner() override; - void ConstructJacobianSolver() override; void ConstructOperator() override; diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 4c4110256..bdb5acee0 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -67,7 +67,7 @@ HCurlFormulation::ConstructOperator() } void -HCurlFormulation::ConstructJacobianPreconditioner() +HCurlFormulation::ConstructJacobianSolver() { auto precond = std::make_unique(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); @@ -76,11 +76,7 @@ HCurlFormulation::ConstructJacobianPreconditioner() precond->SetPrintLevel(-1); GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); -} -void -HCurlFormulation::ConstructJacobianSolver() -{ ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG); } diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp index b1c4b8ba4..dc9e491a1 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 ConstructJacobianPreconditioner() override; - void ConstructJacobianSolver() override; void RegisterGridFunctions() override; diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 3d06d1fe8..25ddd21e5 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -38,7 +38,7 @@ StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string } void -StaticsFormulation::ConstructJacobianPreconditioner() +StaticsFormulation::ConstructJacobianSolver() { auto precond = std::make_unique( GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace()); @@ -47,11 +47,7 @@ StaticsFormulation::ConstructJacobianPreconditioner() precond->SetPrintLevel(-1); GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); -} -void -StaticsFormulation::ConstructJacobianSolver() -{ ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES, {._max_iteration = 100, ._k_dim = 10}); } diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp index cf45d6f8a..b0a901695 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 ConstructJacobianPreconditioner() override; - void ConstructJacobianSolver() override; void ConstructOperator() override; diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index 50d794ccb..d5e4ef8dd 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -216,17 +216,13 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr(); precond->SetPrintLevel(GetGlobalPrintLevel()); GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); -} -void -ProblemBuilder::ConstructJacobianSolver() -{ ConstructJacobianSolverWithOptions(SolverType::HYPRE_GMRES); } @@ -380,7 +376,6 @@ ProblemBuilder::FinalizeProblem(bool build_operator) InitializeSources(); InitializeOperator(); - ConstructJacobianPreconditioner(); ConstructJacobianSolver(); ConstructNonlinearSolver(); diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index cf0a50326..39d85ac62 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 ConstructJacobianPreconditioner(); virtual void ConstructJacobianSolver(); virtual void ConstructNonlinearSolver(); From 91275f421c8b96d35b3ff1f3c7003f8cef9edd7b Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 10:32:41 +0000 Subject: [PATCH 19/20] Removed ConstructJacobianSolverWithOptions. All derived ProblemBuilder classes will create a jacobian solver / preconditioner, set their options and then pass it to the problem's operator. --- .../complex_maxwell_formulation.cpp | 3 +- src/formulations/Dual/dual_formulation.cpp | 11 +- src/formulations/HCurl/hcurl_formulation.cpp | 11 +- .../Magnetostatic/statics_formulation.cpp | 12 ++- src/problem_builders/problem_builder_base.cpp | 102 ++---------------- src/problem_builders/problem_builder_base.hpp | 33 ------ 6 files changed, 39 insertions(+), 133 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp index 238ca0c36..7105c862c 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp @@ -27,7 +27,8 @@ ComplexMaxwellFormulation::ComplexMaxwellFormulation(std::string alpha_coef_name void ComplexMaxwellFormulation::ConstructJacobianSolver() { - ConstructJacobianSolverWithOptions(SolverType::SUPER_LU); + auto solver = std::make_unique(GetProblem()->_comm); + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); } void diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 56a21370c..9cab0391a 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -67,9 +67,16 @@ DualFormulation::ConstructJacobianSolver() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + auto solver = std::make_unique(GetProblem()->_comm); + + solver->SetTol(1e-16); + solver->SetAbsTol(1e-16); + solver->SetMaxIter(1000); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); - ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000}); + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); } void diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index bdb5acee0..45ba244c6 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -75,9 +75,16 @@ HCurlFormulation::ConstructJacobianSolver() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + auto solver = std::make_unique(GetProblem()->_comm); + + solver->SetTol(1e-16); + solver->SetAbsTol(1e-16); + solver->SetMaxIter(1000); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); - ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG); + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); } void diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 25ddd21e5..143536f24 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -46,10 +46,16 @@ StaticsFormulation::ConstructJacobianSolver() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + auto solver = std::make_unique(GetProblem()->_comm); + + solver->SetTol(1e-16); + solver->SetMaxIter(100); + solver->SetKDim(10); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); - ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES, - {._max_iteration = 100, ._k_dim = 10}); + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); } void diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index d5e4ef8dd..fba6de562 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -219,102 +219,20 @@ void ProblemBuilder::ConstructJacobianSolver() { auto precond = std::make_unique(); - 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(); + precond->SetPrintLevel(GetGlobalPrintLevel()); - solver->SetTol(tolerance); - solver->SetMaxIter(max_iter); - solver->SetPrintLevel(print_level); + auto solver = std::make_unique(GetProblem()->_comm); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - case SolverType::SUPER_LU: - { - auto solver = std::make_unique(GetProblem()->_comm); + solver->SetTol(1e-16); + solver->SetAbsTol(1e-16); + solver->SetMaxIter(1000); + solver->SetKDim(10); + solver->SetPrintLevel(GetGlobalPrintLevel()); + solver->SetPreconditioner(*precond); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); - break; - } - default: - { - MFEM_ABORT("Unsupported solver type specified."); - break; - } - } + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); + GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); } 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; From ace7f88f80a0176fbacfc402f3d61737920f5b6a Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Wed, 22 May 2024 12:39:47 +0000 Subject: [PATCH 20/20] Revert "Removed ConstructJacobianSolverWithOptions. All derived ProblemBuilder classes will create a jacobian solver / preconditioner, set their options and then pass it to the problem's operator." This reverts commit 91275f421c8b96d35b3ff1f3c7003f8cef9edd7b. --- .../complex_maxwell_formulation.cpp | 3 +- src/formulations/Dual/dual_formulation.cpp | 11 +- src/formulations/HCurl/hcurl_formulation.cpp | 11 +- .../Magnetostatic/statics_formulation.cpp | 12 +-- src/problem_builders/problem_builder_base.cpp | 102 ++++++++++++++++-- src/problem_builders/problem_builder_base.hpp | 33 ++++++ 6 files changed, 133 insertions(+), 39 deletions(-) diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp index 7105c862c..238ca0c36 100644 --- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp +++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp @@ -27,8 +27,7 @@ ComplexMaxwellFormulation::ComplexMaxwellFormulation(std::string alpha_coef_name void ComplexMaxwellFormulation::ConstructJacobianSolver() { - auto solver = std::make_unique(GetProblem()->_comm); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); + ConstructJacobianSolverWithOptions(SolverType::SUPER_LU); } void diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp index 9cab0391a..56a21370c 100644 --- a/src/formulations/Dual/dual_formulation.cpp +++ b/src/formulations/Dual/dual_formulation.cpp @@ -67,16 +67,9 @@ DualFormulation::ConstructJacobianSolver() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - auto solver = std::make_unique(GetProblem()->_comm); - - solver->SetTol(1e-16); - solver->SetAbsTol(1e-16); - solver->SetMaxIter(1000); - solver->SetPrintLevel(GetGlobalPrintLevel()); - solver->SetPreconditioner(*precond); - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); + + ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000}); } void diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp index 45ba244c6..bdb5acee0 100644 --- a/src/formulations/HCurl/hcurl_formulation.cpp +++ b/src/formulations/HCurl/hcurl_formulation.cpp @@ -75,16 +75,9 @@ HCurlFormulation::ConstructJacobianSolver() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - auto solver = std::make_unique(GetProblem()->_comm); - - solver->SetTol(1e-16); - solver->SetAbsTol(1e-16); - solver->SetMaxIter(1000); - solver->SetPrintLevel(GetGlobalPrintLevel()); - solver->SetPreconditioner(*precond); - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); + + ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG); } void diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp index 143536f24..25ddd21e5 100644 --- a/src/formulations/Magnetostatic/statics_formulation.cpp +++ b/src/formulations/Magnetostatic/statics_formulation.cpp @@ -46,16 +46,10 @@ StaticsFormulation::ConstructJacobianSolver() precond->SetSingularProblem(); precond->SetPrintLevel(-1); - auto solver = std::make_unique(GetProblem()->_comm); - - solver->SetTol(1e-16); - solver->SetMaxIter(100); - solver->SetKDim(10); - solver->SetPrintLevel(GetGlobalPrintLevel()); - solver->SetPreconditioner(*precond); - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); + + ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES, + {._max_iteration = 100, ._k_dim = 10}); } void diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp index fba6de562..d5e4ef8dd 100644 --- a/src/problem_builders/problem_builder_base.cpp +++ b/src/problem_builders/problem_builder_base.cpp @@ -219,20 +219,102 @@ void ProblemBuilder::ConstructJacobianSolver() { auto precond = std::make_unique(); - precond->SetPrintLevel(GetGlobalPrintLevel()); - auto solver = std::make_unique(GetProblem()->_comm); + GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - solver->SetTol(1e-16); - solver->SetAbsTol(1e-16); - solver->SetMaxIter(1000); - solver->SetKDim(10); - solver->SetPrintLevel(GetGlobalPrintLevel()); - solver->SetPreconditioner(*precond); + ConstructJacobianSolverWithOptions(SolverType::HYPRE_GMRES); +} - GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond)); - GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver)); +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; + } + } } void diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp index 64a4a1979..39d85ac62 100644 --- a/src/problem_builders/problem_builder_base.hpp +++ b/src/problem_builders/problem_builder_base.hpp @@ -115,6 +115,39 @@ 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;