From e6fec8aee7bc92297ec493806f7324679735b5a0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:39:17 -0400 Subject: [PATCH 01/13] lukes nicer cleanup --- src/mesh/mesh-amr_loadbalance.cpp | 44 +++++++++++++------------------ 1 file changed, 19 insertions(+), 25 deletions(-) diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index 00d768720988..6852467ffc46 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -502,8 +502,8 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { // collect refinement flags from all the meshblocks // count the number of the blocks to be (de)refined - nref[Globals::my_rank] = 0; - nderef[Globals::my_rank] = 0; + nref[Globals::my_rank] = 0; // n to refine. Local per rank. + nderef[Globals::my_rank] = 0; // n to derefine. Local per rank. for (auto const &pmb : block_list) { if (pmb->pmr->refine_flag_ == 1) nref[Globals::my_rank]++; if (pmb->pmr->refine_flag_ == -1) nderef[Globals::my_rank]++; @@ -516,7 +516,7 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { #endif // count the number of the blocks to be (de)refined and displacement - int tnref = 0, tnderef = 0; + int tnref = 0, tnderef = 0; // totals to refine/de-refine. Not local per rank. for (int n = 0; n < Globals::nranks; n++) { tnref += nref[n]; tnderef += nderef[n]; @@ -542,11 +542,11 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { } // allocate memory for the location arrays - LogicalLocation *lref{}, *lderef{}, *clderef{}; - if (tnref > 0) lref = new LogicalLocation[tnref]; + std::vector lref, lderef, clderef; + if (tnref > 0) lref.resize(tnref); if (tnderef >= nleaf) { - lderef = new LogicalLocation[tnderef]; - clderef = new LogicalLocation[tnderef / nleaf]; + lderef.resize(tnderef); + clderef.resize(tnderef / nleaf); } // collect the locations and costs @@ -558,35 +558,33 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { #ifdef MPI_PARALLEL if (tnref > 0) { PARTHENON_MPI_CHECK(MPI_Allgatherv(MPI_IN_PLACE, bnref[Globals::my_rank], MPI_BYTE, - lref, bnref.data(), brdisp.data(), MPI_BYTE, + lref.data(), bnref.data(), brdisp.data(), MPI_BYTE, MPI_COMM_WORLD)); } if (tnderef >= nleaf) { PARTHENON_MPI_CHECK(MPI_Allgatherv(MPI_IN_PLACE, bnderef[Globals::my_rank], MPI_BYTE, - lderef, bnderef.data(), bddisp.data(), MPI_BYTE, - MPI_COMM_WORLD)); + lderef.data(), bnderef.data(), bddisp.data(), + MPI_BYTE, MPI_COMM_WORLD)); } #endif // calculate the list of the newly derefined blocks int ctnd = 0; if (tnderef >= nleaf) { - int lk = 0, lj = 0; - if (!mesh_size.symmetry(X2DIR)) lj = 1; - if (!mesh_size.symmetry(X3DIR)) lk = 1; + int lk = !mesh_size.symmetry(X3DIR); + int lj = !mesh_size.symmetry(X2DIR); for (int n = 0; n < tnderef; n++) { if ((lderef[n].lx1() & 1LL) == 0LL && (lderef[n].lx2() & 1LL) == 0LL && (lderef[n].lx3() & 1LL) == 0LL) { int r = n, rr = 0; - for (std::int64_t k = 0; k <= lk; k++) { - for (std::int64_t j = 0; j <= lj; j++) { - for (std::int64_t i = 0; i <= 1; i++) { + // offsets from "bottom left" corner of parent block + for (std::int64_t ok = 0; ok <= lk; ok++) { + for (std::int64_t oj = 0; oj <= lj; oj++) { + for (std::int64_t oi = 0; oi <= 1; oi++) { if (r < tnderef) { - if ((lderef[n].lx1() + i) == lderef[r].lx1() && - (lderef[n].lx2() + j) == lderef[r].lx2() && - (lderef[n].lx3() + k) == lderef[r].lx3() && - lderef[n].level() == lderef[r].level()) + if (lderef[n].GetSameLevelNeighbor(oi, oj, ok) == lderef[r]) { rr++; + } r++; } } @@ -601,12 +599,11 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { } // sort the lists by level if (ctnd > 1) { - std::sort(clderef, &(clderef[ctnd - 1]), + std::sort(&(clderef[0]), &(clderef[ctnd - 1]), [](const LogicalLocation &left, const LogicalLocation &right) { return left.level() > right.level(); }); } - if (tnderef >= nleaf) delete[] lderef; // Now the lists of the blocks to be refined and derefined are completed // Start tree manipulation @@ -614,14 +611,11 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { for (int n = 0; n < tnref; n++) { nnew += forest.Refine(lref[n]); } - if (tnref != 0) delete[] lref; // Step 2. perform derefinement for (int n = 0; n < ctnd; n++) { ndel += forest.Derefine(clderef[n]); } - - if (tnderef >= nleaf) delete[] clderef; } //---------------------------------------------------------------------------------------- From 51d5a4aa7b99f2a1099a57adb31c155d466956ee Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:39:38 -0400 Subject: [PATCH 02/13] add compute asymmetry script --- .../parthenon_tools/compute_asymmetry.py | 65 +++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py b/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py new file mode 100644 index 000000000000..7cafb439c920 --- /dev/null +++ b/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py @@ -0,0 +1,65 @@ +#!/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): + var_diff[b] = var[b] - np.flip(var[bb], axis=-2) + var_diff[bb] = var[bb] - np.flip(var[b], axis=-2) + + if np.any(np.abs(var_diff) > 1): + print("var_diff > 1e-8!", np.where(np.abs(var_diff) > 1)) + + 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.", +) +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 From 822698b7f131c1ff49d81af16e51727842951836 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:39:55 -0400 Subject: [PATCH 03/13] Add shebang to movie2d --- .../python/packages/parthenon_tools/parthenon_tools/movie2d.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py b/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py index 02509746604c..83c4be2ab2c1 100644 --- a/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py +++ b/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. # From d622d1c5b2472ceed8859f3c877b29429b19bd77 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:40:13 -0400 Subject: [PATCH 04/13] document how to change output cadence from restart --- doc/sphinx/src/outputs.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index 6162e49242d3..ace250b81411 100644 --- a/doc/sphinx/src/outputs.rst +++ b/doc/sphinx/src/outputs.rst @@ -167,6 +167,17 @@ 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 paremeter +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``. + Postprocessing/native analysis ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ From 88329de9f56b73166f3c33a0b7952b01e3ee9eeb Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:50:23 -0400 Subject: [PATCH 05/13] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d8f6ae2cb6ba..5665db843c4b 100644 --- a/CHANGELOG.md +++ b/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 From 628df9b83ad3f1c2b86371c9efbe924c5360c9f8 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:50:57 -0400 Subject: [PATCH 06/13] CC --- src/mesh/mesh-amr_loadbalance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index 6852467ffc46..ef6a0fda4ee9 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -1,6 +1,6 @@ //======================================================================================== // Parthenon performance portable AMR framework -// Copyright(C) 2020-2024 The Parthenon collaboration +// Copyright(C) 2020-2025 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== // Athena++ astrophysical MHD code From 6facd7ffc70889b5ff84b7b42071f9327d4032d9 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 24 Apr 2025 12:54:32 -0400 Subject: [PATCH 07/13] CC --- src/mesh/mesh-amr_loadbalance.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index ef6a0fda4ee9..c29328d88a59 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -7,7 +7,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-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 @@ -18,7 +18,7 @@ // 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. //======================================================================================== -//! \file mesh_amr.cpp +//! \file mesh-amr_loadbalance.cpp // \brief implementation of Mesh::AdaptiveMeshRefinement() and related utilities #include From 508f0c7ad1d2a4b46b2d4c8bd69d3cedacdee1e1 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 24 Apr 2025 13:29:48 -0600 Subject: [PATCH 08/13] Try to move all LogicalLocation related logic to class methods --- src/defs.hpp | 21 +++++++++++++- src/mesh/forest/logical_location.hpp | 27 +++++++++++++++++ src/mesh/forest/tree.cpp | 36 +++++++++++------------ src/mesh/mesh-amr_loadbalance.cpp | 43 +++++++++------------------- 4 files changed, 79 insertions(+), 48 deletions(-) diff --git a/src/defs.hpp b/src/defs.hpp index 552e9cb32087..0983438c294f 100644 --- a/src/defs.hpp +++ b/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 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}; + case default: + PARTHENON_FAIL("Asking for offsets for an invalid BoundaryFace."); + } + return {0, 0, 0}; +} + //------------------ // strongly typed / scoped enums (C++11): //------------------ diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 1c7a5d7e073b..10a4586ac0d1 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -33,6 +33,7 @@ #include #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 GetLocationInParent() const { + return {(lx1() & 1LL) == 1LL, (lx2() & 1LL) == 1LL, (lx3() & 1LL) == 1LL}; + } + + bool IsOnTreeBoundary(int ox1, int ox2, int ox3) { + const int nup = 1 << std::max(level(), 0) - 1; + const bool bound1 = + (ox1 == 0) || (ox1 == -1 && loc.lx1() == 0) || (ox1 == 1 && loc.lx1() == nup); + const bool bound2 = + (ox2 == 0) || (ox2 == -1 && loc.lx2() == 0) || (ox2 == 1 && loc.lx2() == nup); + const bool bound3 = + (ox3 == 0) || (ox3 == -1 && loc.lx3() == 0) || (ox3 == 1 && loc.lx3() == nup); + return bound1 && bound2 && bound3; + } + + bool IsOnTreeBoundary(BoundaryFace face) { + 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 diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index d1d019ee5bd4..f35248275efb 100644 --- a/src/mesh/forest/tree.cpp +++ b/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 Tree::GetBlockBCs(const LogicalLocation &loc) const { PARTHENON_REQUIRE(loc.IsInTree(), "Probably there is a mistake..."); - std::array 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 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; } diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index c29328d88a59..dddbfbf1b70b 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -54,21 +54,21 @@ namespace parthenon { // tag = local id of destination (remaining bits) + ox1(1 bit) + ox2(1 bit) + ox3(1 bit) // + physics(5 bits) -int CreateAMRMPITag(int lid, int ox1, int ox2, int ox3) { +int CreateAMRMPITag(int lid, std::optional loc = {}) { // the trailing zero is used as "id" to indicate an AMR related tag - return (lid << 8) | (ox1 << 7) | (ox2 << 6) | (ox3 << 5) | 0; + int tag = lid << 8; + if (loc) { + auto offsets = (*loc).GetLocationInParent(); + tag = tag | (offsets[0] << 7) | (offsets[1] << 6) | (offsets[2] << 5) | 0; + } + return tag; } MPI_Request SendCoarseToFine(int lid_recv, int dest_rank, const LogicalLocation &fine_loc, Variable *var, Mesh *pmesh) { MPI_Request req; MPI_Comm comm = pmesh->GetMPIComm(var->label()); - - const int ox1 = ((fine_loc.lx1() & 1LL) == 1LL); - const int ox2 = ((fine_loc.lx2() & 1LL) == 1LL); - const int ox3 = ((fine_loc.lx3() & 1LL) == 1LL); - - int tag = CreateAMRMPITag(lid_recv, ox1, ox2, ox3); + int tag = CreateAMRMPITag(lid_recv, fine_loc); if (var->IsAllocated()) { PARTHENON_MPI_CHECK(MPI_Isend(var->data.data(), var->data.size(), MPI_PARTHENON_REAL, dest_rank, tag, comm, &req)); @@ -83,14 +83,10 @@ MPI_Request SendCoarseToFine(int lid_recv, int dest_rank, const LogicalLocation bool TryRecvCoarseToFine(int lid_recv, int send_rank, const LogicalLocation &fine_loc, Variable *var_in, Variable *var, MeshBlock *pmb, Mesh *pmesh) { - const int ox1 = ((fine_loc.lx1() & 1LL) == 1LL); - const int ox2 = ((fine_loc.lx2() & 1LL) == 1LL); - const int ox3 = ((fine_loc.lx3() & 1LL) == 1LL); - int test = 1; #ifdef MPI_PARALLEL MPI_Comm comm = pmesh->GetMPIComm(var->label()); - int tag = CreateAMRMPITag(lid_recv, ox1, ox2, ox3); + int tag = CreateAMRMPITag(lid_recv, fine_loc); MPI_Status status; PARTHENON_MPI_CHECK(MPI_Iprobe(send_rank, tag, comm, &test, &status)); #endif @@ -154,12 +150,7 @@ MPI_Request SendFineToCoarse(int lid_recv, int dest_rank, const LogicalLocation Variable *var, Mesh *pmesh) { MPI_Request req; MPI_Comm comm = pmesh->GetMPIComm(var->label()); - - const int ox1 = ((fine_loc.lx1() & 1LL) == 1LL); - const int ox2 = ((fine_loc.lx2() & 1LL) == 1LL); - const int ox3 = ((fine_loc.lx3() & 1LL) == 1LL); - - int tag = CreateAMRMPITag(lid_recv, ox1, ox2, ox3); + int tag = CreateAMRMPITag(lid_recv, fine_loc); if (var->IsAllocated()) { PARTHENON_MPI_CHECK(MPI_Isend(var->coarse_s.data(), var->coarse_s.size(), MPI_PARTHENON_REAL, dest_rank, tag, comm, &req)); @@ -175,14 +166,10 @@ bool TryRecvFineToCoarse(int lid_recv, int send_rank, const LogicalLocation &fin Variable *var_in, Variable *var, MeshBlock *pmb, Mesh *pmesh) { const int ndim = pmb->pmy_mesh->ndim; - const int ox1 = ((fine_loc.lx1() & 1LL) == 1LL); - const int ox2 = ((fine_loc.lx2() & 1LL) == 1LL); - const int ox3 = ((fine_loc.lx3() & 1LL) == 1LL); - int test = 1; #ifdef MPI_PARALLEL MPI_Comm comm = pmesh->GetMPIComm(var->label()); - int tag = CreateAMRMPITag(lid_recv, ox1, ox2, ox3); + int tag = CreateAMRMPITag(lid_recv, fine_loc); MPI_Status status; PARTHENON_MPI_CHECK(MPI_Iprobe(send_rank, tag, comm, &test, &status)); #endif @@ -252,7 +239,7 @@ MPI_Request SendSameToSame(int lid_recv, int dest_rank, Variable *var, MeshBlock *pmb, Mesh *pmesh) { MPI_Request req; MPI_Comm comm = pmesh->GetMPIComm(var->label()); - int tag = CreateAMRMPITag(lid_recv, 0, 0, 0); + int tag = CreateAMRMPITag(lid_recv); if (var->IsAllocated()) { // Metadata about this field also needs to be copied from this rank to the // receiving rank (namely the dereference count and the dealloc_count). Not @@ -286,8 +273,7 @@ MPI_Request SendSameToSame(int lid_recv, int dest_rank, Variable *var, bool TryRecvSameToSame(int lid_recv, int send_rank, Variable *var, MeshBlock *pmb, Mesh *pmesh) { MPI_Comm comm = pmesh->GetMPIComm(var->label()); - int tag = CreateAMRMPITag(lid_recv, 0, 0, 0); - + int tag = CreateAMRMPITag(lid_recv); int test; MPI_Status status; PARTHENON_MPI_CHECK(MPI_Iprobe(send_rank, tag, comm, &test, &status)); @@ -574,8 +560,7 @@ void Mesh::UpdateMeshBlockTree(int &nnew, int &ndel) { int lk = !mesh_size.symmetry(X3DIR); int lj = !mesh_size.symmetry(X2DIR); for (int n = 0; n < tnderef; n++) { - if ((lderef[n].lx1() & 1LL) == 0LL && (lderef[n].lx2() & 1LL) == 0LL && - (lderef[n].lx3() & 1LL) == 0LL) { + if (lderef[n].IsLowerLeftCornerOfParent()) { int r = n, rr = 0; // offsets from "bottom left" corner of parent block for (std::int64_t ok = 0; ok <= lk; ok++) { From 2b8ab9ba12076caa6040d77bb462d2a6ffd7db8f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 24 Apr 2025 14:45:24 -0600 Subject: [PATCH 09/13] fix bugs --- src/defs.hpp | 2 +- src/mesh/forest/logical_location.hpp | 12 ++++++------ src/mesh/mesh-amr_loadbalance.cpp | 2 ++ 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/defs.hpp b/src/defs.hpp index 0983438c294f..e6ade682a687 100644 --- a/src/defs.hpp +++ b/src/defs.hpp @@ -181,7 +181,7 @@ inline std::array GetOffsetsFromBoundaryFace(BoundaryFace face) { return {0, 0, -1}; case BoundaryFace::outer_x3: return {0, 0, 1}; - case default: + default: PARTHENON_FAIL("Asking for offsets for an invalid BoundaryFace."); } return {0, 0, 0}; diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 10a4586ac0d1..e1417d5d49ac 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -123,18 +123,18 @@ class LogicalLocation { // aggregate and POD type return {(lx1() & 1LL) == 1LL, (lx2() & 1LL) == 1LL, (lx3() & 1LL) == 1LL}; } - bool IsOnTreeBoundary(int ox1, int ox2, int ox3) { - const int nup = 1 << std::max(level(), 0) - 1; + 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 && loc.lx1() == 0) || (ox1 == 1 && loc.lx1() == nup); + (ox1 == 0) || ((ox1 == -1) && (lx1() == 0)) || ((ox1 == 1) && (lx1() == nup)); const bool bound2 = - (ox2 == 0) || (ox2 == -1 && loc.lx2() == 0) || (ox2 == 1 && loc.lx2() == nup); + (ox2 == 0) || ((ox2 == -1) && (lx2() == 0)) || ((ox2 == 1) && (lx2() == nup)); const bool bound3 = - (ox3 == 0) || (ox3 == -1 && loc.lx3() == 0) || (ox3 == 1 && loc.lx3() == nup); + (ox3 == 0) || ((ox3 == -1) && (lx3() == 0)) || ((ox3 == 1) && (lx3() == nup)); return bound1 && bound2 && bound3; } - bool IsOnTreeBoundary(BoundaryFace face) { + bool IsOnTreeBoundary(BoundaryFace face) const { const auto [ox1, ox2, ox3] = GetOffsetsFromBoundaryFace(face); return IsOnTreeBoundary(ox1, ox2, ox3); } diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index dddbfbf1b70b..9f8039e9d7bd 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -120,6 +120,7 @@ bool TryRecvCoarseToFine(int lid_recv, int send_rank, const LogicalLocation &fin IndexRange jb_int = cellbounds.GetBoundsJ(IndexDomain::interior, te); IndexRange kb_int = cellbounds.GetBoundsK(IndexDomain::interior, te); + const auto [ox1, ox2, ox3] = fine_loc.GetLocationInParent(); const int ks = (ox3 == 0) ? 0 : (kb_int.e - kb_int.s + 1) / 2; const int js = (ox2 == 0) ? 0 : (jb_int.e - jb_int.s + 1) / 2; const int is = (ox1 == 0) ? 0 : (ib_int.e - ib_int.s + 1) / 2; @@ -203,6 +204,7 @@ bool TryRecvFineToCoarse(int lid_recv, int send_rank, const LogicalLocation &fin // space if fine block is on the left side of a direction. I think this // should work fine even if the ownership model is changed elsewhere, since // the fine blocks should be consistent in their shared elements at this point + const auto [ox1, ox2, ox3] = fine_loc.GetLocationInParent(); if (ox3 == 0 && ndim > 2) kb.e -= TopologicalOffsetK(te); if (ox2 == 0 && ndim > 1) jb.e -= TopologicalOffsetJ(te); if (ox1 == 0) ib.e -= TopologicalOffsetI(te); From 48b65fc58e90a2b970f3f38bdafc84efca90b767 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 28 Apr 2025 15:33:49 -0400 Subject: [PATCH 10/13] make sure colorbar label generated --- .../python/packages/parthenon_tools/parthenon_tools/movie2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py b/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py index 83c4be2ab2c1..52abf6181b24 100644 --- a/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py +++ b/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py @@ -296,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) From 37efb1fa76c09928998347100054c999b87b125d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 30 Apr 2025 10:43:11 -0400 Subject: [PATCH 11/13] Update doc/sphinx/src/outputs.rst Co-authored-by: Philipp Grete --- doc/sphinx/src/outputs.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index ace250b81411..d142abe8c170 100644 --- a/doc/sphinx/src/outputs.rst +++ b/doc/sphinx/src/outputs.rst @@ -170,7 +170,7 @@ written. Changing output cadence from command line when restarting from file -------------------------------------------------------------------- -When restarting from a restart file, you can change any paremeter +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 From 2f11fcb6899ce6652c71d7097ad57f55a12fa4df Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 30 Apr 2025 10:43:41 -0400 Subject: [PATCH 12/13] Update doc/sphinx/src/outputs.rst Co-authored-by: Philipp Grete --- doc/sphinx/src/outputs.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index d142abe8c170..780ca70699b8 100644 --- a/doc/sphinx/src/outputs.rst +++ b/doc/sphinx/src/outputs.rst @@ -177,6 +177,8 @@ special care is required, as the code stores a ``next_time`` and a 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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ From 1aaa8b32faebb11f68452a4e5cdd0b6abf02b37d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 30 Apr 2025 11:12:45 -0400 Subject: [PATCH 13/13] add resport_asymmetry --- .../parthenon_tools/compute_asymmetry.py | 13 ++--- .../parthenon_tools/report_asymmetry.py | 54 +++++++++++++++++++ 2 files changed, 61 insertions(+), 6 deletions(-) create mode 100755 scripts/python/packages/parthenon_tools/parthenon_tools/report_asymmetry.py diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py b/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py index 7cafb439c920..a0de84a448e4 100644 --- a/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py +++ b/scripts/python/packages/parthenon_tools/parthenon_tools/compute_asymmetry.py @@ -36,18 +36,19 @@ def compute_asymmetry(f, varname): for b in range(ylocs.shape[0]): bb = matches[b] if np.any(ylocs[b] >= 0): - var_diff[b] = var[b] - np.flip(var[bb], axis=-2) - var_diff[bb] = var[bb] - np.flip(var[b], axis=-2) - - if np.any(np.abs(var_diff) > 1): - print("var_diff > 1e-8!", np.where(np.abs(var_diff) > 1)) + 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.", + 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") diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/report_asymmetry.py b/scripts/python/packages/parthenon_tools/parthenon_tools/report_asymmetry.py new file mode 100755 index 000000000000..614a203702dc --- /dev/null +++ b/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!")