Move toward fixing regression test

parthenon-hpc-lab · Oct 4, 2023 · 0392b56 · 0392b56
1 parent 00a1d6c
commit 0392b56
Show file tree

Hide file tree

Showing 5 changed files with 32 additions and 10 deletions.
diff --git a/example/poisson_gmg/poisson_driver.cpp b/example/poisson_gmg/poisson_driver.cpp
@@ -35,8 +35,24 @@ using namespace parthenon::driver::prelude;
 namespace poisson_example {
 
 parthenon::DriverStatus PoissonDriver::Execute() {
+  using namespace parthenon;
+  using namespace poisson_package;
+  auto pkg = pmesh->packages.Get("poisson_package");
+  auto solver = pkg->Param<std::string>("solver");
+  auto *mg_solver =
+      pkg->MutableParam<parthenon::solvers::MGSolver<u, rhs, PoissonEquation>>(
+          "MGsolver");
+  auto *bicgstab_solver =
+      pkg->MutableParam<parthenon::solvers::BiCGSTABSolver<u, rhs, PoissonEquation>>(
+          "MGBiCGSTABsolver");
+
   pouts->MakeOutputs(pmesh, pinput);
   ConstructAndExecuteTaskLists<>(this);
+  if (solver == "BiCGSTAB") {
+    final_rms_residual = bicgstab_solver->GetFinalResidual();
+  } else if (solver == "MG") {
+    final_rms_residual = mg_solver->GetFinalResidual();
+  }
   pouts->MakeOutputs(pmesh, pinput);
   return DriverStatus::complete;
 }

diff --git a/example/poisson_gmg/poisson_equation.hpp b/example/poisson_gmg/poisson_equation.hpp
@@ -31,10 +31,12 @@ namespace poisson_package {
 // to consistently deal with coarse fine boundaries on the grid.
 class PoissonEquation {
  public:
+  bool do_flux_cor = false;
+
   template <class x_t, class out_t, bool only_md_level = false, class TL_t>
   parthenon::TaskID Ax(TL_t &tl, parthenon::TaskID depends_on,
-                       std::shared_ptr<parthenon::MeshData<Real>> &md, bool only_interior,
-                       bool do_flux_cor = false) {
+                       std::shared_ptr<parthenon::MeshData<Real>> &md,
+                       bool only_interior) {
     auto flux_res = tl.AddTask(depends_on, CalculateFluxes<x_t, only_md_level>, md);
     if (do_flux_cor && !only_md_level) {
       auto start_flxcor =

diff --git a/example/poisson_gmg/poisson_package.cpp b/example/poisson_gmg/poisson_package.cpp
@@ -99,13 +99,16 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   bicgstab_params.max_iters = max_poisson_iterations;
   bicgstab_params.residual_tolerance = res_tol;
   bicgstab_params.precondition = precondition;
-  bicgstab_params.flux_correct = flux_correct;
 
-  parthenon::solvers::MGSolver<u, rhs, PoissonEquation> mg_solver(pkg.get(), mg_params);
+  PoissonEquation eq;
+  eq.do_flux_cor = flux_correct;
+
+  parthenon::solvers::MGSolver<u, rhs, PoissonEquation> mg_solver(pkg.get(), mg_params,
+                                                                  eq);
   pkg->AddParam<>("MGsolver", mg_solver, parthenon::Params::Mutability::Mutable);
 
   parthenon::solvers::BiCGSTABSolver<u, rhs, PoissonEquation> bicg_solver(
-      pkg.get(), bicgstab_params);
+      pkg.get(), bicgstab_params, eq);
   pkg->AddParam<>("MGBiCGSTABsolver", bicg_solver,
                   parthenon::Params::Mutability::Mutable);
 

diff --git a/src/solvers/bicgstab_solver.hpp b/src/solvers/bicgstab_solver.hpp
@@ -36,7 +36,6 @@ struct BiCGSTABParams {
   Real residual_tolerance = 1.e-12;
   Real restart_threshold = -1.0;
   bool precondition = true;
-  bool flux_correct = true;
 };
 
 template <class x, class rhs, class equations>
@@ -54,7 +53,7 @@ class BiCGSTABSolver {
   BiCGSTABSolver(StateDescriptor *pkg, BiCGSTABParams params_in,
                  equations eq_in = equations())
       : preconditioner(pkg, MGParams(), eq_in), params_(params_in), iter_counter(0),
-        eqs_() {
+        eqs_(eq_in) {
     using namespace refinement_ops;
     auto mu = Metadata({Metadata::Cell, Metadata::Independent, Metadata::FillGhost,
                         Metadata::WithFluxes, Metadata::GMGRestrict});
@@ -123,7 +122,7 @@ class BiCGSTABSolver {
 
     // 2. v <- A u
     auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(precon1, itl, md, true);
-    auto get_v = eqs_.template Ax<u, v>(itl, comm, md, false, params_.flux_correct);
+    auto get_v = eqs_.template Ax<u, v>(itl, comm, md, false);
 
     // 3. rhat0v <- (rhat0, v)
     auto get_rhat0v =
@@ -172,7 +171,7 @@ class BiCGSTABSolver {
 
     // 7. t <- A u
     auto pre_t_comm = AddBoundaryExchangeTasks<BoundaryType::any>(precon2, itl, md, true);
-    auto get_t = eqs_.template Ax<u, t>(itl, pre_t_comm, md, false, params_.flux_correct);
+    auto get_t = eqs_.template Ax<u, t>(itl, pre_t_comm, md, false);
 
     // 8. omega <- (t,s) / (t,t)
     auto get_ts = DotProduct<t, s>(get_t, region, itl, i, reg_dep_id, &ts, md);
@@ -262,7 +261,7 @@ class BiCGSTABSolver {
   int GetCurrentIterations() const { return iter_counter; }
 
   Real GetFinalResidual() const { return final_residual; }
-  Real GetFinalIterations() const { return final_iteration; }
+  int GetFinalIterations() const { return final_iteration; }
 
  protected:
   MGSolver<u, rhs, equations> preconditioner;

diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
@@ -126,6 +126,8 @@ class MGSolver {
 
   Real GetSquaredResidualSum() const { return residual.val; }
   int GetCurrentIterations() const { return iter_counter; }
+  Real GetFinalResidual() const { return final_residual; }
+  int GetFinalIterations() const { return final_iteration; }
 
  protected:
   MGParams params_;