diff --git a/CHANGELOG.md b/CHANGELOG.md index 7dc28b51eb00..c2a19d8b7b2a 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 diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index 6162e49242d3..780ca70699b8 100644 --- a/doc/sphinx/src/outputs.rst +++ b/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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 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..a0de84a448e4 --- /dev/null +++ b/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 diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py b/scripts/python/packages/parthenon_tools/parthenon_tools/movie2d.py index 02509746604c..52abf6181b24 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. # @@ -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) 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!") diff --git a/src/defs.hpp b/src/defs.hpp index 552e9cb32087..e6ade682a687 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}; + 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..e1417d5d49ac 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 { + 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 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 00d768720988..9f8039e9d7bd 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -1,13 +1,13 @@ //======================================================================================== // 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 // 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 @@ -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 @@ -124,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; @@ -154,12 +151,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 +167,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 @@ -216,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); @@ -252,7 +241,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 +275,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)); @@ -502,8 +490,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 +504,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 +530,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 +546,32 @@ 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) { + if (lderef[n].IsLowerLeftCornerOfParent()) { 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 +586,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 +598,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; } //----------------------------------------------------------------------------------------