WIP: Multi-D limited prolongation

src/CMakeLists.txt

@@ -12,6 +12,7 @@ add_executable(
         hydro/hydro_driver.cpp
         hydro/hydro.cpp
         hydro/glmmhd/dedner_source.cpp
+        hydro/prolongation/custom_ops.hpp
         hydro/srcterms/gravitational_field.hpp
         hydro/srcterms/tabular_cooling.hpp
         hydro/srcterms/tabular_cooling.cpp

src/hydro/hydro.cpp

@@ -30,6 +30,7 @@
 #include "glmmhd/glmmhd.hpp"
 #include "hydro.hpp"
 #include "outputs/outputs.hpp"
+#include "prolongation/custom_ops.hpp"
 #include "rsolvers/rsolvers.hpp"
 #include "srcterms/tabular_cooling.hpp"
 #include "utils/error_checking.hpp"
@@ -541,6 +542,8 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   Metadata m(
       {Metadata::Cell, Metadata::Independent, Metadata::FillGhost, Metadata::WithFluxes},
       std::vector<int>({nhydro + nscalars}), cons_labels);
+  m.RegisterRefinementOps<refinement_ops::ProlongateCellMinModMultiD,
+                          parthenon::refinement_ops::RestrictCellAverage>();
   pkg->AddField("cons", m);
 
   m = Metadata({Metadata::Cell, Metadata::Derived}, std::vector<int>({nhydro + nscalars}),

src/hydro/hydro_driver.cpp

@@ -12,7 +12,7 @@
 
 // Parthenon headers
 #include "bvals/cc/bvals_cc_in_one.hpp"
-#include "mesh/refinement_cc_in_one.hpp"
+#include "mesh/refinement_in_one.hpp"
 #include "refinement/refinement.hpp"
 #include <parthenon/parthenon.hpp>
 // AthenaPK headers
@@ -275,8 +275,7 @@ TaskCollection HydroDriver::MakeTaskCollection(BlockList_t &blocks, int stage) {
         tl.AddTask(none, parthenon::cell_centered_bvars::ReceiveBoundBufs<any>, mu0);
     auto set = tl.AddTask(recv, parthenon::cell_centered_bvars::SetBounds<any>, mu0);
     if (pmesh->multilevel) {
-      tl.AddTask(set, parthenon::cell_centered_refinement::RestrictPhysicalBounds,
-                 mu0.get());
+      tl.AddTask(set, parthenon::cell_centered_bvars::RestrictGhostHalos, mu0, false);
     }
   }
 

src/hydro/prolongation/custom_ops.hpp

@@ -0,0 +1,146 @@
+//========================================================================================
+// AthenaPK - a performance portable block structured AMR astrophysical MHD code.
+// Copyright (c) 2022, Athena-Parthenon Collaboration. All rights reserved.
+// Licensed under the BSD 3-Clause License (the "LICENSE").
+//========================================================================================
+// Parthenon performance portable AMR framework
+// Copyright(C) 2020-2022 The Parthenon collaboration
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+// (C) (or copyright) 2021. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S. Department of Energy/National Nuclear
+// Security Administration. All rights in the program are reserved by Triad
+// National Security, LLC, and the U.S. Department of Energy/National Nuclear
+// Security Administration. The Government is granted for itself and others
+// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do
+// so.
+//========================================================================================
+
+#ifndef HYDRO_PROLONGATION_CUSTOM_OPS_HPP_
+#define HYDRO_PROLONGATION_CUSTOM_OPS_HPP_
+
+#include <algorithm>
+#include <cstring>
+
+#include "coordinates/coordinates.hpp" // for coordinates
+#include "kokkos_abstraction.hpp"      // ParArray
+#include "mesh/domain.hpp"             // for IndesShape
+#include "mesh/mesh_refinement_ops.hpp"
+
+namespace Hydro {
+namespace refinement_ops {
+
+using parthenon::Coordinates_t;
+using parthenon::ParArray6D;
+
+template <int DIM>
+struct ProlongateCellMinModMultiD {
+  KOKKOS_FORCEINLINE_FUNCTION static void
+  Do(const int l, const int m, const int n, const int k, const int j, const int i,
+     const IndexRange &ckb, const IndexRange &cjb, const IndexRange &cib,
+     const IndexRange &kb, const IndexRange &jb, const IndexRange &ib,
+     const Coordinates_t &coords, const Coordinates_t &coarse_coords,
+     const ParArray6D<Real> *pcoarse, const ParArray6D<Real> *pfine) {
+    using namespace parthenon::refinement_ops::util;
+    auto &coarse = *pcoarse;
+    auto &fine = *pfine;
+
+    const Real fc = coarse(l, m, n, k, j, i);
+
+    int fi;
+    Real dx1fm, dx1fp, dx1m, dx1p;
+    GetGridSpacings<1>(coords, coarse_coords, cib, ib, i, &fi, &dx1m, &dx1p, &dx1fm,
+                       &dx1fp);
+    Real gx1c = GradMinMod(fc, coarse(l, m, n, k, j, i - 1), coarse(l, m, n, k, j, i + 1),
+                           dx1m, dx1p);
+
+    int fj = jb.s; // overwritten as needed
+    Real dx2fm = 0;
+    Real dx2fp = 0;
+    Real gx2c = 0;
+    if constexpr (DIM > 1) {
+      Real dx2m, dx2p;
+      GetGridSpacings<2>(coords, coarse_coords, cjb, jb, j, &fj, &dx2m, &dx2p, &dx2fm,
+                         &dx2fp);
+      gx2c = GradMinMod(fc, coarse(l, m, n, k, j - 1, i), coarse(l, m, n, k, j + 1, i),
+                        dx2m, dx2p);
+    }
+    int fk = kb.s;
+    Real dx3fm = 0;
+    Real dx3fp = 0;
+    Real gx3c = 0;
+    if constexpr (DIM > 2) {
+      Real dx3m, dx3p;
+      GetGridSpacings<3>(coords, coarse_coords, ckb, kb, k, &fk, &dx3m, &dx3p, &dx3fm,
+                         &dx3fp);
+      gx3c = GradMinMod(fc, coarse(l, m, n, k - 1, j, i), coarse(l, m, n, k + 1, j, i),
+                        dx3m, dx3p);
+    }
+
+    // Max. expected total difference. (dx#fm/p are positive by construction)
+    Real dqmax = std::abs(gx1c) * std::max(dx1fm, dx1p);
+    int jlim = 0;
+    int klim = 0;
+    if constexpr (DIM > 1) {
+      dqmax += std::abs(gx2c) * std::max(dx2fm, dx2fp);
+      jlim = 1;
+    }
+    if constexpr (DIM > 2) {
+      dqmax += std::abs(gx3c) * std::max(dx3fm, dx3fp);
+      klim = 1;
+    }
+    // Min/max values of all coarse cells used here
+    Real qmin = fc;
+    Real qmax = fc;
+    for (int koff = -klim; koff <= klim; koff++) {
+      for (int joff = -jlim; joff <= jlim; joff++) {
+        for (int ioff = -1; ioff <= 1; ioff++) {
+          qmin = std::min(qmin, coarse(l, m, n, k + koff, j + joff, i + ioff));
+          qmax = std::max(qmax, coarse(l, m, n, k + koff, j + joff, i + ioff));
+        }
+      }
+    }
+
+    // Scaling factor to limit all slopes simultaneously
+    Real alpha = 1.0;
+    if (dqmax * alpha > (qmax - fc)) {
+      alpha = (qmax - fc) / dqmax;
+    }
+    if (dqmax * alpha > (fc - qmin)) {
+      alpha = (fc - qmin) / dqmax;
+    }
+
+    // Ensure no new extrema are introduced in multi-D
+    gx1c *= alpha;
+    gx2c *= alpha;
+    gx3c *= alpha;
+
+    // KGF: add the off-centered quantities first to preserve FP symmetry
+    // JMM: Extraneous quantities are zero
+    fine(l, m, n, fk, fj, fi) = fc - (gx1c * dx1fm + gx2c * dx2fm + gx3c * dx3fm);
+    fine(l, m, n, fk, fj, fi + 1) = fc + (gx1c * dx1fp - gx2c * dx2fm - gx3c * dx3fm);
+    if constexpr (DIM > 1) {
+      fine(l, m, n, fk, fj + 1, fi) = fc - (gx1c * dx1fm - gx2c * dx2fp + gx3c * dx3fm);
+      fine(l, m, n, fk, fj + 1, fi + 1) =
+          fc + (gx1c * dx1fp + gx2c * dx2fp - gx3c * dx3fm);
+    }
+    if constexpr (DIM > 2) {
+      fine(l, m, n, fk + 1, fj, fi) = fc - (gx1c * dx1fm + gx2c * dx2fm - gx3c * dx3fp);
+      fine(l, m, n, fk + 1, fj, fi + 1) =
+          fc + (gx1c * dx1fp - gx2c * dx2fm + gx3c * dx3fp);
+      fine(l, m, n, fk + 1, fj + 1, fi) =
+          fc - (gx1c * dx1fm - gx2c * dx2fp - gx3c * dx3fp);
+      fine(l, m, n, fk + 1, fj + 1, fi + 1) =
+          fc + (gx1c * dx1fp + gx2c * dx2fp + gx3c * dx3fp);
+    }
+  }
+};
+} // namespace refinement_ops
+} // namespace Hydro
+
+#endif // HYDRO_PROLONGATION_CUSTOM_OPS_HPP_

src/pgen/linear_wave_mhd.cpp

@@ -157,6 +157,8 @@ void InitUserMeshData(Mesh *mesh, ParameterInput *pin) {
 
   Eigensystem(d0, u0, v0, w0, h0, bx0, by0, bz0, xfact, yfact, ev, rem, lem);
 
+  std::cout << "rem=" << rem[0][wave_flag] << std::endl;
+
   // TODO(pgrete) see how to get access to the SimTime object outside the driver
   // if (pin->GetOrAddBoolean("problem/linear_wave", "test", false) && ncycle == 0) {
   if (pin->GetOrAddBoolean("problem/linear_wave", "test", false)) {
@@ -366,9 +368,6 @@ void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
 //========================================================================================
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
-  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
-  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
-  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
   // Initialize the magnetic fields.  Note wavevector, eigenvectors, and other variables
   // are set in InitUserMeshData
 
@@ -392,9 +391,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto &coords = pmb->coords;
 
   // Initialize components of the vector potential
-  for (int k = kb.s - 1; k <= kb.e + 1; k++) {
-    for (int j = jb.s - 1; j <= jb.e + 1; j++) {
-      for (int i = ib.s - 1; i <= ib.e + 1; i++) {
+  IndexRange ib = pmb->cellbounds.GetBoundsI(IndexDomain::entire);
+  IndexRange jb = pmb->cellbounds.GetBoundsJ(IndexDomain::entire);
+  IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire);
+  for (int k = kb.s; k <= kb.e; k++) {
+    for (int j = jb.s; j <= jb.e; j++) {
+      for (int i = ib.s; i <= ib.e; i++) {
         a1(k, j, i) = A1(coords.x1v(i), coords.x2v(j), coords.x3v(k));
         a2(k, j, i) = A2(coords.x1v(i), coords.x2v(j), coords.x3v(k));
         a3(k, j, i) = A3(coords.x1v(i), coords.x2v(j), coords.x3v(k));
@@ -407,6 +409,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto &u_dev = rc->Get("cons").data;
   // initializing on host
   auto u = u_dev.GetHostMirrorAndCopy();
+  ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
+  jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
+  kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
+
+  const bool three_d = pmb->cellbounds.ncellsk(IndexDomain::entire) > 1;
   for (int k = kb.s; k <= kb.e; k++) {
     for (int j = jb.s; j <= jb.e; j++) {
       for (int i = ib.s; i <= ib.e; i++) {
@@ -422,10 +429,12 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
         u(IM2, k, j, i) = mx * cos_a2 * sin_a3 + my * cos_a3 - mz * sin_a2 * sin_a3;
         u(IM3, k, j, i) = mx * sin_a2 + mz * cos_a2;
 
-        u(IB1, k, j, i) = (a3(k, j + 1, i) - a3(k, j - 1, i)) / coords.dx2v(j) / 2.0 -
-                          (a2(k + 1, j, i) - a2(k - 1, j, i)) / coords.dx3v(k) / 2.0;
-        u(IB2, k, j, i) = (a1(k + 1, j, i) - a1(k - 1, j, i)) / coords.dx3v(k) / 2.0 -
-                          (a3(k, j, i + 1) - a3(k, j, i - 1)) / coords.dx1v(i) / 2.0;
+        u(IB1, k, j, i) = (a3(k, j + 1, i) - a3(k, j - 1, i)) / coords.dx2v(j) / 2.0;
+        u(IB1, k, j, i) -=
+            three_d ? (a2(k + 1, j, i) - a2(k - 1, j, i)) / coords.dx3v(k) / 2.0 : 0.0;
+        u(IB2, k, j, i) = -(a3(k, j, i + 1) - a3(k, j, i - 1)) / coords.dx1v(i) / 2.0;
+        u(IB2, k, j, i) +=
+            three_d ? (a1(k + 1, j, i) - a1(k - 1, j, i)) / coords.dx3v(k) / 2.0 : 0.0;
         u(IB3, k, j, i) = (a2(k, j, i + 1) - a2(k, j, i - 1)) / coords.dx1v(i) / 2.0 -
                           (a1(k, j + 1, i) - a1(k, j - 1, i)) / coords.dx2v(j) / 2.0;
 

tst/regression/test_suites/mhd_convergence/mhd_convergence.py

@@ -98,6 +98,7 @@ def Prepare(self, parameters, step):
             mb_nx1 //= 2
 
         parameters.driver_cmd_line_args = [
+            "job/problem_id=linear_wave_mhd",
             "parthenon/mesh/nx1=%d" % (2 * res),
             "parthenon/meshblock/nx1=%d" % mb_nx1,
             "parthenon/mesh/nx2=%d" % res,
@@ -164,7 +165,7 @@ def Analyse(self, parameters):
         )
 
         # quick and dirty test
-        if data[47, 4] > 6.14e-12:
+        if data[47, 4] > 4.05e-10:
             print("QUICK AND DIRTY TEST FAILED")
             analyze_status = False
 
@@ -182,21 +183,18 @@ def Analyse(self, parameters):
                 ),
             )
 
-        plt.plot([32, 512], [7e-7, 7e-7 / (512 / 32)], "--", label="first order")
+        plt.plot([32, 256], [1.1e-6, 1.1e-6 / (256 / 32)], "--", label="first order")
         plt.plot(
-            [32, 512], [1.7e-7, 1.7e-7 / (512 / 32) ** 2], "--", label="second order"
+            [32, 256], [1.7e-7, 1.7e-7 / (256 / 32) ** 2], "--", label="second order"
         )
         plt.plot(
-            [32, 512], [3.7e-8, 3.7e-8 / (512 / 32) ** 2], "--", label="second order"
+            [32, 256], [3.1e-8, 3.1e-8 / (256 / 32) ** 2], "--", label="second order"
         )
         plt.plot(
-            [32, 512], [5.6e-8, 5.6e-8 / (512 / 32) ** 3], "--", label="third order"
-        )
-        plt.plot(
-            [32, 512], [3.6e-9, 3.6e-9 / (512 / 32) ** 3], "--", label="third order"
+            [32, 256], [8.3e-8, 8.3e-8 / (256 / 32) ** 3], "--", label="third order"
         )
 
-        plt.ylim(1e-12, 5e-6)
+        plt.ylim(3e-10, 2e-6)
 
         plt.legend(bbox_to_anchor=(1, 1), loc="upper left")
         plt.xscale("log")