Respond to Philipp's comments

parthenon-hpc-lab · Mar 19, 2024 · dddd0a2 · dddd0a2
1 parent 77057f4
commit dddd0a2
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 8 deletions.
diff --git a/example/poisson_gmg/main.cpp b/example/poisson_gmg/main.cpp
@@ -51,6 +51,8 @@ int main(int argc, char *argv[]) {
     if (driver_status != parthenon::DriverStatus::complete ||
         driver.final_rms_residual > 1.e-10 || driver.final_rms_error > 1.e-12)
       success = false;
+    if (driver.final_rms_residual != driver.final_rms_residual) success = false;
+    if (driver.final_rms_error != driver.final_rms_error) success = false;
   }
   // call MPI_Finalize and Kokkos::finalize if necessary
   pman.ParthenonFinalize();

diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
@@ -111,7 +111,8 @@ class MGSolver {
         &iter_counter);
     auto mg_finest = AddLinearOperatorTasks(itl, none, partition, pmesh);
     auto &md = pmesh->mesh_data.GetOrAdd("base", partition);
-    auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(mg_finest, itl, md, true);
+    auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(mg_finest, itl, md,
+                                                            pmesh->multilevel);
     auto calc_pointwise_res = eqs_.template Ax<u, res_err>(itl, comm, md);
     calc_pointwise_res = itl.AddTask(
         calc_pointwise_res, AddFieldsAndStoreInteriorSelect<rhs, res_err, res_err>, md,
@@ -152,7 +153,7 @@ class MGSolver {
   TaskID AddSetupTasks(TL_t &tl, TaskID dependence, int partition, Mesh *pmesh) {
     using namespace utils;
 
-    int min_level = 0;
+    int min_level = std::max(pmesh->GetGMGMaxLevel() - params_.max_coarsenings, 0);
     int max_level = pmesh->GetGMGMaxLevel();
 
     return AddMultiGridSetupPartitionLevel(tl, dependence, partition, max_level,
@@ -237,8 +238,9 @@ class MGSolver {
               Real *prhs = &pack(b, te, rhs_t(c), k, jb.s, ib.s);
               Real *xo = &pack(b, te, xold_t(c), k, jb.s, ib.s);
               Real *xn = &pack(b, te, xnew_t(c), k, jb.s, ib.s);
-              const int npoints =
-                  (jb.e - jb.s + 1) * (ib.e - ib.s + 1 + 2 * nghost) - 2 * nghost;
+              // Use ptr arithmetic to get the number of points we need to go over
+              // (including ghost zones) to get from (k, jb.s, ib.s) to (k, jb.e, ib.e)
+              const int npoints = &pack(b, te, Axold_t(c), k, jb.e, ib.e) - Ax + 1;
               parthenon::par_for_inner(
                   DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints - 1, [&](const int idx) {
                     const Real off_diag = Ax[idx] - diag[idx] * xo[idx];

diff --git a/src/solvers/solver_utils.hpp b/src/solvers/solver_utils.hpp
@@ -160,16 +160,18 @@ TaskStatus CopyData(const std::shared_ptr<MeshData<Real>> &md) {
   auto pack = desc.GetPack(md.get(), only_fine_on_composite);
   const int scratch_size = 0;
   const int scratch_level = 0;
+  // Warning: This inner loop strategy only works because we are using IndexDomain::entire
+  const int npoints_inner = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
   parthenon::par_for_outer(
       DEFAULT_OUTER_LOOP_PATTERN, "CopyData", DevExecSpace(), scratch_size, scratch_level,
       0, pack.GetNBlocks() - 1, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
         const int nvars =
             pack.GetUpperBound(b, in_t()) - pack.GetLowerBound(b, in_t()) + 1;
-        const int npoints = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
         for (int c = 0; c < nvars; ++c) {
           Real *in = &pack(b, te, in_t(c), kb.s, jb.s, ib.s);
           Real *out = &pack(b, te, out_t(c), kb.s, jb.s, ib.s);
-          parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints - 1,
+          parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0,
+                                   npoints_inner - 1,
                                    [&](const int idx) { out[idx] = in[idx]; });
         }
       });
@@ -198,18 +200,19 @@ TaskStatus AddFieldsAndStoreInteriorSelect(const std::shared_ptr<MeshData<Real>>
   auto pack = desc.GetPack(md.get(), include_block, only_fine_on_composite);
   const int scratch_size = 0;
   const int scratch_level = 0;
+  // Warning: This inner loop strategy only works because we are using IndexDomain::entire
+  const int npoints_inner = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
   parthenon::par_for_outer(
       DEFAULT_OUTER_LOOP_PATTERN, "AddFieldsAndStore", DevExecSpace(), scratch_size,
       scratch_level, 0, pack.GetNBlocks() - 1,
       KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
         const int nvars = pack.GetUpperBound(b, a_t()) - pack.GetLowerBound(b, a_t()) + 1;
-        const int npoints = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
         for (int c = 0; c < nvars; ++c) {
           Real *avar = &pack(b, te, a_t(c), kb.s, jb.s, ib.s);
           Real *bvar = &pack(b, te, b_t(c), kb.s, jb.s, ib.s);
           Real *out = &pack(b, te, out_t(c), kb.s, jb.s, ib.s);
           parthenon::par_for_inner(
-              DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints - 1,
+              DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints_inner - 1,
               [&](const int idx) { out[idx] = wa * avar[idx] + wb * bvar[idx]; });
         }
       });