Fix edge case where meshblocks can be improperly de-refined

CHANGELOG.md

@@ -32,6 +32,7 @@
 - [[PR 1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1248]](https://github.com/parthenon-hpc-lab/parthenon/pull/1248) Fix edge case regarding AMR de-refinement logic
 - [[PR 1240]](https://github.com/parthenon-hpc-lab/parthenon/pull/1240) Fix I/O for CellMemAligned variables
 - [[PR 1229]](https://github.com/parthenon-hpc-lab/parthenon/pull/1229) Ensure builds function on 32 bit architectures
 - [[PR 1230]](https://github.com/parthenon-hpc-lab/parthenon/pull/1230) Add missing using statements in curvilinear coordinates classes

doc/sphinx/src/outputs.rst

@@ -167,6 +167,19 @@ immediately prior to restart files being written with the optional
 callbacks (if provided) will be called in that order before restart files are
 written.
 
+Changing output cadence from command line when restarting from file
+--------------------------------------------------------------------
+
+When restarting from a restart file, you can change any parameter
+input argument. To change the output cadence in this way, however,
+special care is required, as the code stores a ``next_time`` and a
+``next_n`` for the next time to output and next iteration to
+output. Fortunately, these can be overwritten on the command line in
+the usual way. To set the next time to output when restarting, run
+with ``parthenon/output*/next_time=some_number``.
+Note that this is independent of updating ``parthenon/output*/dt=some_number``.
+In other words, to change the cadence and apply the change immediately both ``dt`` (or ``dn``) and ``next_time`` (or ``next_n`` need to be updated simultaneously.
+
 Postprocessing/native analysis
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+# ========================================================================================
+#  (C) (or copyright) 2025. 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.
+# ========================================================================================
+
+import sys
+import numpy as np
+import h5py
+from argparse import ArgumentParser
+
+
+def compute_asymmetry(f, varname):
+    "Computes the asymmetry of var with varname in hdf5 output file object f"
+    xlocs = f["Locations/x"][:]
+    ylocs = f["Locations/y"][:]
+    iylocs = -np.flip(f["Locations/y"][:], axis=1)
+    matches = np.zeros((ylocs.shape[0]), dtype=int)
+    for b in range(ylocs.shape[0]):
+        for bb in range(ylocs.shape[0]):
+            if np.all(np.abs(ylocs[b] - iylocs[bb]) <= 1e-10) and np.all(
+                np.abs(xlocs[b] - xlocs[bb]) <= 1e-10
+            ):
+                matches[b] = bb
+
+    var = f[varname][:]
+    var_diff = np.zeros_like(var)
+    for b in range(ylocs.shape[0]):
+        bb = matches[b]
+        if np.any(ylocs[b] >= 0):
+            for d in range(var_diff.shape[1]):
+                sign = -1 if (var.shape[1] == 3) and d == 1 else 1
+                var_diff[b, d] = var[b, d] - sign * np.flip(var[bb, d], axis=-2)
+                var_diff[bb, d] = var[bb, d] - sign * np.flip(var[b, d], axis=-2)
+
+    return var_diff
+
+
+parser = ArgumentParser(
+    prog="compute_asymmetry.py",
+    description="compute asymmetry in X2 of a field and save it to the output file. "
+    + "Assumes mesh is symmetric about 0 in X2. "
+    + "Only works for cell- and node-centered data.",
+)
+parser.add_argument("field", type=str, help="Variable to compute")
+parser.add_argument("files", type=str, nargs="+", help="Files to compute")
+
+if __name__ == "__main__":
+    args = parser.parse_args()
+    for i, fname in enumerate(args.files):
+        with h5py.File(fname, "a") as f:
+            print(f"Computing asymmetry for {args.field} in {fname}...")
+            var_diff = compute_asymmetry(f, args.field)
+            savename = f"{args.field}_asymmetry"
+            try:
+                f.create_dataset(savename, data=var_diff)
+            except ValueError:
+                f[savename][:] = var_diff

scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py

@@ -1,3 +1,4 @@
+#!/usr/bin/env python
 # =========================================================================================
 # (C) (or copyright) 2020-2025. Triad National Security, LLC. All rights reserved.
 #
@@ -295,7 +296,7 @@ def plot_dump(
     if swarmx is not None and swarmy is not None:
         p.scatter(swarmx, swarmy, s=particlesize, c=swarmcolor)
     if colorbar is not None:
-        plt.colorbar(pm, fraction=0.02, pad=0.04, ax=p)
+        plt.colorbar(pm, label=colorbar, fraction=0.02, pad=0.04, ax=p)
 
     fig.savefig(output_file, dpi=300)
     plt.close(fig=fig)

scripts/python/packages/parthenon_tools/parthenon_tools/report_asymmetry.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env/python
+# ========================================================================================
+#  (C) (or copyright) 2025. 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.
+# ========================================================================================
+
+import sys
+import numpy as np
+from argparse import ArgumentParser
+
+import os
+
+sys.path.insert(
+    0, "../external/parthenon/scripts/python/packages/parthenon_tools/parthenon_tools"
+)
+import h5py
+from compute_asymmetry import compute_asymmetry
+
+parser = ArgumentParser(
+    prog="report_asymmetry.py",
+    description="Report asymmetry in X2 of all fields and save it to the output file. "
+    + "Assumes mesh is symmetric about 0 in X2. "
+    + "Only works for cell- and node-centered data.",
+)
+parser.add_argument("files", type=str, nargs="+", help="Files to compute")
+
+dset_type = h5py._hl.dataset.Dataset
+if __name__ == "__main__":
+    args = parser.parse_args()
+    for fname in args.files:
+        with h5py.File(fname, "r") as f:
+            print(f"Computing asymmetry in {fname} for vars...")
+            for k, v in f.items():
+                # Try to exclude face- and edge-centered data
+                if (type(v) == dset_type) and ("bnd_flux" not in k):
+                    if len(v.shape) > 2:
+                        print(f"\t...{k}:")
+                        try:
+                            var_diff = compute_asymmetry(f, k)
+                            print(
+                                "\t\t{:14e} out of {:14e}".format(
+                                    np.max(np.abs(var_diff)), np.max(np.abs(f[k]))
+                                )
+                            )
+                        except ValueError:
+                            print("\t\tcorrupted!")

src/defs.hpp

@@ -28,7 +28,6 @@
 
 #include "basic_types.hpp"
 #include "config.hpp"
-#include "mesh/forest/logical_location.hpp"
 #include "parthenon_arrays.hpp"
 
 namespace parthenon {
@@ -168,6 +167,26 @@ inline BoundaryFace GetOuterBoundaryFace(CoordinateDirection dir) {
   return BoundaryFace::undef;
 }
 
+inline std::array<int, 3> GetOffsetsFromBoundaryFace(BoundaryFace face) {
+  switch (face) {
+  case BoundaryFace::inner_x1:
+    return {-1, 0, 0};
+  case BoundaryFace::outer_x1:
+    return {1, 0, 0};
+  case BoundaryFace::inner_x2:
+    return {0, -1, 0};
+  case BoundaryFace::outer_x2:
+    return {0, 1, 0};
+  case BoundaryFace::inner_x3:
+    return {0, 0, -1};
+  case BoundaryFace::outer_x3:
+    return {0, 0, 1};
+  default:
+    PARTHENON_FAIL("Asking for offsets for an invalid BoundaryFace.");
+  }
+  return {0, 0, 0};
+}
+
 //------------------
 // strongly typed / scoped enums (C++11):
 //------------------

src/mesh/forest/logical_location.hpp

@@ -33,6 +33,7 @@
 #include <vector>
 
 #include "basic_types.hpp"
+#include "defs.hpp"
 #include "utils/error_checking.hpp"
 #include "utils/morton_number.hpp"
 
@@ -112,6 +113,32 @@ class LogicalLocation { // aggregate and POD type
                            (lx3() << 1) + ox3);
   }
 
+  bool IsLowerLeftCornerOfParent() const {
+    return ((lx1() & 1LL) == 0LL) && ((lx2() & 1LL) == 0LL) && ((lx3() & 1LL) == 0LL);
+  }
+
+  // Get the location in the parent, i.e. the lower left corner of the block
+  // is (0, 0, 0) and the upper right corner of the block is (1, 1, 1)
+  std::array<int, 3> GetLocationInParent() const {
+    return {(lx1() & 1LL) == 1LL, (lx2() & 1LL) == 1LL, (lx3() & 1LL) == 1LL};
+  }
+
+  bool IsOnTreeBoundary(int ox1, int ox2, int ox3) const {
+    const int nup = (1 << std::max(level(), 0)) - 1;
+    const bool bound1 =
+        (ox1 == 0) || ((ox1 == -1) && (lx1() == 0)) || ((ox1 == 1) && (lx1() == nup));
+    const bool bound2 =
+        (ox2 == 0) || ((ox2 == -1) && (lx2() == 0)) || ((ox2 == 1) && (lx2() == nup));
+    const bool bound3 =
+        (ox3 == 0) || ((ox3 == -1) && (lx3() == 0)) || ((ox3 == 1) && (lx3() == nup));
+    return bound1 && bound2 && bound3;
+  }
+
+  bool IsOnTreeBoundary(BoundaryFace face) const {
+    const auto [ox1, ox2, ox3] = GetOffsetsFromBoundaryFace(face);
+    return IsOnTreeBoundary(ox1, ox2, ox3);
+  }
+
   // LFR: This returns the face offsets of fine-coarse neighbor blocks as defined in
   // Athena++, which are stored in the NeighborBlock struct. I believe that these are
   // currently only required for flux correction and can eventually be removed when flux

src/mesh/forest/tree.cpp

@@ -112,22 +112,21 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) {
 
   if (enforce_proper_nesting) {
     LogicalLocation parent = ref_loc.GetParent();
-    int ox1 = ref_loc.lx1() - (parent.lx1() << 1);
-    int ox2 = ref_loc.lx2() - (parent.lx2() << 1);
-    int ox3 = ref_loc.lx3() - (parent.lx3() << 1);
-
-    for (int k = 0; k < (ndim > 2 ? 2 : 1); ++k) {
-      for (int j = 0; j < (ndim > 1 ? 2 : 1); ++j) {
-        for (int i = 0; i < (ndim > 0 ? 2 : 1); ++i) {
-          LogicalLocation neigh = parent.GetSameLevelNeighbor(
-              i + ox1 - 1, j + ox2 - (ndim > 1), k + ox3 - (ndim > 2));
+    const auto [ox1, ox2, ox3] = ref_loc.GetLocationInParent();
+    // Iterate over the 2^ndim possible blocks on level ref_loc.level() - 1 that could
+    // abut the newly created blocks on level ref_loc.level() + 1 and refine them
+    // if necessary
+    for (int k = -(ndim > 2); k < 1; ++k) {
+      for (int j = -(ndim > 1); j < 1; ++j) {
+        for (int i = -1; i < 1; ++i) {
+          LogicalLocation neigh = parent.GetSameLevelNeighbor(i + ox1, j + ox2, k + ox3);
           // Need to communicate this refinement action to possible neighboring tree(s)
           // and trigger refinement there
           int n_idx =
               neigh.NeighborTreeIndex(); // Note that this can point you back to this tree
           for (auto &[neighbor_tree, lcoord_trans] : neighbors[n_idx]) {
             nadded += neighbor_tree->Refine(
-                lcoord_trans.Transform(neigh, neighbor_tree->GetId()));
+                lcoord_trans.Transform(neigh, neighbor_tree->GetId()), true);
           }
         }
       }
@@ -322,14 +321,15 @@ RegionSize Tree::GetBlockDomain(const LogicalLocation &loc) const {
 std::array<BoundaryFlag, BOUNDARY_NFACES>
 Tree::GetBlockBCs(const LogicalLocation &loc) const {
   PARTHENON_REQUIRE(loc.IsInTree(), "Probably there is a mistake...");
-  std::array<BoundaryFlag, BOUNDARY_NFACES> block_bcs = boundary_conditions;
-  const int nblock = 1 << std::max(loc.level(), 0);
-  if (loc.lx1() != 0) block_bcs[BoundaryFace::inner_x1] = BoundaryFlag::block;
-  if (loc.lx1() != nblock - 1) block_bcs[BoundaryFace::outer_x1] = BoundaryFlag::block;
-  if (loc.lx2() != 0) block_bcs[BoundaryFace::inner_x2] = BoundaryFlag::block;
-  if (loc.lx2() != nblock - 1) block_bcs[BoundaryFace::outer_x2] = BoundaryFlag::block;
-  if (loc.lx3() != 0) block_bcs[BoundaryFace::inner_x3] = BoundaryFlag::block;
-  if (loc.lx3() != nblock - 1) block_bcs[BoundaryFace::outer_x3] = BoundaryFlag::block;
+  std::array<BoundaryFlag, BOUNDARY_NFACES> block_bcs;
+  for (auto face :
+       {BoundaryFace::inner_x1, BoundaryFace::outer_x1, BoundaryFace::inner_x2,
+        BoundaryFace::outer_x2, BoundaryFace::inner_x3, BoundaryFace::outer_x3}) {
+    if (loc.IsOnTreeBoundary(face))
+      block_bcs[face] = boundary_conditions[face];
+    else
+      block_bcs[face] = BoundaryFlag::block;
+  }
   return block_bcs;
 }