Merge pull request #1071 from parthenon-hpc-lab/lroberts36/fix-static…

…-refinement

Fix bug in static refinement

CHANGELOG.md

@@ -17,6 +17,7 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1071]](https://github.com/parthenon-hpc-lab/parthenon/pull/1070) Fix bug in static mesh refinement related to redefinition of Mesh::root_level
 - [[PR 1073]](https://github.com/parthenon-hpc-lab/parthenon/pull/1073) Fix bug in AMR and sparse restarts
 - [[PR 1070]](https://github.com/parthenon-hpc-lab/parthenon/pull/1070) Correctly exclude flux vars from searches by default
 - [[PR 1049]](https://github.com/parthenon-hpc-lab/parthenon/pull/1049) Catch task failures from threads

src/mesh/forest/forest.hpp

@@ -39,7 +39,7 @@ class Forest {
 
  public:
   int root_level;
-  int forest_level;
+  std::optional<int> forest_level{};
 
   void AddTree(const std::shared_ptr<Tree> &in) {
     if (trees.count(in->GetId())) {

src/mesh/mesh.cpp

@@ -287,7 +287,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
           ref_size.xmax(X3DIR) = mesh_size.xmax(X3DIR);
         }
         int ref_lev = pin->GetInteger(pib->block_name, "level");
-        int lrlev = ref_lev + root_level;
+        int lrlev = ref_lev + GetLegacyTreeRootLevel();
         // range check
         if (ref_lev < 1) {
           msg << "### FATAL ERROR in Mesh constructor" << std::endl
@@ -324,9 +324,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
         for (auto dir : {X1DIR, X2DIR, X3DIR}) {
           if (!mesh_size.symmetry(dir)) {
             l_region_min[dir - 1] =
-                GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir));
+                GetLegacyLLFromMeshCoordinate(dir, lrlev, ref_size.xmin(dir));
             l_region_max[dir - 1] =
-                GetLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir));
+                GetLegacyLLFromMeshCoordinate(dir, lrlev, ref_size.xmax(dir));
             l_region_min[dir - 1] =
                 std::max(l_region_min[dir - 1], static_cast<std::int64_t>(0));
             l_region_max[dir - 1] =
@@ -335,7 +335,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
             auto current_loc =
                 LogicalLocation(lrlev, l_region_max[0], l_region_max[1], l_region_max[2]);
             // Remove last block if it just it's boundary overlaps with the region
-            if (GetMeshCoordinate(dir, BlockLocation::Left, current_loc) ==
+            if (GetLegacyMeshCoordinate(dir, BlockLocation::Left, current_loc) ==
                 ref_size.xmax(dir))
               l_region_max[dir - 1]--;
             if (l_region_min[dir - 1] % 2 == 1) l_region_min[dir - 1]--;

src/mesh/mesh.hpp

@@ -192,8 +192,8 @@ class Mesh {
       PostStepUserDiagnosticsInLoop = PostStepUserDiagnosticsInLoopDefault;
 
   int GetRootLevel() const noexcept { return root_level; }
-  int GetLegacyTreeRootLevel() const noexcept {
-    return forest.root_level + forest.forest_level;
+  int GetLegacyTreeRootLevel() const {
+    return forest.root_level + forest.forest_level.value();
   }
 
   int GetMaxLevel() const noexcept { return max_level; }
@@ -324,17 +324,19 @@ class Mesh {
 
   // Transform from logical location coordinates to uniform mesh coordinates accounting
   // for root grid
-  Real GetMeshCoordinate(CoordinateDirection dir, BlockLocation bloc,
-                         const LogicalLocation &loc) const {
+  Real GetLegacyMeshCoordinate(CoordinateDirection dir, BlockLocation bloc,
+                               const LogicalLocation &loc) const {
     auto xll = loc.LLCoord(dir, bloc);
-    auto root_fac = static_cast<Real>(1 << root_level) / static_cast<Real>(nrbx[dir - 1]);
+    auto root_fac = static_cast<Real>(1 << GetLegacyTreeRootLevel()) /
+                    static_cast<Real>(nrbx[dir - 1]);
     xll *= root_fac;
     return mesh_size.xmin(dir) * (1.0 - xll) + mesh_size.xmax(dir) * xll;
   }
 
-  std::int64_t GetLLFromMeshCoordinate(CoordinateDirection dir, int level,
-                                       Real xmesh) const {
-    auto root_fac = static_cast<Real>(1 << root_level) / static_cast<Real>(nrbx[dir - 1]);
+  std::int64_t GetLegacyLLFromMeshCoordinate(CoordinateDirection dir, int level,
+                                             Real xmesh) const {
+    auto root_fac = static_cast<Real>(1 << GetLegacyTreeRootLevel()) /
+                    static_cast<Real>(nrbx[dir - 1]);
     auto xLL = (xmesh - mesh_size.xmin(dir)) /
                (mesh_size.xmax(dir) - mesh_size.xmin(dir)) / root_fac;
     return static_cast<std::int64_t>((1 << std::max(level, 0)) * xLL);

tst/regression/CMakeLists.txt

@@ -125,7 +125,8 @@ if (ENABLE_HDF5)
   list(APPEND TEST_DIRS sparse_advection)
   list(APPEND TEST_PROCS ${NUM_MPI_PROC_TESTING})
   list(APPEND TEST_ARGS "--driver ${PROJECT_BINARY_DIR}/example/sparse_advection/sparse_advection-example \
-    --driver_input ${CMAKE_CURRENT_SOURCE_DIR}/test_suites/sparse_advection/parthinput.sparse_advection")
+    --driver_input ${CMAKE_CURRENT_SOURCE_DIR}/test_suites/sparse_advection/parthinput.sparse_advection \
+    --num_steps 3")
   list(APPEND EXTRA_TEST_LABELS "")
 
   list(APPEND TEST_DIRS particle_tracers)

tst/regression/test_suites/sparse_advection/sparse_advection.py

@@ -33,6 +33,27 @@ def Prepare(self, parameters, step):
                 "parthenon/sparse/enable_sparse=false",
             ]
 
+        # Run a test with two trees
+        if step == 2:
+            parameters.driver_cmd_line_args = [
+                "parthenon/mesh/nx2=32",
+                "parthenon/job/problem_id=sparse_twotree",
+            ]
+
+        # Run a test with two trees and a statically refined region
+        if step == 3:
+            parameters.driver_cmd_line_args = [
+                "parthenon/mesh/nx2=32",
+                "parthenon/time/nlim=50",
+                "parthenon/job/problem_id=sparse_twotree_static",
+                "parthenon/mesh/refinement=static",
+                "parthenon/static_refinement0/x1min=-0.75",
+                "parthenon/static_refinement0/x1max=-0.5",
+                "parthenon/static_refinement0/x2min=-0.75",
+                "parthenon/static_refinement0/x2max=-0.5",
+                "parthenon/static_refinement0/level=3",
+            ]
+
         return parameters
 
     def Analyse(self, parameters):
@@ -77,5 +98,35 @@ def Analyse(self, parameters):
                 tol=1e-12,
                 check_metadata=False,
             )
+            if delta != 0:
+                print("Sparse advection failed for standard AMR grid setup.")
+                return False
+
+            delta = compare(
+                [
+                    "sparse_twotree.out0.final.phdf",
+                    parameters.parthenon_path
+                    + "/tst/regression/gold_standard/sparse_twotree.out0.final.phdf",
+                ],
+                one=True,
+                tol=1e-12,
+                check_metadata=False,
+            )
+            if delta != 0:
+                print("Sparse advection failed for two-tree AMR grid setup.")
+                return False
+
+            delta = compare(
+                [
+                    "sparse_twotree_static.out0.final.phdf",
+                    parameters.parthenon_path
+                    + "/tst/regression/gold_standard/sparse_twotree_static.out0.final.phdf",
+                ],
+                one=True,
+                tol=1e-12,
+                check_metadata=False,
+            )
+            if delta != 0:
+                print("Sparse advection failed for two-tree SMR grid setup.")
 
         return delta == 0