Fix bug in static refinement

src/mesh/forest/forest.cpp

@@ -180,7 +180,7 @@ Forest Forest::HyperRectangular(RegionSize mesh_size, RegionSize block_size,
   // Sort trees by their logical location in the tree mesh
   Forest fout;
   fout.root_level = ref_level;
-  fout.forest_level = level;
+  *(fout.forest_level) = level;
   for (auto &[loc, p] : ll_map)
     fout.AddTree(p.second);
   return fout;

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,7 +192,10 @@ class Mesh {
 
   int GetRootLevel() const noexcept { return root_level; }
   int GetLegacyTreeRootLevel() const noexcept {
-    return forest.root_level + forest.forest_level;
+    PARTHENON_REQUIRE(
+        forest.forest_level,
+        "The level of the forest must be defined to reference the legacy tree.");
+    return forest.root_level + *(forest.forest_level);
   }
 
   int GetMaxLevel() const noexcept { return max_level; }
@@ -323,17 +326,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);