Merge branch 'develop' into lroberts36/fix-forest-multigrid

CHANGELOG.md

@@ -11,11 +11,13 @@
 - [[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 1053]](https://github.com/parthenon-hpc-lab/parthenon/pull/1053) Set the correct root level on restart
 - [[PR 1024]](https://github.com/parthenon-hpc-lab/parthenon/pull/1024) Add features to history output
 - [[PR 1031]](https://github.com/parthenon-hpc-lab/parthenon/pull/1031) Fix bug in non-cell centered AMR
 
 ### Infrastructure (changes irrelevant to downstream codes)
 - [[PR 1035]](https://github.com/parthenon-hpc-lab/parthenon/pull/1035) Fix multigrid infrastructure to work with forest
+- [[PR 1048]](https://github.com/parthenon-hpc-lab/parthenon/pull/1048) Tiny fixes to custom coords logic
 - [[PR 1028]](https://github.com/parthenon-hpc-lab/parthenon/pull/1028) Internal reorganization of LogicalLocation files
 - [[PR 1009]](https://github.com/parthenon-hpc-lab/parthenon/pull/1009) Move from a single octree to a forest of octrees
 

src/mesh/mesh.cpp

@@ -567,6 +567,10 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr,
   // Load balancing flag and parameters
   RegisterLoadBalancing_(pin);
 
+  // Initialize the forest
+  forest = forest::Forest::HyperRectangular(mesh_size, block_size, mesh_bcs);
+  root_level = forest.root_level;
+
   // SMR / AMR
   if (adaptive) {
     // read from file or from input?  input for now.
@@ -584,27 +588,29 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr,
 
   InitUserMeshData(this, pin);
 
-  // Populate logical locations
+  // Populate legacy logical locations
   auto lx123 = mesh_info.lx123;
   auto locLevelGidLidCnghostGflag = mesh_info.level_gid_lid_cnghost_gflag;
   current_level = -1;
   for (int i = 0; i < nbtotal; i++) {
     loclist[i] = LogicalLocation(locLevelGidLidCnghostGflag[5 * i], lx123[3 * i],
                                  lx123[3 * i + 1], lx123[3 * i + 2]);
-
-    if (loclist[i].level() > current_level) {
-      current_level = loclist[i].level();
-    }
   }
 
   // rebuild the Block Tree
-  forest = forest::Forest::HyperRectangular(mesh_size, block_size, mesh_bcs);
+
   for (int i = 0; i < nbtotal; i++) {
     forest.AddMeshBlock(forest.GetForestLocationFromLegacyTreeLocation(loclist[i]),
                         false);
   }
 
+  // Update the location list and levels to agree with forest levels
   loclist = forest.GetMeshBlockListAndResolveGids();
+
+  current_level = std::numeric_limits<int>::min();
+  for (const auto &loc : loclist)
+    current_level = std::max(current_level, loc.level());
+
   int nnb = loclist.size();
   if (nnb != nbtotal) {
     msg << "### FATAL ERROR in Mesh constructor" << std::endl
@@ -1278,4 +1284,24 @@ void Mesh::SetupMPIComms() {
 #endif
 }
 
+// Return list of locations and levels for the legacy tree
+// TODO(LFR): It doesn't make sense to offset the level by the
+//   legacy tree root level since the location indices are defined
+//   for loc.level(). It seems this level offset is required for
+//   the output to agree with the legacy output though.
+std::pair<std::vector<std::int64_t>, std::vector<std::int64_t>>
+Mesh::GetLevelsAndLogicalLocationsFlat() const noexcept {
+  std::vector<std::int64_t> levels, logicalLocations;
+  levels.reserve(nbtotal);
+  logicalLocations.reserve(nbtotal * 3);
+  for (auto loc : loclist) {
+    loc = forest.GetLegacyTreeLocation(loc);
+    levels.push_back(loc.level() - GetLegacyTreeRootLevel());
+    logicalLocations.push_back(loc.lx1());
+    logicalLocations.push_back(loc.lx2());
+    logicalLocations.push_back(loc.lx3());
+  }
+  return std::make_pair(levels, logicalLocations);
+}
+
 } // namespace parthenon

src/mesh/mesh.hpp

@@ -198,20 +198,8 @@ class Mesh {
   std::vector<int> GetNbList() const noexcept { return nblist; }
   std::vector<LogicalLocation> GetLocList() const noexcept { return loclist; }
 
-  // TODO(JMM): Put in implementation file?
-  auto GetLevelsAndLogicalLocationsFlat() const noexcept {
-    std::vector<std::int64_t> levels, logicalLocations;
-    levels.reserve(nbtotal);
-    logicalLocations.reserve(nbtotal * 3);
-    for (auto loc : loclist) {
-      loc = forest.GetLegacyTreeLocation(loc);
-      levels.push_back(loc.level() - GetLegacyTreeRootLevel());
-      logicalLocations.push_back(loc.lx1());
-      logicalLocations.push_back(loc.lx2());
-      logicalLocations.push_back(loc.lx3());
-    }
-    return std::make_pair(levels, logicalLocations);
-  }
+  std::pair<std::vector<std::int64_t>, std::vector<std::int64_t>>
+  GetLevelsAndLogicalLocationsFlat() const noexcept;
 
   void OutputMeshStructure(const int dim, const bool dump_mesh_structure = true);
 

src/outputs/parthenon_hdf5.cpp

@@ -244,6 +244,9 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
     const VariableVector<Real> &var_vec = pmb->meshblock_data.Get()->GetVariableVector();
     VariableVector<Real> coords_vars =
         GetAnyVariables(var_vec, {parthenon::Metadata::CoordinatesVec});
+    PARTHENON_DEBUG_REQUIRE(coords_vars.size() <= 1,
+                            "There can be at most one coordinates vector");
+
     VariableVector<Real> out;
     if (restart_) {
       // get all vars with flag Independent OR restart

src/outputs/parthenon_xdmf.cpp

@@ -133,7 +133,7 @@ void genXDMF(std::string hdfFile, Mesh *pm, SimTime *tm, IndexDomain domain, int
       mesh_type = "2DSMesh";
       dimstring = StringPrintf("%d %d", nx2 + n2_offset, nx1 + 1);
     } else {
-      PARTHENON_FAIL("1D curvilinear meshes not supported.");
+      PARTHENON_FAIL("Custom coordinates not supported in XDMF for 1D meshes.");
     }
   } else {
     mesh_type = "3DRectMesh";