From 1dafe264a0f43c4f5171963482417b7bb5a0ddb4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 20 Mar 2024 15:37:17 -0600 Subject: [PATCH 01/39] move some stuff around --- src/CMakeLists.txt | 4 ++-- src/bvals/neighbor_block.cpp | 2 +- src/defs.hpp | 2 +- src/mesh/forest/cell_center_offsets.hpp | 2 +- src/mesh/forest/forest.cpp | 2 +- src/mesh/forest/forest.hpp | 2 +- src/mesh/{ => forest}/logical_location.cpp | 27 +--------------------- src/mesh/{ => forest}/logical_location.hpp | 10 ++++---- src/mesh/forest/relative_orientation.cpp | 2 +- src/mesh/forest/relative_orientation.hpp | 2 +- src/mesh/forest/tree.cpp | 2 +- src/mesh/forest/tree.hpp | 2 +- 12 files changed, 16 insertions(+), 43 deletions(-) rename src/mesh/{ => forest}/logical_location.cpp (94%) rename src/mesh/{ => forest}/logical_location.hpp (97%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0a028442e651..55a26c3ac925 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -165,8 +165,8 @@ add_library(parthenon mesh/forest/relative_orientation.cpp mesh/forest/tree.hpp mesh/forest/tree.cpp - mesh/logical_location.cpp - mesh/logical_location.hpp + mesh/forest/logical_location.cpp + mesh/forest/logical_location.hpp mesh/mesh_refinement.cpp mesh/mesh_refinement.hpp mesh/mesh-gmg.cpp diff --git a/src/bvals/neighbor_block.cpp b/src/bvals/neighbor_block.cpp index 22b6b6dce80d..2d0c0b13a597 100644 --- a/src/bvals/neighbor_block.cpp +++ b/src/bvals/neighbor_block.cpp @@ -30,7 +30,7 @@ #include #include "globals.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "mesh/mesh.hpp" #include "utils/buffer_utils.hpp" #include "utils/error_checking.hpp" diff --git a/src/defs.hpp b/src/defs.hpp index 1403f985f3a3..e2ea6db20d6b 100644 --- a/src/defs.hpp +++ b/src/defs.hpp @@ -25,7 +25,7 @@ #include "basic_types.hpp" #include "config.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "parthenon_arrays.hpp" namespace parthenon { diff --git a/src/mesh/forest/cell_center_offsets.hpp b/src/mesh/forest/cell_center_offsets.hpp index 98e18d8f2d31..a1237a514201 100644 --- a/src/mesh/forest/cell_center_offsets.hpp +++ b/src/mesh/forest/cell_center_offsets.hpp @@ -25,7 +25,7 @@ #include "basic_types.hpp" #include "defs.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index 9ee4280a4205..60cc149af22b 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -29,7 +29,7 @@ #include "mesh/forest/forest.hpp" #include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index 4b9dd443e71d..a31f17dcd288 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -26,7 +26,7 @@ #include "basic_types.hpp" #include "defs.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/logical_location.cpp b/src/mesh/forest/logical_location.cpp similarity index 94% rename from src/mesh/logical_location.cpp rename to src/mesh/forest/logical_location.cpp index 672a4845a549..71559a069dbd 100644 --- a/src/mesh/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -30,7 +30,7 @@ #include #include -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/error_checking.hpp" #include "utils/morton_number.hpp" @@ -53,31 +53,6 @@ bool LogicalLocation::Contains(const LogicalLocation &containee) const { return (shifted_lx1 == lx1()) && (shifted_lx2 == lx2()) && (shifted_lx3 == lx3()); } -std::array LogicalLocation::GetOffset(const LogicalLocation &neighbor, - const RootGridInfo &rg_info) const { - std::array offset; - const int level_diff_1 = std::max(neighbor.level() - level(), 0); - const int level_diff_2 = std::max(level() - neighbor.level(), 0); - const int n_per_root_block = 1 << (std::min(level(), neighbor.level()) - rg_info.level); - const int root_block_per_n = - 1 << std::max(rg_info.level - std::min(level(), neighbor.level()), 0); - for (int i = 0; i < 3; ++i) { - offset[i] = (neighbor.l(i) >> level_diff_1) - (l(i) >> level_diff_2); - if (rg_info.periodic[i]) { - int n_cells_level = std::max(n_per_root_block * rg_info.n[i], 1); - if (root_block_per_n > 1) - n_cells_level = - rg_info.n[i] / root_block_per_n + (rg_info.n[i] % root_block_per_n != 0); - if (std::abs(offset[i]) > (n_cells_level / 2)) { - offset[i] %= n_cells_level; - offset[i] += offset[i] > 0 ? -n_cells_level : n_cells_level; - } - } - } - - return offset; -} - std::array LogicalLocation::GetSameLevelOffsetsForest(const LogicalLocation &neighbor) const { std::array offsets; diff --git a/src/mesh/logical_location.hpp b/src/mesh/forest/logical_location.hpp similarity index 97% rename from src/mesh/logical_location.hpp rename to src/mesh/forest/logical_location.hpp index 45cae77056b7..f2f169c2457b 100644 --- a/src/mesh/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -49,6 +49,8 @@ struct std::hash { namespace parthenon { +// TODO(LFR): This can go away once MG is fixed for forests, and probably any routine that +// depends on it. struct RootGridInfo { int level; std::array n; @@ -127,8 +129,6 @@ class LogicalLocation { // aggregate and POD type bool Contains(const LogicalLocation &containee) const; - std::array GetOffset(const LogicalLocation &neighbor, - const RootGridInfo &rg_info = RootGridInfo()) const; // TODO(LFR): Remove the corresponding non-forest routine once GMG is working std::array GetSameLevelOffsetsForest(const LogicalLocation &neighbor) const; std::array, 3> GetSameLevelOffsets(const LogicalLocation &neighbor, @@ -152,8 +152,7 @@ class LogicalLocation { // aggregate and POD type } LogicalLocation - GetSameLevelNeighbor(int ox1, int ox2, int ox3, - const RootGridInfo &rg_info = RootGridInfo()) const { + GetSameLevelNeighbor(int ox1, int ox2, int ox3) const { return LogicalLocation(tree(), level(), lx1() + ox1, lx2() + ox2, lx3() + ox3); } @@ -175,8 +174,7 @@ class LogicalLocation { // aggregate and POD type // 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 // correction is combined with boundary communication. - auto GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, int ox2, int ox3, - const RootGridInfo &rg_info = RootGridInfo()) const { + auto GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, int ox2, int ox3) const { // The neighbor block struct should only use the first two, but we have three to allow // for this being a parent of neighbor, this should be checked for elsewhere std::array f{0, 0, 0}; diff --git a/src/mesh/forest/relative_orientation.cpp b/src/mesh/forest/relative_orientation.cpp index e6b2b531dec2..6250ef5cb3f9 100644 --- a/src/mesh/forest/relative_orientation.cpp +++ b/src/mesh/forest/relative_orientation.cpp @@ -28,7 +28,7 @@ #include "defs.hpp" #include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/relative_orientation.hpp b/src/mesh/forest/relative_orientation.hpp index 7bca5006b666..7bb711da47bc 100644 --- a/src/mesh/forest/relative_orientation.hpp +++ b/src/mesh/forest/relative_orientation.hpp @@ -26,7 +26,7 @@ #include "basic_types.hpp" #include "defs.hpp" #include "mesh/forest/cell_center_offsets.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index 06c59cd9346b..fc1d6ea1013c 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -28,7 +28,7 @@ #include "defs.hpp" #include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 4509e2848dd2..46192e25e09d 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -27,7 +27,7 @@ #include "defs.hpp" #include "mesh/forest/cell_center_offsets.hpp" #include "mesh/forest/relative_orientation.hpp" -#include "mesh/logical_location.hpp" +#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" From 022901200078c143e9a1ba74967e3258b416055e Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 20 Mar 2024 16:00:09 -0600 Subject: [PATCH 02/39] split into more, smaller files --- src/CMakeLists.txt | 2 + src/bvals/neighbor_block.hpp | 2 + src/mesh/forest/block_ownership.cpp | 183 +++++++++++++++++++++++++++ src/mesh/forest/block_ownership.hpp | 108 ++++++++++++++++ src/mesh/forest/forest.hpp | 8 ++ src/mesh/forest/logical_location.cpp | 143 --------------------- src/mesh/forest/logical_location.hpp | 64 ---------- src/mesh/mesh-gmg.cpp | 6 +- src/utils/indexer.hpp | 1 + 9 files changed, 306 insertions(+), 211 deletions(-) create mode 100644 src/mesh/forest/block_ownership.cpp create mode 100644 src/mesh/forest/block_ownership.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 55a26c3ac925..8f8aeb47479d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -158,6 +158,8 @@ add_library(parthenon mesh/amr_loadbalance.cpp mesh/domain.hpp + mesh/forest/block_ownership.cpp + mesh/forest/block_ownership.hpp mesh/forest/cell_center_offsets.hpp mesh/forest/forest.cpp mesh/forest/forest.hpp diff --git a/src/bvals/neighbor_block.hpp b/src/bvals/neighbor_block.hpp index b6ee753a6758..99b849985142 100644 --- a/src/bvals/neighbor_block.hpp +++ b/src/bvals/neighbor_block.hpp @@ -29,6 +29,8 @@ #include "parthenon_mpi.hpp" #include "defs.hpp" +#include "mesh/forest/block_ownership.hpp" +#include "mesh/forest/logical_location.hpp" #include "parthenon_arrays.hpp" #include "utils/error_checking.hpp" diff --git a/src/mesh/forest/block_ownership.cpp b/src/mesh/forest/block_ownership.cpp new file mode 100644 index 000000000000..967ef5829659 --- /dev/null +++ b/src/mesh/forest/block_ownership.cpp @@ -0,0 +1,183 @@ +//======================================================================================== +// 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 +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2020-2023 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 2020-2023. 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. +//======================================================================================== + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "mesh/forest/block_ownership.hpp" +#include "mesh/forest/logical_location.hpp" +#include "utils/error_checking.hpp" +#include "utils/morton_number.hpp" + +namespace parthenon { + +// TODO(LFR): Remove this +block_ownership_t +DetermineOwnership(const LogicalLocation &main_block, + const std::unordered_set &allowed_neighbors, + const RootGridInfo &rg_info, + const std::unordered_set &newly_refined) { + block_ownership_t main_owns; + + auto ownership_level = [&](const LogicalLocation &a) { + // Newly-refined blocks are treated as higher-level than blocks at their + // parent level, but lower-level than previously-refined blocks at their + // current level. + if (newly_refined.count(a)) return 2 * a.level() - 1; + return 2 * a.level(); + }; + + auto ownership_less_than = [ownership_level](const LogicalLocation &a, + const LogicalLocation &b) { + // Ownership is first determined by block with the highest level, then by maximum + // Morton number this is reversed in precedence from the normal comparators where + // Morton number takes precedence + if (ownership_level(a) == ownership_level(b)) return a.morton() < b.morton(); + return ownership_level(a) < ownership_level(b); + }; + + for (int ox1 : {-1, 0, 1}) { + for (int ox2 : {-1, 0, 1}) { + for (int ox3 : {-1, 0, 1}) { + main_owns(ox1, ox2, ox3) = true; + for (auto &n : allowed_neighbors) { + if (ownership_less_than(main_block, n) && + main_block.IsNeighborOfTE(n, ox1, ox2, ox3, rg_info)) { + main_owns(ox1, ox2, ox3) = false; + break; + } + } + } + } + } + return main_owns; +} + +block_ownership_t +DetermineOwnershipForest(const LogicalLocation &main_block, + const std::vector &allowed_neighbors, + const std::unordered_set &newly_refined) { + block_ownership_t main_owns; + + auto ownership_level = [&](const LogicalLocation &a) { + // Newly-refined blocks are treated as higher-level than blocks at their + // parent level, but lower-level than previously-refined blocks at their + // current level. + if (newly_refined.count(a)) return 2 * a.level() - 1; + return 2 * a.level(); + }; + + auto ownership_less_than = [ownership_level](const LogicalLocation &a, + const LogicalLocation &b) { + // Ownership is first determined by block with the highest level, then by maximum + // (tree, Morton) number this is reversed in precedence from the normal comparators + // where (tree, Morton) number takes precedence + if (ownership_level(a) != ownership_level(b)) + return ownership_level(a) < ownership_level(b); + if (a.tree() != b.tree()) return a.tree() < b.tree(); + return a.morton() < b.morton(); + }; + + for (int ox1 : {-1, 0, 1}) { + for (int ox2 : {-1, 0, 1}) { + for (int ox3 : {-1, 0, 1}) { + main_owns(ox1, ox2, ox3) = true; + for (const auto &n : allowed_neighbors) { + if (ownership_less_than(main_block, n.global_loc) && + main_block.IsNeighborOfTEForest(n.origin_loc, {ox1, ox2, ox3})) { + main_owns(ox1, ox2, ox3) = false; + break; + } + } + } + } + } + return main_owns; +} + +// Given a topological element, ownership array of the sending block, and offset indices +// defining the location of an index region within the block (i.e. the ghost zones passed +// across the x-face or the ghost zones passed across the z-edge), return the index range +// masking array required for masking out unowned regions of the index space. ox? defines +// buffer location on the owner block +block_ownership_t +GetIndexRangeMaskFromOwnership(TopologicalElement el, + const block_ownership_t &sender_ownership, int ox1, + int ox2, int ox3) { + using vp_t = std::vector>; + + // Transform general block ownership to element ownership over entire block. For + // instance, x-faces only care about block ownership in the x-direction First index of + // the pair is the element index and the second index is the block index that is copied + // to that element index + block_ownership_t element_ownership = sender_ownership; + auto x1_idxs = TopologicalOffsetI(el) ? vp_t{{-1, -1}, {0, 0}, {1, 1}} + : vp_t{{-1, 0}, {0, 0}, {1, 0}}; + auto x2_idxs = TopologicalOffsetJ(el) ? vp_t{{-1, -1}, {0, 0}, {1, 1}} + : vp_t{{-1, 0}, {0, 0}, {1, 0}}; + auto x3_idxs = TopologicalOffsetK(el) ? vp_t{{-1, -1}, {0, 0}, {1, 1}} + : vp_t{{-1, 0}, {0, 0}, {1, 0}}; + for (auto [iel, ibl] : x1_idxs) { + for (auto [jel, jbl] : x2_idxs) { + for (auto [kel, kbl] : x3_idxs) { + element_ownership(iel, jel, kel) = sender_ownership(ibl, jbl, kbl); + } + } + } + + // Now, the ownership status is correct for the entire interior index range of the + // block, but the offsets ox? define a subset of these indices (e.g. one edge of the + // interior). Therefore, we need to set the index ownership to true for edges of the + // index range that are contained in the interior of the sending block + if (ox1 != 0) { + for (auto j : {-1, 0, 1}) { + for (auto k : {-1, 0, 1}) { + element_ownership(-ox1, j, k) = element_ownership(0, j, k); + } + } + } + if (ox2 != 0) { + for (auto i : {-1, 0, 1}) { + for (auto k : {-1, 0, 1}) { + element_ownership(i, -ox2, k) = element_ownership(i, 0, k); + } + } + } + if (ox3 != 0) { + for (auto i : {-1, 0, 1}) { + for (auto j : {-1, 0, 1}) { + element_ownership(i, j, -ox3) = element_ownership(i, j, 0); + } + } + } + + return element_ownership; +} + +} // namespace parthenon diff --git a/src/mesh/forest/block_ownership.hpp b/src/mesh/forest/block_ownership.hpp new file mode 100644 index 000000000000..90c20ff1d057 --- /dev/null +++ b/src/mesh/forest/block_ownership.hpp @@ -0,0 +1,108 @@ +//======================================================================================== +// 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 +//======================================================================================== +// Parthenon performance portable AMR framework +// Copyright(C) 2020-2024 The Parthenon collaboration +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +// (C) (or copyright) 2020-2024. 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. +//======================================================================================== +#ifndef MESH_FOREST_BLOCK_OWNERSHIP_HPP_ +#define MESH_FOREST_BLOCK_OWNERSHIP_HPP_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "basic_types.hpp" +#include "mesh/forest/logical_location.hpp" +#include "utils/error_checking.hpp" +#include "utils/morton_number.hpp" + +namespace parthenon { + +struct block_ownership_t { + public: + KOKKOS_FORCEINLINE_FUNCTION + const bool &operator()(int ox1, int ox2, int ox3) const { + return ownership[ox1 + 1][ox2 + 1][ox3 + 1]; + } + KOKKOS_FORCEINLINE_FUNCTION + bool &operator()(int ox1, int ox2, int ox3) { + return ownership[ox1 + 1][ox2 + 1][ox3 + 1]; + } + + KOKKOS_FORCEINLINE_FUNCTION + block_ownership_t() : block_ownership_t(false) {} + + KOKKOS_FORCEINLINE_FUNCTION + explicit block_ownership_t(bool value) : initialized(false) { + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + ownership[i][j][k] = value; + } + } + } + } + + bool initialized; + + bool operator==(const block_ownership_t &rhs) const { + bool same = initialized == rhs.initialized; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + same = same && (ownership[i][j][k] == rhs.ownership[i][j][k]); + } + } + } + return same; + } + + private: + bool ownership[3][3][3]; +}; + +block_ownership_t +DetermineOwnership(const LogicalLocation &main_block, + const std::unordered_set &allowed_neighbors, + const RootGridInfo &rg_info = RootGridInfo(), + const std::unordered_set &newly_refined = {}); + +block_ownership_t +DetermineOwnershipForest(const LogicalLocation &main_block, + const std::vector &allowed_neighbors, + const std::unordered_set &newly_refined = {}); + +// Given a topological element, ownership array of the sending block, and offset indices +// defining the location of an index region within the block (i.e. the ghost zones passed +// across the x-face or the ghost zones passed across the z-edge), return the index range +// masking array required for masking out unowned regions of the index space. ox? defines +// buffer location on the owner block +block_ownership_t +GetIndexRangeMaskFromOwnership(TopologicalElement el, + const block_ownership_t &sender_ownership, int ox1, + int ox2, int ox3); + +} // namespace parthenon + +#endif // MESH_FOREST_BLOCK_OWNERSHIP_HPP_ diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index a31f17dcd288..1c6d29534c54 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -145,6 +145,14 @@ class Forest { std::array mesh_bcs); }; +struct NeighborLocation { + NeighborLocation(const LogicalLocation &g, const LogicalLocation &o) + : global_loc(g), origin_loc(o) {} + LogicalLocation global_loc; // Global location of neighboring block + LogicalLocation + origin_loc; // Logical location of neighboring block in index space of origin block +}; + } // namespace forest } // namespace parthenon diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index 71559a069dbd..7c85d801a492 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -333,147 +333,4 @@ std::unordered_set LogicalLocation::GetPossibleNeighborsImpl( return unique_locs; } -// TODO(LFR): Remove this -block_ownership_t -DetermineOwnership(const LogicalLocation &main_block, - const std::unordered_set &allowed_neighbors, - const RootGridInfo &rg_info, - const std::unordered_set &newly_refined) { - block_ownership_t main_owns; - - auto ownership_level = [&](const LogicalLocation &a) { - // Newly-refined blocks are treated as higher-level than blocks at their - // parent level, but lower-level than previously-refined blocks at their - // current level. - if (newly_refined.count(a)) return 2 * a.level() - 1; - return 2 * a.level(); - }; - - auto ownership_less_than = [ownership_level](const LogicalLocation &a, - const LogicalLocation &b) { - // Ownership is first determined by block with the highest level, then by maximum - // Morton number this is reversed in precedence from the normal comparators where - // Morton number takes precedence - if (ownership_level(a) == ownership_level(b)) return a.morton() < b.morton(); - return ownership_level(a) < ownership_level(b); - }; - - for (int ox1 : {-1, 0, 1}) { - for (int ox2 : {-1, 0, 1}) { - for (int ox3 : {-1, 0, 1}) { - main_owns(ox1, ox2, ox3) = true; - for (auto &n : allowed_neighbors) { - if (ownership_less_than(main_block, n) && - main_block.IsNeighborOfTE(n, ox1, ox2, ox3, rg_info)) { - main_owns(ox1, ox2, ox3) = false; - break; - } - } - } - } - } - return main_owns; -} - -block_ownership_t -DetermineOwnershipForest(const LogicalLocation &main_block, - const std::vector &allowed_neighbors, - const std::unordered_set &newly_refined) { - block_ownership_t main_owns; - - auto ownership_level = [&](const LogicalLocation &a) { - // Newly-refined blocks are treated as higher-level than blocks at their - // parent level, but lower-level than previously-refined blocks at their - // current level. - if (newly_refined.count(a)) return 2 * a.level() - 1; - return 2 * a.level(); - }; - - auto ownership_less_than = [ownership_level](const LogicalLocation &a, - const LogicalLocation &b) { - // Ownership is first determined by block with the highest level, then by maximum - // (tree, Morton) number this is reversed in precedence from the normal comparators - // where (tree, Morton) number takes precedence - if (ownership_level(a) != ownership_level(b)) - return ownership_level(a) < ownership_level(b); - if (a.tree() != b.tree()) return a.tree() < b.tree(); - return a.morton() < b.morton(); - }; - - for (int ox1 : {-1, 0, 1}) { - for (int ox2 : {-1, 0, 1}) { - for (int ox3 : {-1, 0, 1}) { - main_owns(ox1, ox2, ox3) = true; - for (const auto &n : allowed_neighbors) { - if (ownership_less_than(main_block, n.global_loc) && - main_block.IsNeighborOfTEForest(n.origin_loc, {ox1, ox2, ox3})) { - main_owns(ox1, ox2, ox3) = false; - break; - } - } - } - } - } - return main_owns; -} - -// Given a topological element, ownership array of the sending block, and offset indices -// defining the location of an index region within the block (i.e. the ghost zones passed -// across the x-face or the ghost zones passed across the z-edge), return the index range -// masking array required for masking out unowned regions of the index space. ox? defines -// buffer location on the owner block -block_ownership_t -GetIndexRangeMaskFromOwnership(TopologicalElement el, - const block_ownership_t &sender_ownership, int ox1, - int ox2, int ox3) { - using vp_t = std::vector>; - - // Transform general block ownership to element ownership over entire block. For - // instance, x-faces only care about block ownership in the x-direction First index of - // the pair is the element index and the second index is the block index that is copied - // to that element index - block_ownership_t element_ownership = sender_ownership; - auto x1_idxs = TopologicalOffsetI(el) ? vp_t{{-1, -1}, {0, 0}, {1, 1}} - : vp_t{{-1, 0}, {0, 0}, {1, 0}}; - auto x2_idxs = TopologicalOffsetJ(el) ? vp_t{{-1, -1}, {0, 0}, {1, 1}} - : vp_t{{-1, 0}, {0, 0}, {1, 0}}; - auto x3_idxs = TopologicalOffsetK(el) ? vp_t{{-1, -1}, {0, 0}, {1, 1}} - : vp_t{{-1, 0}, {0, 0}, {1, 0}}; - for (auto [iel, ibl] : x1_idxs) { - for (auto [jel, jbl] : x2_idxs) { - for (auto [kel, kbl] : x3_idxs) { - element_ownership(iel, jel, kel) = sender_ownership(ibl, jbl, kbl); - } - } - } - - // Now, the ownership status is correct for the entire interior index range of the - // block, but the offsets ox? define a subset of these indices (e.g. one edge of the - // interior). Therefore, we need to set the index ownership to true for edges of the - // index range that are contained in the interior of the sending block - if (ox1 != 0) { - for (auto j : {-1, 0, 1}) { - for (auto k : {-1, 0, 1}) { - element_ownership(-ox1, j, k) = element_ownership(0, j, k); - } - } - } - if (ox2 != 0) { - for (auto i : {-1, 0, 1}) { - for (auto k : {-1, 0, 1}) { - element_ownership(i, -ox2, k) = element_ownership(i, 0, k); - } - } - } - if (ox3 != 0) { - for (auto i : {-1, 0, 1}) { - for (auto j : {-1, 0, 1}) { - element_ownership(i, j, -ox3) = element_ownership(i, j, 0); - } - } - } - - return element_ownership; -} - } // namespace parthenon diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index f2f169c2457b..d0c9894284cb 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -227,49 +227,6 @@ inline bool operator!=(const LogicalLocation &lhs, const LogicalLocation &rhs) { return !(lhs == rhs); } -struct block_ownership_t { - public: - KOKKOS_FORCEINLINE_FUNCTION - const bool &operator()(int ox1, int ox2, int ox3) const { - return ownership[ox1 + 1][ox2 + 1][ox3 + 1]; - } - KOKKOS_FORCEINLINE_FUNCTION - bool &operator()(int ox1, int ox2, int ox3) { - return ownership[ox1 + 1][ox2 + 1][ox3 + 1]; - } - - KOKKOS_FORCEINLINE_FUNCTION - block_ownership_t() : block_ownership_t(false) {} - - KOKKOS_FORCEINLINE_FUNCTION - explicit block_ownership_t(bool value) : initialized(false) { - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - for (int k = 0; k < 3; ++k) { - ownership[i][j][k] = value; - } - } - } - } - - bool initialized; - - bool operator==(const block_ownership_t &rhs) const { - bool same = initialized == rhs.initialized; - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - for (int k = 0; k < 3; ++k) { - same = same && (ownership[i][j][k] == rhs.ownership[i][j][k]); - } - } - } - return same; - } - - private: - bool ownership[3][3][3]; -}; - struct NeighborLocation { NeighborLocation(const LogicalLocation &g, const LogicalLocation &o) : global_loc(g), origin_loc(o) {} @@ -278,27 +235,6 @@ struct NeighborLocation { origin_loc; // Logical location of neighboring block in index space of origin block }; -block_ownership_t -DetermineOwnership(const LogicalLocation &main_block, - const std::unordered_set &allowed_neighbors, - const RootGridInfo &rg_info = RootGridInfo(), - const std::unordered_set &newly_refined = {}); - -block_ownership_t -DetermineOwnershipForest(const LogicalLocation &main_block, - const std::vector &allowed_neighbors, - const std::unordered_set &newly_refined = {}); - -// Given a topological element, ownership array of the sending block, and offset indices -// defining the location of an index region within the block (i.e. the ghost zones passed -// across the x-face or the ghost zones passed across the z-edge), return the index range -// masking array required for masking out unowned regions of the index space. ox? defines -// buffer location on the owner block -block_ownership_t -GetIndexRangeMaskFromOwnership(TopologicalElement el, - const block_ownership_t &sender_ownership, int ox1, - int ox2, int ox3); - } // namespace parthenon inline std::size_t std::hash::operator()( diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 84398a90237f..ada986fcdc25 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -164,10 +164,8 @@ void Mesh::SetSameLevelNeighbors( } else if (connect_indicator == 3) { nc = NeighborConnect::corner; } - auto f = loc.GetAthenaXXFaceOffsets(pos_neighbor_location, ox1, ox2, ox3, - root_grid); - auto fn = pos_neighbor_location.GetAthenaXXFaceOffsets(loc, -ox1, -ox2, - -ox3, root_grid); + auto f = loc.GetAthenaXXFaceOffsets(pos_neighbor_location, ox1, ox2, ox3); + auto fn = pos_neighbor_location.GetAthenaXXFaceOffsets(loc, -ox1, -ox2, -ox3); int bid = buffer_id.GetID(ox1, ox2, ox3, f[0], f[1]); int tid = buffer_id.GetID(-ox1, -ox2, -ox3, fn[0], fn[1]); neighbor_list->emplace_back( diff --git a/src/utils/indexer.hpp b/src/utils/indexer.hpp index 0902019ff4c2..21951fc1ceab 100644 --- a/src/utils/indexer.hpp +++ b/src/utils/indexer.hpp @@ -19,6 +19,7 @@ #include #include +#include "mesh/forest/block_ownership.hpp" #include "utils/concepts_lite.hpp" #include "utils/utils.hpp" From 7a43c76c44e151723effe447ba88d55bf2698e7f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 20 Mar 2024 16:01:10 -0600 Subject: [PATCH 03/39] format and lint --- src/mesh/forest/forest.cpp | 2 +- src/mesh/forest/forest.hpp | 2 +- src/mesh/forest/logical_location.hpp | 14 +++++++------- src/mesh/forest/relative_orientation.cpp | 2 +- src/mesh/forest/tree.cpp | 2 +- src/mesh/forest/tree.hpp | 2 +- src/mesh/mesh-gmg.cpp | 3 ++- 7 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index 60cc149af22b..0b4f191450da 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -27,9 +27,9 @@ #include "basic_types.hpp" #include "defs.hpp" #include "mesh/forest/forest.hpp" +#include "mesh/forest/logical_location.hpp" #include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index 1c6d29534c54..2a24e74c439b 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -25,8 +25,8 @@ #include "basic_types.hpp" #include "defs.hpp" -#include "mesh/forest/tree.hpp" #include "mesh/forest/logical_location.hpp" +#include "mesh/forest/tree.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index d0c9894284cb..a81a07a0d952 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -18,8 +18,8 @@ // 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. //======================================================================================== -#ifndef MESH_LOGICAL_LOCATION_HPP_ -#define MESH_LOGICAL_LOCATION_HPP_ +#ifndef MESH_FOREST_LOGICAL_LOCATION_HPP_ +#define MESH_FOREST_LOGICAL_LOCATION_HPP_ #include #include @@ -49,7 +49,7 @@ struct std::hash { namespace parthenon { -// TODO(LFR): This can go away once MG is fixed for forests, and probably any routine that +// TODO(LFR): This can go away once MG is fixed for forests, and probably any routine that // depends on it. struct RootGridInfo { int level; @@ -151,8 +151,7 @@ class LogicalLocation { // aggregate and POD type return NeighborFindingImpl(in, std::array{ox1, ox2, ox3}, rg_info); } - LogicalLocation - GetSameLevelNeighbor(int ox1, int ox2, int ox3) const { + LogicalLocation GetSameLevelNeighbor(int ox1, int ox2, int ox3) const { return LogicalLocation(tree(), level(), lx1() + ox1, lx2() + ox2, lx3() + ox3); } @@ -174,7 +173,8 @@ class LogicalLocation { // aggregate and POD type // 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 // correction is combined with boundary communication. - auto GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, int ox2, int ox3) const { + auto GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, int ox2, + int ox3) const { // The neighbor block struct should only use the first two, but we have three to allow // for this being a parent of neighbor, this should be checked for elsewhere std::array f{0, 0, 0}; @@ -245,4 +245,4 @@ inline std::size_t std::hash::operator()( return key.morton().bits[0]; } -#endif // MESH_LOGICAL_LOCATION_HPP_ +#endif // MESH_FOREST_LOGICAL_LOCATION_HPP_ diff --git a/src/mesh/forest/relative_orientation.cpp b/src/mesh/forest/relative_orientation.cpp index 6250ef5cb3f9..38492705861b 100644 --- a/src/mesh/forest/relative_orientation.cpp +++ b/src/mesh/forest/relative_orientation.cpp @@ -26,9 +26,9 @@ #include "basic_types.hpp" #include "defs.hpp" +#include "mesh/forest/logical_location.hpp" #include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index fc1d6ea1013c..06a935ccb178 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -26,9 +26,9 @@ #include "basic_types.hpp" #include "defs.hpp" +#include "mesh/forest/logical_location.hpp" #include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/tree.hpp" -#include "mesh/forest/logical_location.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 46192e25e09d..2b46493d6798 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -26,8 +26,8 @@ #include "basic_types.hpp" #include "defs.hpp" #include "mesh/forest/cell_center_offsets.hpp" -#include "mesh/forest/relative_orientation.hpp" #include "mesh/forest/logical_location.hpp" +#include "mesh/forest/relative_orientation.hpp" #include "utils/bit_hacks.hpp" #include "utils/indexer.hpp" diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index ada986fcdc25..05446c15b023 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -165,7 +165,8 @@ void Mesh::SetSameLevelNeighbors( nc = NeighborConnect::corner; } auto f = loc.GetAthenaXXFaceOffsets(pos_neighbor_location, ox1, ox2, ox3); - auto fn = pos_neighbor_location.GetAthenaXXFaceOffsets(loc, -ox1, -ox2, -ox3); + auto fn = + pos_neighbor_location.GetAthenaXXFaceOffsets(loc, -ox1, -ox2, -ox3); int bid = buffer_id.GetID(ox1, ox2, ox3, f[0], f[1]); int tid = buffer_id.GetID(-ox1, -ox2, -ox3, fn[0], fn[1]); neighbor_list->emplace_back( From a9810f1e09aec595909d397ab8c67861e8117f63 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 20 Mar 2024 16:46:09 -0600 Subject: [PATCH 04/39] clarify that amr_loadbalance implements mesh members --- src/CMakeLists.txt | 2 +- src/mesh/{amr_loadbalance.cpp => mesh-amr_loadbalance.cpp} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/mesh/{amr_loadbalance.cpp => mesh-amr_loadbalance.cpp} (100%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f8aeb47479d..9c332bf2713c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -156,7 +156,6 @@ add_library(parthenon interface/variable.cpp interface/variable.hpp - mesh/amr_loadbalance.cpp mesh/domain.hpp mesh/forest/block_ownership.cpp mesh/forest/block_ownership.hpp @@ -171,6 +170,7 @@ add_library(parthenon mesh/forest/logical_location.hpp mesh/mesh_refinement.cpp mesh/mesh_refinement.hpp + mesh/mesh-amr_loadbalance.cpp mesh/mesh-gmg.cpp mesh/mesh.cpp mesh/mesh.hpp diff --git a/src/mesh/amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp similarity index 100% rename from src/mesh/amr_loadbalance.cpp rename to src/mesh/mesh-amr_loadbalance.cpp From 21372e56d8035618ecdd256f21f2b0388ac7bc68 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 10:18:52 -0600 Subject: [PATCH 05/39] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fa0687d2a16f..e6f590fe6f70 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ ### Infrastructure (changes irrelevant to downstream codes) +- [[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 From dd6fbdca3f895cbf19feb0922a0b3bb1452f8c0f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 14:04:42 -0600 Subject: [PATCH 06/39] start on adding two-level composite grids --- src/mesh/forest/tree.cpp | 20 +++++++++++++++++++- src/mesh/forest/tree.hpp | 8 ++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index 06a935ccb178..a3a376ce80fb 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -37,7 +37,7 @@ namespace forest { Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSize domain, std::array bcs) - : my_id(id), ndim(ndim), domain(domain), boundary_conditions(bcs) { + : my_id(id), ndim(ndim), domain(domain), boundary_conditions(bcs), gmg_tlc_grids(50) { // Add internal and leaf nodes of the initial tree for (int l = 0; l <= root_level; ++l) { for (int k = 0; k < (ndim > 2 ? (1LL << l) : 1); ++k) { @@ -49,10 +49,14 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz } else { internal_nodes.emplace(my_id, l, i, j, k); } + gmg_tlc_grids[l].emplace(std::make_pair(LogicalLocation(my_id, l, i, j, k), + std::make_pair(-1, -1))) } } } } + // Pre-populate the next finest two-level composite grid with coarser blocks + gmg_tlc_grid[l + 1] = gmg_tlc_grid[l]; } int Tree::AddMeshBlock(const LogicalLocation &loc, bool enforce_proper_nesting) { @@ -95,6 +99,17 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { leaves.insert(std::make_pair(d, std::make_pair(gid_parent, -1))); } int nadded = daughters.size() - 1; + + // Update multigrid levels + auto& gmg_grid_1 = gmg_tlc_grids[ref_loc.level() + 1]; + auto& gmg_grid_2 = gmg_tlc_grids[ref_loc.level() + 2]; + if (gmg_grid_1.count(ref_loc)) gmg_grid_1.erase(ref_loc); + for (auto &d : daughters) { + // Insert as fine leaf blocks on coarser two-level grid + gmg_grid_1.insert(std::make_pair(d, std::make_pair(gid_parent, -1))); + // Insert as coarse leaf blocks on finer two-level grid + gmg_grid_2.insert(std::make_pair(d, std::make_pair(gid_parent, -1))); + } if (enforce_proper_nesting) { LogicalLocation parent = ref_loc.GetParent(); @@ -218,13 +233,16 @@ int Tree::Derefine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) } // Derefinement is ok + auto& gmg_grid = gmg_tlc_grids[ref_loc.level() + 1]; std::int64_t dgid = std::numeric_limits::max(); for (auto &d : daughters) { + gmg_grid.erase(d); auto node = leaves.extract(d); dgid = std::min(dgid, node.mapped().first); } internal_nodes.erase(ref_loc); leaves.insert(std::make_pair(ref_loc, std::make_pair(dgid, -1))); + gmg_grid.insert(std::make_pair(ref_loc, std::make_pair(dgid, -1))); return daughters.size() - 1; } diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 2b46493d6798..ddc69d465a03 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -109,8 +109,12 @@ class Tree : public std::enable_shared_from_this { int ndim; const std::uint64_t my_id; - // Location of block in this tree, current gid, previous gid - std::unordered_map> leaves; + // Structure mapping location of block in this tree to current gid and previous gid + using LocMap_t = std::unordered_map>; + LocMap_t leaves; + // Two-level composite grids for geometric multigrid + std::vector gmg_tlc_grids; + std::unordered_set internal_nodes; // This contains all of the neighbor information for this tree, for each of the From a5e19a6d4304d84136a0ed1f3afc674c894b5197 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 14:48:47 -0600 Subject: [PATCH 07/39] more work --- doc/sphinx/src/mesh/mesh.rst | 5 ++++- src/mesh/forest/tree.cpp | 4 ++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/doc/sphinx/src/mesh/mesh.rst b/doc/sphinx/src/mesh/mesh.rst index 88ffa29af2ed..c76896d1e2fe 100644 --- a/doc/sphinx/src/mesh/mesh.rst +++ b/doc/sphinx/src/mesh/mesh.rst @@ -71,7 +71,10 @@ level that are neighbors to finer blocks, which implies that below the root grid level the blocks may not cover the entire mesh. For levels above the root grid, blocks may change shape so that they only cover the domain of the root grid. Note that leaf blocks may be contained in multiple blocklists, and the lists all point -to the same block (not a separate copy). To be explicit, when ``parthenon/mesh/multigrid`` is set to ``true`` blocks corresponding to *all* internal nodes of the refinement tree are created, in addition to the leaf node blocks that are normally created. +to the same block (not a separate copy). To be explicit, when +``parthenon/mesh/multigrid`` is set to ``true`` blocks corresponding to *all* +internal nodes of the refinement tree are created, in addition to the leaf node blocks +that are normally created. *GMG Implementation Note:* The reason for including two levels in the GMG block lists is for dealing with diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index a3a376ce80fb..2f4671f713ed 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -50,13 +50,13 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz internal_nodes.emplace(my_id, l, i, j, k); } gmg_tlc_grids[l].emplace(std::make_pair(LogicalLocation(my_id, l, i, j, k), - std::make_pair(-1, -1))) + std::make_pair(-1, -1))); } } } } // Pre-populate the next finest two-level composite grid with coarser blocks - gmg_tlc_grid[l + 1] = gmg_tlc_grid[l]; + gmg_tlc_grids[root_level + 1] = gmg_tlc_grids[root_level]; } int Tree::AddMeshBlock(const LogicalLocation &loc, bool enforce_proper_nesting) { From 3345de8aa21091b86a24ab7a63ff66c920a5cfab Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 16:19:41 -0600 Subject: [PATCH 08/39] simplify level scheme for gmg --- src/driver/driver.cpp | 4 ++-- src/interface/data_collection.cpp | 15 ++++++------- src/interface/mesh_data.cpp | 6 ++--- src/mesh/mesh-gmg.cpp | 37 ++++++++++++++----------------- src/mesh/mesh.cpp | 8 +++---- src/mesh/mesh.hpp | 9 ++++---- src/solvers/mg_solver.hpp | 6 +++-- 7 files changed, 42 insertions(+), 43 deletions(-) diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp index f384a5503eda..372ba30f0964 100644 --- a/src/driver/driver.cpp +++ b/src/driver/driver.cpp @@ -186,8 +186,8 @@ void EvolutionDriver::InitializeBlockTimeStepsAndBoundaries() { auto &mbase = pmesh->mesh_data.GetOrAdd("base", i); Update::EstimateTimestep(mbase.get()); BuildBoundaryBuffers(mbase); - for (int gmg_level = 0; gmg_level < pmesh->gmg_mesh_data.size(); ++gmg_level) { - auto &mdg = pmesh->gmg_mesh_data[gmg_level].GetOrAdd(gmg_level, "base", i); + for (auto &[gmg_level, mdc] : pmesh->gmg_mesh_data) { + auto &mdg = mdc.GetOrAdd(gmg_level, "base", i); BuildBoundaryBuffers(mdg); BuildGMGBoundaryBuffers(mdg); } diff --git a/src/interface/data_collection.cpp b/src/interface/data_collection.cpp index 3ac429fd6df4..4a28e0bce993 100644 --- a/src/interface/data_collection.cpp +++ b/src/interface/data_collection.cpp @@ -58,9 +58,9 @@ std::shared_ptr> & GetOrAdd_impl(Mesh *pmy_mesh_, std::map>> &containers_, BlockList_t &block_list, const std::string &mbd_label, - const int &partition_id, const int gmg_level) { + const int &partition_id, const std::optional gmg_level) { std::string label = mbd_label + "_part-" + std::to_string(partition_id); - if (gmg_level >= 0) label = label + "_gmg-" + std::to_string(gmg_level); + if (gmg_level) label = label + "_gmg-" + std::to_string(*gmg_level); auto it = containers_.find(label); if (it == containers_.end()) { // TODO(someone) add caching of partitions to Mesh at some point @@ -70,13 +70,12 @@ GetOrAdd_impl(Mesh *pmy_mesh_, if (partitions.size() == 0) partitions = std::vector(1); for (auto i = 0; i < partitions.size(); i++) { std::string md_label = mbd_label + "_part-" + std::to_string(i); - if (gmg_level >= 0) md_label = md_label + "_gmg-" + std::to_string(gmg_level); + if (gmg_level) md_label = md_label + "_gmg-" + std::to_string(*gmg_level); containers_[md_label] = std::make_shared>(mbd_label); containers_[md_label]->Set(partitions[i], pmy_mesh_); - if (gmg_level >= 0) { - int min_gmg_logical_level = pmy_mesh_->GetGMGMinLogicalLevel(); - containers_[md_label]->grid = GridIdentifier{GridType::two_level_composite, - gmg_level + min_gmg_logical_level}; + if (gmg_level) { + containers_[md_label]->grid = + GridIdentifier{GridType::two_level_composite, *gmg_level}; } else { containers_[md_label]->grid = GridIdentifier{GridType::leaf, 0}; } @@ -90,7 +89,7 @@ std::shared_ptr> & DataCollection>::GetOrAdd(const std::string &mbd_label, const int &partition_id) { return GetOrAdd_impl(pmy_mesh_, containers_, pmy_mesh_->block_list, mbd_label, - partition_id, -1); + partition_id, {}); } template <> diff --git a/src/interface/mesh_data.cpp b/src/interface/mesh_data.cpp index 27a4dc520180..9f326c16e222 100644 --- a/src/interface/mesh_data.cpp +++ b/src/interface/mesh_data.cpp @@ -28,10 +28,10 @@ void MeshData::Initialize(const MeshData *src, grid = src->grid; if (grid.type == GridType::two_level_composite) { - int gmg_level = src->grid.logical_level - pmy_mesh_->GetGMGMinLogicalLevel(); for (int i = 0; i < nblocks; i++) { - block_data_[i] = pmy_mesh_->gmg_block_lists[gmg_level][i]->meshblock_data.Add( - stage_name_, src->GetBlockData(i), names, shallow); + block_data_[i] = + pmy_mesh_->gmg_block_lists[src->grid.logical_level][i]->meshblock_data.Add( + stage_name_, src->GetBlockData(i), names, shallow); } } else { for (int i = 0; i < nblocks; i++) { diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 429cac855543..55a59b4ebbe5 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -157,39 +157,39 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app const int gmg_min_level = root_level - gmg_level_offset; gmg_min_logical_level_ = gmg_min_level; - const int gmg_levels = current_level - gmg_min_level + 1; - gmg_grid_locs = std::vector(gmg_levels); - gmg_block_lists = std::vector(gmg_levels); + for (int level = gmg_min_level; level <= current_level; ++level) { + gmg_grid_locs[level] = LogicalLocMap_t(); + gmg_block_lists[level] = BlockList_t(); + gmg_mesh_data[level] = DataCollection>(); + } // Create MeshData objects for GMG - gmg_mesh_data = std::vector>>(gmg_levels); - for (auto &mdc : gmg_mesh_data) + for (auto &[l, mdc] : gmg_mesh_data) mdc.SetMeshPointer(this); // Add leaf grid locations to GMG grid levels int gmg_gid = 0; for (auto loc : loclist) { - const int gmg_level = gmg_levels - 1 + loc.level() - current_level; + const int gmg_level = loc.level(); gmg_grid_locs[gmg_level].insert( {loc, std::pair(gmg_gid, ranklist[gmg_gid])}); - if (gmg_level < gmg_levels - 1) { + if (gmg_level < current_level) { gmg_grid_locs[gmg_level + 1].insert( {loc, std::pair(gmg_gid, ranklist[gmg_gid])}); } if (ranklist[gmg_gid] == Globals::my_rank) { const int lid = gmg_gid - nslist[Globals::my_rank]; gmg_block_lists[gmg_level].push_back(block_list[lid]); - if (gmg_level < gmg_levels - 1) + if (gmg_level < current_level) gmg_block_lists[gmg_level + 1].push_back(block_list[lid]); } gmg_gid++; } // Fill in internal nodes for GMG grid levels from levels on finer GMG grid - for (int gmg_level = gmg_levels - 2; gmg_level >= 0; --gmg_level) { - int grid_logical_level = gmg_level - gmg_levels + 1 + current_level; + for (int gmg_level = current_level - 1; gmg_level >= gmg_min_level; --gmg_level) { for (auto &[loc, gid_rank] : gmg_grid_locs[gmg_level + 1]) { - if (loc.level() == grid_logical_level + 1) { + if (loc.level() == gmg_level + 1) { auto parent = loc.GetParent(); if (parent.morton() == loc.morton()) { gmg_grid_locs[gmg_level].insert( @@ -210,17 +210,15 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app // Find same level neighbors on all GMG levels auto root_grid = this->GetRootGridInfo(); - for (int gmg_level = 0; gmg_level < gmg_levels; ++gmg_level) { - int grid_logical_level = gmg_level - gmg_levels + 1 + current_level; + for (int gmg_level = gmg_min_level; gmg_level <= current_level; ++gmg_level) { SetSameLevelNeighbors(gmg_block_lists[gmg_level], gmg_grid_locs[gmg_level], root_grid, - nbs, true, grid_logical_level); + nbs, true, gmg_level); } // Now find GMG coarser neighbor - for (int gmg_level = 1; gmg_level < gmg_levels; ++gmg_level) { - int grid_logical_level = gmg_level - gmg_levels + 1 + current_level; + for (int gmg_level = gmg_min_level + 1; gmg_level <= current_level; ++gmg_level) { for (auto &pmb : gmg_block_lists[gmg_level]) { - if (pmb->loc.level() != grid_logical_level) continue; + if (pmb->loc.level() != gmg_level) continue; auto parent_loc = pmb->loc.GetParent(); auto loc = pmb->loc; auto gid = pmb->gid; @@ -239,10 +237,9 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } // Now find finer GMG neighbors - for (int gmg_level = 0; gmg_level < gmg_levels - 1; ++gmg_level) { - int grid_logical_level = gmg_level - gmg_levels + 1 + current_level; + for (int gmg_level = gmg_min_level; gmg_level <= current_level - 1; ++gmg_level) { for (auto &pmb : gmg_block_lists[gmg_level]) { - if (pmb->loc.level() != grid_logical_level) continue; + if (pmb->loc.level() != gmg_level) continue; auto daughter_locs = pmb->loc.GetDaughters(); for (auto &daughter_loc : daughter_locs) { if (gmg_grid_locs[gmg_level + 1].count(daughter_loc) > 0) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index cee6af316fc8..1bb6b2e27d1d 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -997,8 +997,8 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * for (int i = 0; i < num_partitions; i++) { auto &md = mesh_data.GetOrAdd("base", i); tag_map.AddMeshDataToMap(md); - for (int gmg_level = 0; gmg_level < gmg_mesh_data.size(); ++gmg_level) { - auto &mdg = gmg_mesh_data[gmg_level].GetOrAdd(gmg_level, "base", i); + for (auto &[gmg_level, mdc] : gmg_mesh_data) { + auto &mdg = mdc.GetOrAdd(gmg_level, "base", i); // tag_map.AddMeshDataToMap(mdg); tag_map.AddMeshDataToMap(mdg); tag_map.AddMeshDataToMap(mdg); @@ -1036,8 +1036,8 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput * for (int i = 0; i < num_partitions; i++) { auto &md = mesh_data.GetOrAdd("base", i); BuildBoundaryBuffers(md); - for (int gmg_level = 0; gmg_level < gmg_mesh_data.size(); ++gmg_level) { - auto &mdg = gmg_mesh_data[gmg_level].GetOrAdd(gmg_level, "base", i); + for (auto &[gmg_level, mdc] : gmg_mesh_data) { + auto &mdg = mdc.GetOrAdd(gmg_level, "base", i); BuildBoundaryBuffers(mdg); BuildGMGBoundaryBuffers(mdg); } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 6963b30b046f..9aaf9e894e87 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -121,10 +121,11 @@ class Mesh { DataCollection> mesh_data; LogicalLocMap_t leaf_grid_locs; - std::vector gmg_grid_locs; - std::vector gmg_block_lists; - std::vector>> gmg_mesh_data; - int GetGMGMaxLevel() { return gmg_grid_locs.size() - 1; } + std::map gmg_grid_locs; + std::map gmg_block_lists; + std::map>> gmg_mesh_data; + int GetGMGMaxLevel() { return current_level; } + int GetGMGMinLevel() { return gmg_min_logical_level_; } int GetGMGMinLogicalLevel() { return gmg_min_logical_level_; } // functions diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp index 153220a4a4cf..55ac23f43eb7 100644 --- a/src/solvers/mg_solver.hpp +++ b/src/solvers/mg_solver.hpp @@ -142,7 +142,8 @@ class MGSolver { using namespace utils; iter_counter = 0; - int min_level = std::max(pmesh->GetGMGMaxLevel() - params_.max_coarsenings, 0); + int min_level = std::max(pmesh->GetGMGMaxLevel() - params_.max_coarsenings, + pmesh->GetGMGMinLevel()); int max_level = pmesh->GetGMGMaxLevel(); return AddMultiGridTasksPartitionLevel(tl, dependence, partition, max_level, @@ -153,7 +154,8 @@ class MGSolver { TaskID AddSetupTasks(TL_t &tl, TaskID dependence, int partition, Mesh *pmesh) { using namespace utils; - int min_level = std::max(pmesh->GetGMGMaxLevel() - params_.max_coarsenings, 0); + int min_level = std::max(pmesh->GetGMGMaxLevel() - params_.max_coarsenings, + pmesh->GetGMGMinLevel()); int max_level = pmesh->GetGMGMaxLevel(); return AddMultiGridSetupPartitionLevel(tl, dependence, partition, max_level, From 3a4da84b2870838dc4036cf94475282312c6f80f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 17:44:21 -0600 Subject: [PATCH 09/39] provide gids for the whole mesh --- src/mesh/forest/forest.cpp | 11 +++++++- src/mesh/forest/tree.cpp | 33 +++++++++++++++--------- src/mesh/forest/tree.hpp | 51 +++++++++++++++++++++++++++++++------- 3 files changed, 73 insertions(+), 22 deletions(-) diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index 0b4f191450da..4df3fd84f93e 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -41,13 +41,22 @@ std::vector Forest::GetMeshBlockListAndResolveGids() { std::uint64_t gid{0}; for (auto &[id, tree] : trees) { std::size_t start = mb_list.size(); - auto tree_mbs = tree->GetMeshBlockList(); + auto tree_mbs = tree->GetSortedMeshBlockList(); mb_list.insert(mb_list.end(), std::make_move_iterator(tree_mbs.begin()), std::make_move_iterator(tree_mbs.end())); std::size_t end = mb_list.size(); for (int i = start; i < end; ++i) tree->InsertGid(mb_list[i], gid++); } + + // Assign gids to the internal nodes + for (auto &[id, tree] : trees) { + std::size_t start = mb_list.size(); + auto tree_int_locs = tree->GetSortedInternalNodeList(); + for (auto &loc : tree_int_locs) + tree->InsertGid(loc, gid++); + } + // The index of blocks in this list corresponds to their gid gids_resolved = true; return mb_list; diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index 2f4671f713ed..3b7f406bd3af 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -43,14 +43,13 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz for (int k = 0; k < (ndim > 2 ? (1LL << l) : 1); ++k) { for (int j = 0; j < (ndim > 1 ? (1LL << l) : 1); ++j) { for (int i = 0; i < (ndim > 0 ? (1LL << l) : 1); ++i) { + LogicalLocation loc(my_id, l, i, j, k); if (l == root_level) { - leaves.emplace(std::make_pair(LogicalLocation(my_id, l, i, j, k), - std::make_pair(-1, -1))); + leaves.emplace(LocMapEntry(loc, -1, -1)); } else { - internal_nodes.emplace(my_id, l, i, j, k); + internal_nodes.emplace(LocMapEntry(loc, -1, -1)); } - gmg_tlc_grids[l].emplace(std::make_pair(LogicalLocation(my_id, l, i, j, k), - std::make_pair(-1, -1))); + gmg_tlc_grids[l].emplace(LocMapEntry(loc, -1, -1)); } } } @@ -94,9 +93,9 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { std::vector daughters = ref_loc.GetDaughters(ndim); auto gid_parent = leaves[ref_loc].first; leaves.erase(ref_loc); - internal_nodes.insert(ref_loc); + internal_nodes.insert(LocMapEntry(ref_loc, -1, -1)); for (auto &d : daughters) { - leaves.insert(std::make_pair(d, std::make_pair(gid_parent, -1))); + leaves.insert(LocMapEntry(d, gid_parent, -1)); } int nadded = daughters.size() - 1; @@ -106,9 +105,9 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { if (gmg_grid_1.count(ref_loc)) gmg_grid_1.erase(ref_loc); for (auto &d : daughters) { // Insert as fine leaf blocks on coarser two-level grid - gmg_grid_1.insert(std::make_pair(d, std::make_pair(gid_parent, -1))); + gmg_grid_1.insert(LocMapEntry(d, gid_parent, -1)); // Insert as coarse leaf blocks on finer two-level grid - gmg_grid_2.insert(std::make_pair(d, std::make_pair(gid_parent, -1))); + gmg_grid_2.insert(LocMapEntry(d, gid_parent, -1)); } if (enforce_proper_nesting) { @@ -241,12 +240,12 @@ int Tree::Derefine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) dgid = std::min(dgid, node.mapped().first); } internal_nodes.erase(ref_loc); - leaves.insert(std::make_pair(ref_loc, std::make_pair(dgid, -1))); - gmg_grid.insert(std::make_pair(ref_loc, std::make_pair(dgid, -1))); + leaves.insert(LocMapEntry(ref_loc, dgid, -1)); + gmg_grid.insert(LocMapEntry(ref_loc, dgid, -1)); return daughters.size() - 1; } -std::vector Tree::GetMeshBlockList() const { +std::vector Tree::GetSortedMeshBlockList() const { std::vector mb_list; mb_list.reserve(leaves.size()); for (auto &[loc, gid] : leaves) @@ -256,6 +255,16 @@ std::vector Tree::GetMeshBlockList() const { return mb_list; } +std::vector Tree::GetSortedInternalNodeList() const { + std::vector mb_list; + mb_list.reserve(internal_nodes.size()); + for (auto &[loc, gid] : internal_nodes) + mb_list.push_back(loc); + std::sort(mb_list.begin(), mb_list.end(), + [](const auto &a, const auto &b) { return a < b; }); + return mb_list; +} + RegionSize Tree::GetBlockDomain(const LogicalLocation &loc) const { PARTHENON_REQUIRE(loc.IsInTree(), "Probably there is a mistake..."); RegionSize out = domain; diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index ddc69d465a03..600b4e56e784 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -65,7 +65,8 @@ class Tree : public std::enable_shared_from_this { // Methods for getting block properties int count(const LogicalLocation &loc) const { return leaves.count(loc); } - std::vector GetMeshBlockList() const; + std::vector GetSortedMeshBlockList() const; + std::vector GetSortedInternalNodeList() const; RegionSize GetBlockDomain(const LogicalLocation &loc) const; std::array GetBlockBCs(const LogicalLocation &loc) const; std::vector FindNeighbors(const LogicalLocation &loc) const; @@ -89,15 +90,41 @@ class Tree : public std::enable_shared_from_this { } void InsertGid(const LogicalLocation &loc, std::int64_t gid) { - PARTHENON_REQUIRE(leaves.count(loc) == 1, - "Trying to add gid for non-existent location."); - leaves[loc].second = leaves[loc].first; - leaves[loc].first = gid; + if (leaves.count(loc)) { + leaves[loc].second = leaves[loc].first; + leaves[loc].first = gid; + } else if (internal_nodes.count(loc)) { + internal_nodes[loc].second = internal_nodes[loc].first; + internal_nodes[loc].first = gid; + } + PARTHENON_FAIL("Tried to assign gid to non-existent block."); } - std::int64_t GetGid(const LogicalLocation &loc) const { return leaves.at(loc).first; } - std::int64_t GetOldGid(const LogicalLocation &loc) const { - return leaves.at(loc).second; + std::int64_t GetGid(const LogicalLocation &loc) const { + if (leaves.count(loc)) { + return leaves.at(loc).first; + } else if (internal_nodes.count(loc)) { + return internal_nodes.at(loc).first; + } + PARTHENON_FAIL("Asking for GID of non-existent location."); + } + + std::int64_t GetLeafGid(const LogicalLocation &loc) const { + if (leaves.count(loc)) { + return leaves.at(loc).first; + } else if (internal_nodes.count(loc)) { + return GetLeafGid(loc.GetDaughter(0, 0, 0)); + } + PARTHENON_FAIL("Asking for GID of non-existent location."); + } + + std::int64_t GetOldGid(const LogicalLocation &loc) const { + if (leaves.count(loc)) { + return leaves.at(loc).second; + } else if (internal_nodes.count(loc)) { + return internal_nodes.at(loc).second; + } + PARTHENON_FAIL("Asking for GID of non-existent location."); } // TODO(LFR): Eventually remove this. @@ -111,11 +138,17 @@ class Tree : public std::enable_shared_from_this { const std::uint64_t my_id; // Structure mapping location of block in this tree to current gid and previous gid using LocMap_t = std::unordered_map>; + static std::pair> + LocMapEntry(const LogicalLocation& loc, const int gid, const int gid_old) { + return std::make_pair(loc, std::make_pair(gid, gid_old)); + } + LocMap_t leaves; + LocMap_t internal_nodes; + // Two-level composite grids for geometric multigrid std::vector gmg_tlc_grids; - std::unordered_set internal_nodes; // This contains all of the neighbor information for this tree, for each of the // 3^3 possible neighbor connections. Since an edge or node connection can have From 102b012730a6f689739e6c33fc39858fdf42acd1 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 18:20:06 -0600 Subject: [PATCH 10/39] bugfix --- src/mesh/forest/tree.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 600b4e56e784..64b136e178a9 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -96,8 +96,9 @@ class Tree : public std::enable_shared_from_this { } else if (internal_nodes.count(loc)) { internal_nodes[loc].second = internal_nodes[loc].first; internal_nodes[loc].first = gid; + } else { + PARTHENON_FAIL("Tried to assign gid to non-existent block."); } - PARTHENON_FAIL("Tried to assign gid to non-existent block."); } std::int64_t GetGid(const LogicalLocation &loc) const { From e3c6e81f71d42aa29b7bfa2fbb1ffcaff3c205aa Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 18:20:59 -0600 Subject: [PATCH 11/39] format --- src/mesh/forest/forest.cpp | 4 ++-- src/mesh/forest/tree.cpp | 10 +++++----- src/mesh/forest/tree.hpp | 28 ++++++++++++++-------------- src/mesh/mesh.cpp | 2 +- 4 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/mesh/forest/forest.cpp b/src/mesh/forest/forest.cpp index 4df3fd84f93e..1a7c39f3441d 100644 --- a/src/mesh/forest/forest.cpp +++ b/src/mesh/forest/forest.cpp @@ -48,8 +48,8 @@ std::vector Forest::GetMeshBlockListAndResolveGids() { for (int i = start; i < end; ++i) tree->InsertGid(mb_list[i], gid++); } - - // Assign gids to the internal nodes + + // Assign gids to the internal nodes for (auto &[id, tree] : trees) { std::size_t start = mb_list.size(); auto tree_int_locs = tree->GetSortedInternalNodeList(); diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index 3b7f406bd3af..a65d12a980ca 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -98,11 +98,11 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { leaves.insert(LocMapEntry(d, gid_parent, -1)); } int nadded = daughters.size() - 1; - + // Update multigrid levels - auto& gmg_grid_1 = gmg_tlc_grids[ref_loc.level() + 1]; - auto& gmg_grid_2 = gmg_tlc_grids[ref_loc.level() + 2]; - if (gmg_grid_1.count(ref_loc)) gmg_grid_1.erase(ref_loc); + auto &gmg_grid_1 = gmg_tlc_grids[ref_loc.level() + 1]; + auto &gmg_grid_2 = gmg_tlc_grids[ref_loc.level() + 2]; + if (gmg_grid_1.count(ref_loc)) gmg_grid_1.erase(ref_loc); for (auto &d : daughters) { // Insert as fine leaf blocks on coarser two-level grid gmg_grid_1.insert(LocMapEntry(d, gid_parent, -1)); @@ -232,7 +232,7 @@ int Tree::Derefine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) } // Derefinement is ok - auto& gmg_grid = gmg_tlc_grids[ref_loc.level() + 1]; + auto &gmg_grid = gmg_tlc_grids[ref_loc.level() + 1]; std::int64_t dgid = std::numeric_limits::max(); for (auto &d : daughters) { gmg_grid.erase(d); diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 64b136e178a9..ef07060b2aa8 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -93,38 +93,38 @@ class Tree : public std::enable_shared_from_this { if (leaves.count(loc)) { leaves[loc].second = leaves[loc].first; leaves[loc].first = gid; - } else if (internal_nodes.count(loc)) { - internal_nodes[loc].second = internal_nodes[loc].first; + } else if (internal_nodes.count(loc)) { + internal_nodes[loc].second = internal_nodes[loc].first; internal_nodes[loc].first = gid; } else { PARTHENON_FAIL("Tried to assign gid to non-existent block."); } } - std::int64_t GetGid(const LogicalLocation &loc) const { + std::int64_t GetGid(const LogicalLocation &loc) const { if (leaves.count(loc)) { return leaves.at(loc).first; } else if (internal_nodes.count(loc)) { return internal_nodes.at(loc).first; - } + } PARTHENON_FAIL("Asking for GID of non-existent location."); } - std::int64_t GetLeafGid(const LogicalLocation &loc) const { + std::int64_t GetLeafGid(const LogicalLocation &loc) const { if (leaves.count(loc)) { return leaves.at(loc).first; } else if (internal_nodes.count(loc)) { return GetLeafGid(loc.GetDaughter(0, 0, 0)); - } + } PARTHENON_FAIL("Asking for GID of non-existent location."); } - std::int64_t GetOldGid(const LogicalLocation &loc) const { + std::int64_t GetOldGid(const LogicalLocation &loc) const { if (leaves.count(loc)) { return leaves.at(loc).second; } else if (internal_nodes.count(loc)) { return internal_nodes.at(loc).second; - } + } PARTHENON_FAIL("Asking for GID of non-existent location."); } @@ -138,18 +138,18 @@ class Tree : public std::enable_shared_from_this { int ndim; const std::uint64_t my_id; // Structure mapping location of block in this tree to current gid and previous gid - using LocMap_t = std::unordered_map>; - static std::pair> - LocMapEntry(const LogicalLocation& loc, const int gid, const int gid_old) { + using LocMap_t = + std::unordered_map>; + static std::pair> + LocMapEntry(const LogicalLocation &loc, const int gid, const int gid_old) { return std::make_pair(loc, std::make_pair(gid, gid_old)); } LocMap_t leaves; LocMap_t internal_nodes; - - // Two-level composite grids for geometric multigrid - std::vector gmg_tlc_grids; + // Two-level composite grids for geometric multigrid + std::vector gmg_tlc_grids; // This contains all of the neighbor information for this tree, for each of the // 3^3 possible neighbor connections. Since an edge or node connection can have diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 8570ec5f8912..401e134e4efe 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -918,7 +918,7 @@ void Mesh::BuildTagMapAndBoundaryBuffers() { auto &md = mesh_data.GetOrAdd("base", i); tag_map.AddMeshDataToMap(md); for (auto &[gmg_level, mdc] : gmg_mesh_data) { - auto &mdg = mdc.GetOrAdd(gmg_level, "base", i); + auto &mdg = mdc.GetOrAdd(gmg_level, "base", i); tag_map.AddMeshDataToMap(mdg); tag_map.AddMeshDataToMap(mdg); tag_map.AddMeshDataToMap(mdg); From 71c67920c04b092f6056329d75fa9ff795fa9eed Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 19:18:55 -0600 Subject: [PATCH 12/39] Start on refactor of building gmg block lists --- src/mesh/forest/forest.hpp | 4 +++ src/mesh/forest/tree.hpp | 6 ++--- src/mesh/mesh-gmg.cpp | 55 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 60 insertions(+), 5 deletions(-) diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index 04bb60ae2494..b1c6aa483d1b 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -133,6 +133,10 @@ class Forest { PARTHENON_REQUIRE(gids_resolved, "Asking for GID in invalid state."); return trees.at(loc.tree())->GetGid(loc); } + std::int64_t GetLeafGid(const LogicalLocation &loc) const { + PARTHENON_REQUIRE(gids_resolved, "Asking for GID in invalid state."); + return trees.at(loc.tree())->GetLeafGid(loc); + } std::int64_t GetOldGid(const LogicalLocation &loc) const { PARTHENON_REQUIRE(gids_resolved, "Asking for GID in invalid state."); diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index ef07060b2aa8..520a7baa9fd0 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -107,7 +107,7 @@ class Tree : public std::enable_shared_from_this { } else if (internal_nodes.count(loc)) { return internal_nodes.at(loc).first; } - PARTHENON_FAIL("Asking for GID of non-existent location."); + return -1; } std::int64_t GetLeafGid(const LogicalLocation &loc) const { @@ -116,7 +116,7 @@ class Tree : public std::enable_shared_from_this { } else if (internal_nodes.count(loc)) { return GetLeafGid(loc.GetDaughter(0, 0, 0)); } - PARTHENON_FAIL("Asking for GID of non-existent location."); + return -1; } std::int64_t GetOldGid(const LogicalLocation &loc) const { @@ -125,7 +125,7 @@ class Tree : public std::enable_shared_from_this { } else if (internal_nodes.count(loc)) { return internal_nodes.at(loc).second; } - PARTHENON_FAIL("Asking for GID of non-existent location."); + return -1; } // TODO(LFR): Eventually remove this. diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 961912e3187e..d4489d7da33d 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -221,7 +221,57 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app // Create MeshData objects for GMG for (auto &[l, mdc] : gmg_mesh_data) mdc.SetMeshPointer(this); + + // Fill/create gmg block lists based on this ranks block list + for (auto &pmb : block_list) { + const int level = pmb->loc.level(); + // Add the leaf block to its level + gmg_block_lists[level].push_back(pmb); + + // Add the leaf block to the next finer level if required + if (level < current_level) { + gmg_block_lists[level + 1].push_back(pmb); + } + + // Create internal blocks that share a Morton number with this block + // and add them to gmg two-level composite grid block lists + auto loc = pmb->loc.GetParent(); + while (loc.level() <= gmg_min_level && loc.morton() == pmb->loc.morton()) { + RegionSize block_size = GetBlockSize(); + BoundaryFlag block_bcs[6]; + SetBlockSizeAndBoundaries(loc, block_size, block_bcs); + gmg_block_lists[loc.level()].push_back( + MeshBlock::Make(forest.GetGid(loc), -1, loc, block_size, block_bcs, this, pin, + app_in, packages, resolved_packages, gflag)); + loc = loc.GetParent(); + } + } + // Sort the gmg block lists by gid and find coarser and finer neighbors + for (auto &[level, bl] : gmg_block_lists) { + std::sort(bl.begin(), bl.end(), [](auto &a, auto &b){ return a->gid < b->gid;}); + for (auto &pmb : bl) { + auto ploc = pmb->loc.GetParent(); + int gid = forest.GetGid(ploc); + if (gid >= 0) { + int leaf_gid = forest.GetLeafGid(ploc); + pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, ploc, ranklist[leaf_gid], gid, + std::array{0, 0, 0}, 0, 0, 0, 0); + } + + auto dlocs = pmb->loc.GetDaughters(ndim); + for (auto &d : dlocs) { + int gid = forest.GetGid(d); + if (gid >= 0) { + int leaf_gid = forest.GetLeafGid(d); + pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], gid, + std::array{0, 0, 0}, 0, 0, 0, 0); + } + } + } + } + +/* // Add leaf grid locations to GMG grid levels int gmg_gid = 0; for (auto loc : loclist) { @@ -262,14 +312,14 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } } } - +*/ // Find same level neighbors on all GMG levels auto root_grid = this->GetRootGridInfo(); for (int gmg_level = gmg_min_level; gmg_level <= current_level; ++gmg_level) { SetSameLevelNeighbors(gmg_block_lists[gmg_level], gmg_grid_locs[gmg_level], root_grid, nbs, true, gmg_level); } - + /* // Now find GMG coarser neighbor for (int gmg_level = gmg_min_level + 1; gmg_level <= current_level; ++gmg_level) { for (auto &pmb : gmg_block_lists[gmg_level]) { @@ -307,5 +357,6 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } } } + */ } } // namespace parthenon From bd733161fe71e26d399cd1991208abe9a9233b46 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 26 Mar 2024 20:18:59 -0600 Subject: [PATCH 13/39] generalize neighbor finding to two-level composite grids --- src/basic_types.hpp | 9 +++++ src/interface/data_collection.cpp | 5 +-- src/interface/mesh_data.hpp | 6 --- src/mesh/forest/forest.hpp | 4 +- src/mesh/forest/tree.cpp | 48 +++++++++++++++++++----- src/mesh/forest/tree.hpp | 5 ++- src/mesh/mesh-amr_loadbalance.cpp | 4 +- src/mesh/mesh-gmg.cpp | 61 +++++++++++++++++++------------ src/mesh/mesh.cpp | 4 +- src/mesh/mesh.hpp | 2 +- 10 files changed, 97 insertions(+), 51 deletions(-) diff --git a/src/basic_types.hpp b/src/basic_types.hpp index b36383aa2827..1557c07f981b 100644 --- a/src/basic_types.hpp +++ b/src/basic_types.hpp @@ -66,6 +66,15 @@ enum class BoundaryType : int { gmg_prolongate_recv }; +enum class GridType { none, leaf, two_level_composite, single_level_with_internal }; +struct GridIdentifier { + GridType type = GridType::none; + int logical_level = 0; + + static GridIdentifier leaf() {return GridIdentifier{GridType::leaf, 0};} + static GridIdentifier two_level_composite(int level) {return GridIdentifier{GridType::two_level_composite, level};} +}; + constexpr bool IsSender(BoundaryType btype) { if (btype == BoundaryType::flxcor_recv) return false; if (btype == BoundaryType::gmg_restrict_recv) return false; diff --git a/src/interface/data_collection.cpp b/src/interface/data_collection.cpp index 4a28e0bce993..589c2c7128a6 100644 --- a/src/interface/data_collection.cpp +++ b/src/interface/data_collection.cpp @@ -74,10 +74,9 @@ GetOrAdd_impl(Mesh *pmy_mesh_, containers_[md_label] = std::make_shared>(mbd_label); containers_[md_label]->Set(partitions[i], pmy_mesh_); if (gmg_level) { - containers_[md_label]->grid = - GridIdentifier{GridType::two_level_composite, *gmg_level}; + containers_[md_label]->grid = GridIdentifier::two_level_composite(*gmg_level); } else { - containers_[md_label]->grid = GridIdentifier{GridType::leaf, 0}; + containers_[md_label]->grid = GridIdentifier::leaf(); } } } diff --git a/src/interface/mesh_data.hpp b/src/interface/mesh_data.hpp index 07726e864f56..4df20bf26d88 100644 --- a/src/interface/mesh_data.hpp +++ b/src/interface/mesh_data.hpp @@ -181,12 +181,6 @@ const MeshBlockPack

&PackOnMesh(M &map, BlockDataList_t &block_data_, } // namespace pack_on_mesh_impl -enum class GridType { none, leaf, two_level_composite, single_level_with_internal }; -struct GridIdentifier { - GridType type = GridType::none; - int logical_level = 0; -}; - /// The MeshData class is a container for cached MeshBlockPacks, i.e., it /// contains both the pointers to the MeshBlockData of the MeshBlocks contained /// in the object as well as maps to the cached MeshBlockPacks of VariablePacks or diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index b1c6aa483d1b..52a5f4763d57 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -84,8 +84,8 @@ class Forest { return trees.at(loc.tree())->FindNeighbors(loc, ox1, ox2, ox3); } - std::vector FindNeighbors(const LogicalLocation &loc) const { - return trees.at(loc.tree())->FindNeighbors(loc); + std::vector FindNeighbors(const LogicalLocation &loc, GridIdentifier grid_id = GridIdentifier::leaf()) const { + return trees.at(loc.tree())->FindNeighbors(loc, grid_id); } std::size_t CountMeshBlock() const { std::size_t count{0}; diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index a65d12a980ca..84d1c54a2365 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -136,7 +136,7 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { return nadded; } -std::vector Tree::FindNeighbors(const LogicalLocation &loc) const { +std::vector Tree::FindNeighbors(const LogicalLocation &loc, GridIdentifier grid_id) const { const Indexer3D offsets({ndim > 0 ? -1 : 0, ndim > 0 ? 1 : 0}, {ndim > 1 ? -1 : 0, ndim > 1 ? 1 : 0}, {ndim > 2 ? -1 : 0, ndim > 2 ? 1 : 0}); @@ -144,7 +144,7 @@ std::vector Tree::FindNeighbors(const LogicalLocation &loc) co for (int o = 0; o < offsets.size(); ++o) { auto [ox1, ox2, ox3] = offsets(o); if (std::abs(ox1) + std::abs(ox2) + std::abs(ox3) == 0) continue; - FindNeighborsImpl(loc, ox1, ox2, ox3, &neighbor_locs); + FindNeighborsImpl(loc, ox1, ox2, ox3, &neighbor_locs, grid_id); } const int clev = loc.level() - 1; @@ -155,32 +155,60 @@ std::vector Tree::FindNeighbors(const LogicalLocation &loc) co std::vector Tree::FindNeighbors(const LogicalLocation &loc, int ox1, int ox2, int ox3) const { std::vector neighbor_locs; - FindNeighborsImpl(loc, ox1, ox2, ox3, &neighbor_locs); + FindNeighborsImpl(loc, ox1, ox2, ox3, &neighbor_locs, GridIdentifier::leaf()); return neighbor_locs; } void Tree::FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int ox3, - std::vector *neighbor_locs) const { + std::vector *neighbor_locs, + GridIdentifier grid_id) const { PARTHENON_REQUIRE( loc.tree() == my_id, "Trying to find neighbors in a tree with a LogicalLocation on a different tree."); PARTHENON_REQUIRE(leaves.count(loc) == 1, "Location must be a leaf to find neighbors."); auto neigh = loc.GetSameLevelNeighbor(ox1, ox2, ox3); int n_idx = neigh.NeighborTreeIndex(); + + bool include_same, include_fine, include_internal, include_coarse; + if (grid_id.type == GridType::leaf) { + include_same = true; + include_fine = true; + include_internal = false; + include_coarse = true; + } else if (grid_id.type == GridType::two_level_composite) { + if (loc.level() == grid_id.logical_level) { + include_same = true; + include_fine = false; + include_internal = true; + include_coarse = true; + } else if (loc.level() == grid_id.logical_level - 1) { + include_same = false; + include_fine = true; + include_internal = false; + include_coarse = false; + } else { + PARTHENON_FAIL("Logic is wrong somewhere."); + } + } + for (auto &[neighbor_tree, orientation] : neighbors[n_idx]) { auto tneigh = orientation.Transform(neigh, neighbor_tree->GetId()); auto tloc = orientation.Transform(loc, neighbor_tree->GetId()); PARTHENON_REQUIRE(orientation.TransformBack(tloc, GetId()) == loc, "Inverse transform not working."); - if (neighbor_tree->leaves.count(tneigh)) { + if (neighbor_tree->leaves.count(tneigh) && include_same) { neighbor_locs->push_back({tneigh, orientation.TransformBack(tneigh, GetId())}); } else if (neighbor_tree->internal_nodes.count(tneigh)) { - auto daughters = tneigh.GetDaughters(neighbor_tree->ndim); - for (auto &n : daughters) { - if (tloc.IsNeighborForest(n)) - neighbor_locs->push_back({n, orientation.TransformBack(n, GetId())}); + if (include_fine) { + auto daughters = tneigh.GetDaughters(neighbor_tree->ndim); + for (auto &n : daughters) { + if (tloc.IsNeighborForest(n)) + neighbor_locs->push_back({n, orientation.TransformBack(n, GetId())}); + } + } else if (include_internal) { + neighbor_locs->push_back({tneigh, orientation.TransformBack(tneigh, GetId())}); } - } else if (neighbor_tree->leaves.count(tneigh.GetParent())) { + } else if (neighbor_tree->leaves.count(tneigh.GetParent()) && include_coarse) { auto neighp = orientation.TransformBack(tneigh.GetParent(), GetId()); // Since coarser neighbors can cover multiple elements of the origin block and // because our communication algorithm packs this extra data by hand, we do not wish diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 520a7baa9fd0..ea5b3c286cb2 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -69,7 +69,7 @@ class Tree : public std::enable_shared_from_this { std::vector GetSortedInternalNodeList() const; RegionSize GetBlockDomain(const LogicalLocation &loc) const; std::array GetBlockBCs(const LogicalLocation &loc) const; - std::vector FindNeighbors(const LogicalLocation &loc) const; + std::vector FindNeighbors(const LogicalLocation &loc, GridIdentifier grid_id = GridIdentifier::leaf()) const; std::vector FindNeighbors(const LogicalLocation &loc, int ox1, int ox2, int ox3) const; @@ -133,7 +133,8 @@ class Tree : public std::enable_shared_from_this { private: void FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int ox3, - std::vector *neighbor_locs) const; + std::vector *neighbor_locs, + GridIdentifier grid_type) const; int ndim; const std::uint64_t my_id; diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index 6359dbee3e89..32d98a0d9231 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -954,7 +954,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput // in order to maintain a consistent global state. // Thus we rebuild and synchronize the mesh now, but using a unique // neighbor precedence favoring the "old" fine blocks over "new" ones - SetMeshBlockNeighbors(block_list, nbs, ranklist, newly_refined); + SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist, newly_refined); BuildTagMapAndBoundaryBuffers(); std::string noncc = "mesh_internal_noncc"; for (int i = 0; i < DefaultNumPartitions(); ++i) { @@ -973,7 +973,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput // Rebuild just the ownership model, this time weighting the "new" fine blocks just // like any other blocks at their level. - SetMeshBlockNeighbors(block_list, nbs, ranklist); + SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); BuildGMGHierarchy(nbs, pin, app_in); // Ownership does not impact anything about the buffers, so we don't need to // rebuild them if they were built above diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index d4489d7da33d..561dee7f948c 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -51,8 +51,8 @@ void Mesh::PopulateLeafLocationMap() { } } -void Mesh::SetMeshBlockNeighbors( - BlockList_t &block_list, int nbs, const std::vector &ranklist, +void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, + BlockList_t &block_list, const std::vector &ranklist, const std::unordered_set &newly_refined) { Indexer3D offsets({ndim > 0 ? -1 : 0, ndim > 0 ? 1 : 0}, {ndim > 1 ? -1 : 0, ndim > 1 ? 1 : 0}, @@ -62,7 +62,7 @@ void Mesh::SetMeshBlockNeighbors( for (auto &pmb : block_list) { std::vector all_neighbors; const auto &loc = pmb->loc; - auto neighbors = forest.FindNeighbors(loc); + auto neighbors = forest.FindNeighbors(loc, grid_id); // Build NeighborBlocks for unique neighbors int buf_id = 0; @@ -79,8 +79,8 @@ void Mesh::SetMeshBlockNeighbors( auto fn = nloc.origin_loc.GetAthenaXXFaceOffsets(loc, -offsets[0], -offsets[1], -offsets[2]); int tid = buffer_id.GetID(-offsets[0], -offsets[1], -offsets[2], fn[0], fn[1]); - - all_neighbors.emplace_back(pmb->pmy_mesh, nloc.global_loc, ranklist[gid], gid, + int lgid = forest.GetLeafGid(nloc.global_loc); + all_neighbors.emplace_back(pmb->pmy_mesh, nloc.global_loc, ranklist[lgid], gid, offsets, bid, tid, f[0], f[1]); // Set neighbor block ownership @@ -91,8 +91,16 @@ void Mesh::SetMeshBlockNeighbors( DetermineOwnershipForest(nloc.global_loc, neighbor_neighbors, newly_refined); nb.ownership.initialized = true; } - - pmb->neighbors = all_neighbors; + + if (grid_id.type == GridType::leaf) { + pmb->neighbors = all_neighbors; + } else if (grid_id.type == GridType::two_level_composite + && pmb->loc.level() == grid_id.logical_level) { + pmb->gmg_same_neighbors = all_neighbors; + } else if (grid_id.type == GridType::two_level_composite + && pmb->loc.level() == grid_id.logical_level - 1) { + pmb->gmg_composite_finer_neighbors = all_neighbors; + } } } @@ -247,27 +255,34 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } } - // Sort the gmg block lists by gid and find coarser and finer neighbors + // Sort the gmg block lists by gid and find neighbors for (auto &[level, bl] : gmg_block_lists) { std::sort(bl.begin(), bl.end(), [](auto &a, auto &b){ return a->gid < b->gid;}); - for (auto &pmb : bl) { - auto ploc = pmb->loc.GetParent(); - int gid = forest.GetGid(ploc); - if (gid >= 0) { - int leaf_gid = forest.GetLeafGid(ploc); - pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, ploc, ranklist[leaf_gid], gid, - std::array{0, 0, 0}, 0, 0, 0, 0); - } - auto dlocs = pmb->loc.GetDaughters(ndim); - for (auto &d : dlocs) { - int gid = forest.GetGid(d); + for (auto &pmb : bl) { + if (pmb->loc.level() > gmg_min_level) { + auto ploc = pmb->loc.GetParent(); + int gid = forest.GetGid(ploc); if (gid >= 0) { - int leaf_gid = forest.GetLeafGid(d); - pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], gid, + int leaf_gid = forest.GetLeafGid(ploc); + pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, ploc, ranklist[leaf_gid], gid, std::array{0, 0, 0}, 0, 0, 0, 0); } } + + if (pmb->loc.level() < current_level) { + auto dlocs = pmb->loc.GetDaughters(ndim); + for (auto &d : dlocs) { + int gid = forest.GetGid(d); + if (gid >= 0) { + int leaf_gid = forest.GetLeafGid(d); + pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], gid, + std::array{0, 0, 0}, 0, 0, 0, 0); + } + } + } + + SetMeshBlockNeighbors(GridIdentifier::two_level_composite(level), bl, ranklist); } } @@ -312,14 +327,14 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } } } -*/ + // Find same level neighbors on all GMG levels auto root_grid = this->GetRootGridInfo(); for (int gmg_level = gmg_min_level; gmg_level <= current_level; ++gmg_level) { SetSameLevelNeighbors(gmg_block_lists[gmg_level], gmg_grid_locs[gmg_level], root_grid, nbs, true, gmg_level); } - /* + // Now find GMG coarser neighbor for (int gmg_level = gmg_min_level + 1; gmg_level <= current_level; ++gmg_level) { for (auto &pmb : gmg_block_lists[gmg_level]) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 401e134e4efe..59df71c4787f 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -427,7 +427,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, packages, resolved_packages, gflag); } - SetMeshBlockNeighbors(block_list, nbs, ranklist); + SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); BuildGMGHierarchy(nbs, pin, app_in); ResetLoadBalanceVariables(); } @@ -685,7 +685,7 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, packages, resolved_packages, gflag, costlist[i]); } - SetMeshBlockNeighbors(block_list, nbs, ranklist); + SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); BuildGMGHierarchy(nbs, pin, app_in); ResetLoadBalanceVariables(); } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 11208407c120..860a7bc79220 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -320,7 +320,7 @@ class Mesh { int ntot); void BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app_in); void - SetMeshBlockNeighbors(BlockList_t &block_list, int nbs, + SetMeshBlockNeighbors(GridIdentifier grid_id, BlockList_t &block_list, const std::vector &ranklist, const std::unordered_set &newly_refined = {}); void From 37b5557508f3add76a1da2b646d20b9598b04f33 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 10:59:11 -0600 Subject: [PATCH 14/39] seemingly working --- src/mesh/forest/logical_location.hpp | 3 ++- src/mesh/forest/tree.cpp | 12 +++++++++--- src/mesh/mesh-gmg.cpp | 21 +++++++++++++-------- src/mesh/mesh.hpp | 2 +- 4 files changed, 25 insertions(+), 13 deletions(-) diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index a81a07a0d952..95871fc9ee3a 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -102,7 +102,7 @@ class LogicalLocation { // aggregate and POD type // possibly including a ghost block halo around the tree bool IsInTree(int nghost = 0) const { const int low = -nghost; - const int up = (1LL << level()) + nghost; + const int up = (1LL << std::max(level(), 0)) + nghost; return (l_[0] >= low) && (l_[0] < up) && (l_[1] >= low) && (l_[1] < up) && (l_[2] >= low) && (l_[2] < up); } @@ -111,6 +111,7 @@ class LogicalLocation { // aggregate and POD type bool IsInHalo(int nghost) const { return IsInTree(nghost) && !IsInTree(0); } int NeighborTreeIndex() const { + int up = 1LL << std::max(level(), 0); int i1 = (l_[0] >= 0) - (l_[0] < (1LL << level())) + 1; int i2 = (l_[1] >= 0) - (l_[1] < (1LL << level())) + 1; int i3 = (l_[2] >= 0) - (l_[2] < (1LL << level())) + 1; diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index 84d1c54a2365..284b54e2a2f8 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -54,6 +54,12 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz } } } + + // Build in negative levels + for (int l = -20; l < 0; ++l) { + internal_nodes.emplace(LocMapEntry(LogicalLocation(my_id, l, 0, 0, 0), -1, -1)); + } + // Pre-populate the next finest two-level composite grid with coarser blocks gmg_tlc_grids[root_level + 1] = gmg_tlc_grids[root_level]; } @@ -165,7 +171,7 @@ void Tree::FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int o PARTHENON_REQUIRE( loc.tree() == my_id, "Trying to find neighbors in a tree with a LogicalLocation on a different tree."); - PARTHENON_REQUIRE(leaves.count(loc) == 1, "Location must be a leaf to find neighbors."); + PARTHENON_REQUIRE((leaves.count(loc) == 1 || internal_nodes.count(loc) == 1), "Location must be in the tree to find neighbors."); auto neigh = loc.GetSameLevelNeighbor(ox1, ox2, ox3); int n_idx = neigh.NeighborTreeIndex(); @@ -307,9 +313,9 @@ RegionSize Tree::GetBlockDomain(const LogicalLocation &loc) const { // Negative logical levels correspond to reduced block sizes covering the entire // domain. auto reduction_fac = 1LL << (-loc.level()); - out.nx(dir) = domain.nx(dir) / reduction_fac; - PARTHENON_REQUIRE(out.nx(dir) % reduction_fac == 0, + PARTHENON_REQUIRE(domain.nx(dir) % reduction_fac == 0, "Trying to go to too large of a negative level."); + out.nx(dir) = domain.nx(dir) / reduction_fac; } } // If this is a translational symmetry direction, set the cell to cover the entire diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 561dee7f948c..54ea4f3a9848 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -212,16 +212,16 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (!mesh_size.symmetry(dir)) { int dir_allowed_levels = - NumberOfBinaryTrailingZeros(block_size_default.nx(dir) * nrbx[dir - 1]); + NumberOfBinaryTrailingZeros(block_size_default.nx(dir)); gmg_level_offset = std::min(dir_allowed_levels, gmg_level_offset); } } - const int gmg_min_level = root_level - gmg_level_offset; + const int gmg_min_level = -gmg_level_offset; gmg_min_logical_level_ = gmg_min_level; - + printf("gmg_min_level = %i gmg_max_level = %i\n", gmg_min_level, current_level); for (int level = gmg_min_level; level <= current_level; ++level) { - gmg_grid_locs[level] = LogicalLocMap_t(); + //gmg_grid_locs[level] = LogicalLocMap_t(); gmg_block_lists[level] = BlockList_t(); gmg_mesh_data[level] = DataCollection>(); } @@ -244,7 +244,8 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app // Create internal blocks that share a Morton number with this block // and add them to gmg two-level composite grid block lists auto loc = pmb->loc.GetParent(); - while (loc.level() <= gmg_min_level && loc.morton() == pmb->loc.morton()) { + while (loc.level() >= gmg_min_level && loc.morton() == pmb->loc.morton()) { + printf("adding level %i\n", loc.level()); RegionSize block_size = GetBlockSize(); BoundaryFlag block_bcs[6]; SetBlockSizeAndBoundaries(loc, block_size, block_bcs); @@ -259,7 +260,8 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app for (auto &[level, bl] : gmg_block_lists) { std::sort(bl.begin(), bl.end(), [](auto &a, auto &b){ return a->gid < b->gid;}); - for (auto &pmb : bl) { + for (auto &pmb : bl) { + // Coarser neighbor if (pmb->loc.level() > gmg_min_level) { auto ploc = pmb->loc.GetParent(); int gid = forest.GetGid(ploc); @@ -270,19 +272,22 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } } + // Finer neighbor(s) if (pmb->loc.level() < current_level) { auto dlocs = pmb->loc.GetDaughters(ndim); for (auto &d : dlocs) { int gid = forest.GetGid(d); if (gid >= 0) { int leaf_gid = forest.GetLeafGid(d); - pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], gid, - std::array{0, 0, 0}, 0, 0, 0, 0); + pmb->gmg_finer_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], gid, + std::array{0, 0, 0}, 0, 0, 0, 0); } } } + // Same level neighbors SetMeshBlockNeighbors(GridIdentifier::two_level_composite(level), bl, ranklist); + printf("gid: %i %s coarse: %i fine: %i neighbors: %i\n", pmb->gid, pmb->loc.label().c_str(), pmb->gmg_coarser_neighbors.size(), pmb->gmg_finer_neighbors.size(), pmb->gmg_same_neighbors.size()); } } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 860a7bc79220..bb6f8753adc6 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -120,7 +120,7 @@ class Mesh { DataCollection> mesh_data; LogicalLocMap_t leaf_grid_locs; - std::map gmg_grid_locs; + // std::map gmg_grid_locs; std::map gmg_block_lists; std::map>> gmg_mesh_data; int GetGMGMaxLevel() { return current_level; } From e1c99d80b9efe58f95f891b45f6db09af1fbee5a Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:09:27 -0600 Subject: [PATCH 15/39] remove old gmg stuff --- src/mesh/mesh-gmg.cpp | 201 ++---------------------------------------- 1 file changed, 6 insertions(+), 195 deletions(-) diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 54ea4f3a9848..6c6c2a460b3d 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -104,111 +104,13 @@ void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, } } -void Mesh::SetSameLevelNeighbors( - BlockList_t &block_list, const LogicalLocMap_t &loc_map, RootGridInfo root_grid, - int nbs, bool gmg_neighbors, int composite_logical_level, - const std::unordered_set &newly_refined) { - BufferID buffer_id(ndim, multilevel); - - for (auto &pmb : block_list) { - auto loc = pmb->loc; - auto gid = pmb->gid; - auto *neighbor_list = gmg_neighbors ? &(pmb->gmg_same_neighbors) : &(pmb->neighbors); - if (gmg_neighbors && loc.level() == composite_logical_level - 1) { - neighbor_list = &(pmb->gmg_composite_finer_neighbors); - } else if (gmg_neighbors && loc.level() == composite_logical_level) { - neighbor_list = &(pmb->gmg_same_neighbors); - } else if (gmg_neighbors) { - PARTHENON_FAIL("GMG grid was build incorrectly."); - } - - *neighbor_list = {}; - - auto possible_neighbors = loc.GetPossibleNeighbors(root_grid); - int buf_id = 0; - for (auto &pos_neighbor_location : possible_neighbors) { - if (gmg_neighbors && loc.level() == composite_logical_level - 1 && - loc.level() == pos_neighbor_location.level()) - continue; - if (loc_map.count(pos_neighbor_location) > 0) { - const auto &gid_rank = loc_map.at(pos_neighbor_location); - auto offsets = loc.GetSameLevelOffsets(pos_neighbor_location, root_grid); - // This inner loop is necessary in case a block pair has multiple neighbor - // connections due to periodic boundaries - for (auto ox1 : offsets[0]) { - for (auto ox2 : offsets[1]) { - for (auto ox3 : offsets[2]) { - NeighborConnect nc; - if (pos_neighbor_location.level() != loc.level()) { - // Check that the two blocks are in fact neighbors in this offset - // direction, since we only checked that they are actually neighbors - // when they have both been derefined to the coarser of their levels - // (this should only play a role in small meshes with periodic - // bounday conditions) - auto &fine_loc = pos_neighbor_location.level() > loc.level() - ? pos_neighbor_location - : loc; - int mult = loc.level() - pos_neighbor_location.level(); - std::array ox{ox1, ox2, ox3}; - bool not_neighbor = false; - for (int dir = 0; dir < 3; ++dir) { - if (ox[dir] != 0) { - // temp should be +1 if a block is to the right within its parent - // block and -1 if it is to the left. - const int temp = - 2 * (fine_loc.l(dir) - 2 * (fine_loc.l(dir) >> 1)) - 1; - PARTHENON_DEBUG_REQUIRE(temp * temp == 1, "Bad Offset"); - if (temp != mult * ox[dir]) not_neighbor = true; - } - } - if (not_neighbor) continue; - } - int connect_indicator = std::abs(ox1) + std::abs(ox2) + std::abs(ox3); - if (connect_indicator == 0) continue; - if (connect_indicator == 1) { - nc = NeighborConnect::face; - } else if (connect_indicator == 2) { - nc = NeighborConnect::edge; - } else if (connect_indicator == 3) { - nc = NeighborConnect::corner; - } - auto f = loc.GetAthenaXXFaceOffsets(pos_neighbor_location, ox1, ox2, ox3); - auto fn = - pos_neighbor_location.GetAthenaXXFaceOffsets(loc, -ox1, -ox2, -ox3); - int bid = buffer_id.GetID(ox1, ox2, ox3, f[0], f[1]); - int tid = buffer_id.GetID(-ox1, -ox2, -ox3, fn[0], fn[1]); - neighbor_list->emplace_back( - pmb->pmy_mesh, pos_neighbor_location, gid_rank.second, gid_rank.first, - gid_rank.first - nbs, std::array{ox1, ox2, ox3}, nc, bid, tid, - f[0], f[1]); - } - } - } - } - } - - std::sort(neighbor_list->begin(), neighbor_list->end(), - [](auto lhs, auto rhs) { return lhs.bufid < rhs.bufid; }); - - // Set neighbor block ownership - std::unordered_set allowed_neighbors; - allowed_neighbors.insert(pmb->loc); - for (auto &nb : *neighbor_list) - allowed_neighbors.insert(nb.loc); - for (auto &nb : *neighbor_list) { - nb.ownership = - DetermineOwnership(nb.loc, allowed_neighbors, root_grid, newly_refined); - nb.ownership.initialized = true; - } - } -} - void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app_in) { if (!multigrid) return; - // Create GMG logical location lists, first just copy coarsest grid - auto block_size_default = GetBlockSize(); - + + // See how many times we can go below logical level zero based on the + // number of times a blocks zones can be reduced by 2^D int gmg_level_offset = std::numeric_limits::max(); + auto block_size_default = GetBlockSize(); for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (!mesh_size.symmetry(dir)) { int dir_allowed_levels = @@ -219,9 +121,7 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app const int gmg_min_level = -gmg_level_offset; gmg_min_logical_level_ = gmg_min_level; - printf("gmg_min_level = %i gmg_max_level = %i\n", gmg_min_level, current_level); for (int level = gmg_min_level; level <= current_level; ++level) { - //gmg_grid_locs[level] = LogicalLocMap_t(); gmg_block_lists[level] = BlockList_t(); gmg_mesh_data[level] = DataCollection>(); } @@ -242,10 +142,10 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } // Create internal blocks that share a Morton number with this block - // and add them to gmg two-level composite grid block lists + // and add them to gmg two-level composite grid block lists. This + // determines which process internal blocks live on auto loc = pmb->loc.GetParent(); while (loc.level() >= gmg_min_level && loc.morton() == pmb->loc.morton()) { - printf("adding level %i\n", loc.level()); RegionSize block_size = GetBlockSize(); BoundaryFlag block_bcs[6]; SetBlockSizeAndBoundaries(loc, block_size, block_bcs); @@ -287,96 +187,7 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app // Same level neighbors SetMeshBlockNeighbors(GridIdentifier::two_level_composite(level), bl, ranklist); - printf("gid: %i %s coarse: %i fine: %i neighbors: %i\n", pmb->gid, pmb->loc.label().c_str(), pmb->gmg_coarser_neighbors.size(), pmb->gmg_finer_neighbors.size(), pmb->gmg_same_neighbors.size()); - } - } - -/* - // Add leaf grid locations to GMG grid levels - int gmg_gid = 0; - for (auto loc : loclist) { - const int gmg_level = loc.level(); - gmg_grid_locs[gmg_level].insert( - {loc, std::pair(gmg_gid, ranklist[gmg_gid])}); - if (gmg_level < current_level) { - gmg_grid_locs[gmg_level + 1].insert( - {loc, std::pair(gmg_gid, ranklist[gmg_gid])}); - } - if (ranklist[gmg_gid] == Globals::my_rank) { - const int lid = gmg_gid - nslist[Globals::my_rank]; - gmg_block_lists[gmg_level].push_back(block_list[lid]); - if (gmg_level < current_level) - gmg_block_lists[gmg_level + 1].push_back(block_list[lid]); - } - gmg_gid++; - } - - // Fill in internal nodes for GMG grid levels from levels on finer GMG grid - for (int gmg_level = current_level - 1; gmg_level >= gmg_min_level; --gmg_level) { - for (auto &[loc, gid_rank] : gmg_grid_locs[gmg_level + 1]) { - if (loc.level() == gmg_level + 1) { - auto parent = loc.GetParent(); - if (parent.morton() == loc.morton()) { - gmg_grid_locs[gmg_level].insert( - {parent, std::make_pair(gmg_gid, gid_rank.second)}); - if (gid_rank.second == Globals::my_rank) { - BoundaryFlag block_bcs[6]; - auto block_size = block_size_default; - SetBlockSizeAndBoundaries(parent, block_size, block_bcs); - gmg_block_lists[gmg_level].push_back( - MeshBlock::Make(gmg_gid, -1, parent, block_size, block_bcs, this, pin, - app_in, packages, resolved_packages, gflag)); - } - gmg_gid++; - } - } - } - } - - // Find same level neighbors on all GMG levels - auto root_grid = this->GetRootGridInfo(); - for (int gmg_level = gmg_min_level; gmg_level <= current_level; ++gmg_level) { - SetSameLevelNeighbors(gmg_block_lists[gmg_level], gmg_grid_locs[gmg_level], root_grid, - nbs, true, gmg_level); - } - - // Now find GMG coarser neighbor - for (int gmg_level = gmg_min_level + 1; gmg_level <= current_level; ++gmg_level) { - for (auto &pmb : gmg_block_lists[gmg_level]) { - if (pmb->loc.level() != gmg_level) continue; - auto parent_loc = pmb->loc.GetParent(); - auto loc = pmb->loc; - auto gid = pmb->gid; - auto rank = Globals::my_rank; - if (gmg_grid_locs[gmg_level - 1].count(parent_loc) > 0) { - loc = parent_loc; - gid = gmg_grid_locs[gmg_level - 1][parent_loc].first; - rank = gmg_grid_locs[gmg_level - 1][parent_loc].second; - } else { - PARTHENON_FAIL("There is something wrong with GMG block list."); - } - pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, loc, rank, gid, gid - nbs, - std::array{0, 0, 0}, - NeighborConnect::none, 0, 0, 0, 0); - } - } - - // Now find finer GMG neighbors - for (int gmg_level = gmg_min_level; gmg_level <= current_level - 1; ++gmg_level) { - for (auto &pmb : gmg_block_lists[gmg_level]) { - if (pmb->loc.level() != gmg_level) continue; - auto daughter_locs = pmb->loc.GetDaughters(); - for (auto &daughter_loc : daughter_locs) { - if (gmg_grid_locs[gmg_level + 1].count(daughter_loc) > 0) { - auto &gid_rank = gmg_grid_locs[gmg_level + 1][daughter_loc]; - pmb->gmg_finer_neighbors.emplace_back( - pmb->pmy_mesh, daughter_loc, gid_rank.second, gid_rank.first, - gid_rank.first - nbs, std::array{0, 0, 0}, NeighborConnect::none, 0, - 0, 0, 0); - } - } } } - */ } } // namespace parthenon From 81e5a1858d74fd05ad76bdd1252ecfc407a20f0f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:23:04 -0600 Subject: [PATCH 16/39] remove a bunch of uneccesary stuff --- src/mesh/forest/block_ownership.cpp | 42 ------ src/mesh/forest/block_ownership.hpp | 6 - src/mesh/forest/logical_location.cpp | 197 --------------------------- src/mesh/forest/logical_location.hpp | 31 +---- src/mesh/mesh-amr_loadbalance.cpp | 1 - src/mesh/mesh-gmg.cpp | 8 -- src/mesh/mesh.cpp | 2 - src/mesh/mesh.hpp | 4 - tst/unit/test_logical_location.cpp | 68 --------- 9 files changed, 1 insertion(+), 358 deletions(-) diff --git a/src/mesh/forest/block_ownership.cpp b/src/mesh/forest/block_ownership.cpp index 967ef5829659..047496593b9b 100644 --- a/src/mesh/forest/block_ownership.cpp +++ b/src/mesh/forest/block_ownership.cpp @@ -37,48 +37,6 @@ namespace parthenon { -// TODO(LFR): Remove this -block_ownership_t -DetermineOwnership(const LogicalLocation &main_block, - const std::unordered_set &allowed_neighbors, - const RootGridInfo &rg_info, - const std::unordered_set &newly_refined) { - block_ownership_t main_owns; - - auto ownership_level = [&](const LogicalLocation &a) { - // Newly-refined blocks are treated as higher-level than blocks at their - // parent level, but lower-level than previously-refined blocks at their - // current level. - if (newly_refined.count(a)) return 2 * a.level() - 1; - return 2 * a.level(); - }; - - auto ownership_less_than = [ownership_level](const LogicalLocation &a, - const LogicalLocation &b) { - // Ownership is first determined by block with the highest level, then by maximum - // Morton number this is reversed in precedence from the normal comparators where - // Morton number takes precedence - if (ownership_level(a) == ownership_level(b)) return a.morton() < b.morton(); - return ownership_level(a) < ownership_level(b); - }; - - for (int ox1 : {-1, 0, 1}) { - for (int ox2 : {-1, 0, 1}) { - for (int ox3 : {-1, 0, 1}) { - main_owns(ox1, ox2, ox3) = true; - for (auto &n : allowed_neighbors) { - if (ownership_less_than(main_block, n) && - main_block.IsNeighborOfTE(n, ox1, ox2, ox3, rg_info)) { - main_owns(ox1, ox2, ox3) = false; - break; - } - } - } - } - } - return main_owns; -} - block_ownership_t DetermineOwnershipForest(const LogicalLocation &main_block, const std::vector &allowed_neighbors, diff --git a/src/mesh/forest/block_ownership.hpp b/src/mesh/forest/block_ownership.hpp index 90c20ff1d057..d03fa86c6d41 100644 --- a/src/mesh/forest/block_ownership.hpp +++ b/src/mesh/forest/block_ownership.hpp @@ -82,12 +82,6 @@ struct block_ownership_t { bool ownership[3][3][3]; }; -block_ownership_t -DetermineOwnership(const LogicalLocation &main_block, - const std::unordered_set &allowed_neighbors, - const RootGridInfo &rg_info = RootGridInfo(), - const std::unordered_set &newly_refined = {}); - block_ownership_t DetermineOwnershipForest(const LogicalLocation &main_block, const std::vector &allowed_neighbors, diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index 7c85d801a492..c33936fd99d3 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -65,36 +65,6 @@ LogicalLocation::GetSameLevelOffsetsForest(const LogicalLocation &neighbor) cons return offsets; } -std::array, 3> -LogicalLocation::GetSameLevelOffsets(const LogicalLocation &neighbor, - const RootGridInfo &rg_info) const { - std::array, 3> offsets; - const int level_diff_1 = std::max(neighbor.level() - level(), 0); - const int level_diff_2 = std::max(level() - neighbor.level(), 0); - const int n_per_root_block = - 1 << std::max((std::min(level(), neighbor.level()) - rg_info.level), 0); - const int root_block_per_n = - 1 << std::max(rg_info.level - std::min(level(), neighbor.level()), 0); - for (int i = 0; i < 3; ++i) { - const auto idxt = l(i) >> level_diff_2; - const auto idxn = neighbor.l(i) >> level_diff_1; - if (std::abs(idxn - idxt) <= 1) offsets[i].push_back(idxn - idxt); - - int n_blocks_level = std::max(n_per_root_block * rg_info.n[i], 1); - if (root_block_per_n > 1) - n_blocks_level = - rg_info.n[i] / root_block_per_n + (rg_info.n[i] % root_block_per_n != 0); - if (rg_info.periodic[i]) { - if (std::abs(idxn - n_blocks_level - idxt) <= 1) - offsets[i].push_back(idxn - n_blocks_level - idxt); - if (std::abs(idxn + n_blocks_level - idxt) <= 1) - offsets[i].push_back(idxn + n_blocks_level - idxt); - } - } - - return offsets; -} - bool LogicalLocation::IsNeighborForest(const LogicalLocation &in) const { PARTHENON_REQUIRE(tree() == in.tree(), "Trying to compare locations not in the same octree."); @@ -145,72 +115,6 @@ bool LogicalLocation::IsNeighborOfTEForest(const LogicalLocation &in, return neighbors; } -template -bool LogicalLocation::NeighborFindingImpl(const LogicalLocation &in, - const std::array &te_offset, - const RootGridInfo &rg_info) const { - if (in.level() >= level() && Contains(in)) return false; // You share a volume - if (in.level() < level() && in.Contains(*this)) return false; // You share a volume - - // We work on the finer level of in.level() and this->level() - const int max_level = std::max(in.level(), level()); - const int level_shift_1 = max_level - level(); - const int level_shift_2 = max_level - in.level(); - const auto block_size_1 = 1 << level_shift_1; - const auto block_size_2 = 1 << level_shift_2; - - // TODO(LFR): Think about what this should do when we are above the root level - const int n_per_root_block = 1 << std::max(max_level - rg_info.level, 0); - const int root_block_per_n = 1 << std::max(rg_info.level - max_level, 0); - std::array b; - - for (int i = 0; i < 3; ++i) { - // Index range of daughters of this block on current level plus a one block halo on - // either side - auto low = (l(i) << level_shift_1) - 1; - auto hi = low + block_size_1 + 1; - // Indexing for topological offset calculation - if constexpr (TENeighbor) { - if (te_offset[i] == -1) { - // Left side offset, so only two possible block indices are allowed in this - // direction - hi -= block_size_1; - } else if (te_offset[i] == 0) { - // No offset in this direction, so only interior - low += 1; - hi -= 1; - } else { - // Right side offset, so only two possible block indices are allowed in this - // direction - low += block_size_1; - } - } - // Index range of daughters of possible neighbor block on current level - const auto in_low = in.l(i) << level_shift_2; - const auto in_hi = in_low + block_size_2 - 1; - // Check if these two ranges overlap at all - b[i] = in_hi >= low && in_low <= hi; - if (rg_info.periodic[i]) { - int n_cells_level = std::max(n_per_root_block * rg_info.n[i], 1); - if (root_block_per_n > 1) - n_cells_level = - rg_info.n[i] / root_block_per_n + (rg_info.n[i] % root_block_per_n != 0); - b[i] = b[i] || (in_hi + n_cells_level >= low && in_low + n_cells_level <= hi); - b[i] = b[i] || (in_hi - n_cells_level >= low && in_low - n_cells_level <= hi); - } - } - - return b[0] && b[1] && b[2]; -} -template bool -LogicalLocation::NeighborFindingImpl(const LogicalLocation &in, - const std::array &te_offset, - const RootGridInfo &rg_info) const; -template bool -LogicalLocation::NeighborFindingImpl(const LogicalLocation &in, - const std::array &te_offset, - const RootGridInfo &rg_info) const; - std::vector LogicalLocation::GetDaughters(int ndim) const { std::vector daughters; if (level() < 0) { @@ -231,106 +135,5 @@ std::vector LogicalLocation::GetDaughters(int ndim) const { return daughters; } -std::unordered_set -LogicalLocation::GetPossibleNeighbors(const RootGridInfo &rg_info) { - const std::vector irange{-1, 0, 1}; - const std::vector jrange{-1, 0, 1}; - const std::vector krange{-1, 0, 1}; - const std::vector daughter_irange{0, 1}; - const std::vector daughter_jrange{0, 1}; - const std::vector daughter_krange{0, 1}; - - return GetPossibleNeighborsImpl(irange, jrange, krange, daughter_irange, - daughter_jrange, daughter_krange, rg_info); -} - -std::unordered_set -LogicalLocation::GetPossibleBlocksSurroundingTopologicalElement( - int ox1, int ox2, int ox3, const RootGridInfo &rg_info) const { - const auto irange = - (std::abs(ox1) == 1) ? std::vector{0, ox1} : std::vector{0}; - const auto jrange = - (std::abs(ox2) == 1) ? std::vector{0, ox2} : std::vector{0}; - const auto krange = - (std::abs(ox3) == 1) ? std::vector{0, ox3} : std::vector{0}; - const auto daughter_irange = - (std::abs(ox1) == 1) ? std::vector{ox1 > 0} : std::vector{0, 1}; - const auto daughter_jrange = - (std::abs(ox2) == 1) ? std::vector{ox2 > 0} : std::vector{0, 1}; - const auto daughter_krange = - (std::abs(ox3) == 1) ? std::vector{ox3 > 0} : std::vector{0, 1}; - - return GetPossibleNeighborsImpl(irange, jrange, krange, daughter_irange, - daughter_jrange, daughter_krange, rg_info); -} - -std::unordered_set LogicalLocation::GetPossibleNeighborsImpl( - const std::vector &irange, const std::vector &jrange, - const std::vector &krange, const std::vector &daughter_irange, - const std::vector &daughter_jrange, const std::vector &daughter_krange, - const RootGridInfo &rg_info) const { - std::vector locs; - - auto AddNeighbors = [&](const LogicalLocation &loc, bool include_parents) { - const int n_per_root_block = 1 << std::max(loc.level() - rg_info.level, 0); - const int down_shift = 1 << std::max(rg_info.level - loc.level(), 0); - // Account for the fact that the root grid may be overhanging into a partial block - const int extra1 = (rg_info.n[0] % down_shift > 0); - const int extra2 = (rg_info.n[1] % down_shift > 0); - const int extra3 = (rg_info.n[2] % down_shift > 0); - int n1_cells_level = - std::max(n_per_root_block * (rg_info.n[0] / down_shift + extra1), 1); - int n2_cells_level = - std::max(n_per_root_block * (rg_info.n[1] / down_shift + extra2), 1); - int n3_cells_level = - std::max(n_per_root_block * (rg_info.n[2] / down_shift + extra3), 1); - for (int i : irange) { - for (int j : jrange) { - for (int k : krange) { - auto lx1 = loc.lx1() + i; - auto lx2 = loc.lx2() + j; - auto lx3 = loc.lx3() + k; - // This should include blocks that are connected by periodic boundaries - if (rg_info.periodic[0]) lx1 = (lx1 + n1_cells_level) % n1_cells_level; - if (rg_info.periodic[1]) lx2 = (lx2 + n2_cells_level) % n2_cells_level; - if (rg_info.periodic[2]) lx3 = (lx3 + n3_cells_level) % n3_cells_level; - if (0 <= lx1 && lx1 < n1_cells_level && 0 <= lx2 && lx2 < n2_cells_level && - 0 <= lx3 && lx3 < n3_cells_level) { - if (loc.level() > level()) { - const int s = loc.level() - level(); - if ((lx1 >> s) != this->lx1() || (lx2 >> s) != this->lx2() || - (lx3 >> s) != this->lx3()) { - locs.emplace_back(loc.level(), lx1, lx2, lx3); - } - } else { - locs.emplace_back(loc.level(), lx1, lx2, lx3); - } - if (include_parents) { - auto parent = locs.back().GetParent(); - if (IsNeighbor(parent, rg_info)) locs.push_back(parent); - } - } - } - } - } - }; - - // Find the same level and lower level blocks of this block - AddNeighbors(*this, true); - - // Iterate over daughters of this block that share the same topological element - for (int l : daughter_irange) { - for (int m : daughter_jrange) { - for (int n : daughter_krange) { - AddNeighbors(GetDaughter(l, m, n), false); - } - } - } - // The above procedure likely duplicated some blocks, so put them in a set - std::unordered_set unique_locs; - for (auto &loc : locs) - unique_locs.emplace(std::move(loc)); - return unique_locs; -} } // namespace parthenon diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 95871fc9ee3a..65fde3d5de02 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -132,26 +132,14 @@ class LogicalLocation { // aggregate and POD type // TODO(LFR): Remove the corresponding non-forest routine once GMG is working std::array GetSameLevelOffsetsForest(const LogicalLocation &neighbor) const; - std::array, 3> GetSameLevelOffsets(const LogicalLocation &neighbor, - const RootGridInfo &rg_info) const; + // Being a neighbor implies that you share a face, edge, or node and don't share a // volume - bool IsNeighbor(const LogicalLocation &in, - const RootGridInfo &rg_info = RootGridInfo()) const { - return NeighborFindingImpl(in, std::array(), rg_info); - } - - // TODO(LFR): Remove the corresponding non-forest routine once GMG is working bool IsNeighborForest(const LogicalLocation &in) const; // TODO(LFR): Remove the corresponding non-forest routine once GMG is working bool IsNeighborOfTEForest(const LogicalLocation &in, const std::array &te_offset) const; - bool IsNeighborOfTE(const LogicalLocation &in, int ox1, int ox2, int ox3, - const RootGridInfo &rg_info = RootGridInfo()) const { - return NeighborFindingImpl(in, std::array{ox1, ox2, ox3}, rg_info); - } - LogicalLocation GetSameLevelNeighbor(int ox1, int ox2, int ox3) const { return LogicalLocation(tree(), level(), lx1() + ox1, lx2() + ox2, lx3() + ox3); } @@ -187,23 +175,6 @@ class LogicalLocation { // aggregate and POD type } return f; } - - std::unordered_set - GetPossibleNeighbors(const RootGridInfo &rg_info = RootGridInfo()); - - std::unordered_set GetPossibleBlocksSurroundingTopologicalElement( - int ox1, int ox2, int ox3, const RootGridInfo &rg_info = RootGridInfo()) const; - - private: - template - bool NeighborFindingImpl(const LogicalLocation &in, const std::array &te_offset, - const RootGridInfo &rg_info = RootGridInfo()) const; - - std::unordered_set GetPossibleNeighborsImpl( - const std::vector &irange, const std::vector &jrange, - const std::vector &krange, const std::vector &daughter_irange, - const std::vector &daughter_jrange, const std::vector &daughter_krange, - const RootGridInfo &rg_info = RootGridInfo()) const; }; inline bool operator<(const LogicalLocation &lhs, const LogicalLocation &rhs) { diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index 32d98a0d9231..9e5776977927 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -927,7 +927,6 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput loclist = std::move(newloc); ranklist = std::move(newrank); costlist = std::move(newcost); - PopulateLeafLocationMap(); // Make sure all old sends/receives are done before we reconfigure the mesh #ifdef MPI_PARALLEL diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 6c6c2a460b3d..815d19a7e704 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -43,14 +43,6 @@ namespace parthenon { -void Mesh::PopulateLeafLocationMap() { - const int nbtot = ranklist.size(); - leaf_grid_locs.clear(); - for (int ib = 0; ib < nbtot; ++ib) { - leaf_grid_locs[loclist[ib]] = std::make_pair(ib, ranklist[ib]); - } -} - void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, BlockList_t &block_list, const std::vector &ranklist, const std::unordered_set &newly_refined) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 59df71c4787f..b6ea2ba53e81 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -394,7 +394,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, costlist = std::vector(nbtotal, 1.0); CalculateLoadBalance(costlist, ranklist, nslist, nblist); - PopulateLeafLocationMap(); // Output some diagnostic information to terminal @@ -648,7 +647,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, } CalculateLoadBalance(costlist, ranklist, nslist, nblist); - PopulateLeafLocationMap(); // Output MeshBlock list and quit (mesh test only); do not create meshes if (mesh_test > 0) { diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index bb6f8753adc6..8e5225c16479 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -119,13 +119,10 @@ class Mesh { DataCollection> mesh_data; - LogicalLocMap_t leaf_grid_locs; - // std::map gmg_grid_locs; std::map gmg_block_lists; std::map>> gmg_mesh_data; int GetGMGMaxLevel() { return current_level; } int GetGMGMinLevel() { return gmg_min_logical_level_; } - int GetGMGMinLogicalLevel() { return gmg_min_logical_level_; } // functions void Initialize(bool init_problem, ParameterInput *pin, ApplicationInput *app_in); @@ -339,7 +336,6 @@ class Mesh { void RegisterLoadBalancing_(ParameterInput *pin); void SetupMPIComms(); - void PopulateLeafLocationMap(); void BuildTagMapAndBoundaryBuffers(); void CommunicateBoundaries(std::string md_name = "base"); void PreCommFillDerived(); diff --git a/tst/unit/test_logical_location.cpp b/tst/unit/test_logical_location.cpp index 4b17f73bb5fd..045719898fc1 100644 --- a/tst/unit/test_logical_location.cpp +++ b/tst/unit/test_logical_location.cpp @@ -186,74 +186,6 @@ TEST_CASE("Logical Location", "[Logical Location]") { } } - THEN("We can correctly find the blocks adjacent to a face") { - auto base_loc = LogicalLocation(2, 2, 3, 3); - - auto possible_neighbors = - base_loc.GetPossibleBlocksSurroundingTopologicalElement(1, 0, 0); - std::unordered_set by_hand_elements, automatic_elements; - // There should be five total neighboring blocks of this face since one neighbor is - // refined - by_hand_elements.insert(LogicalLocation(2, 2, 3, 3)); - by_hand_elements.insert(LogicalLocation(3, 6, 6, 6)); - by_hand_elements.insert(LogicalLocation(3, 6, 6, 7)); - by_hand_elements.insert(LogicalLocation(3, 6, 7, 6)); - by_hand_elements.insert(LogicalLocation(3, 6, 7, 7)); - - for (auto &n : possible_neighbors) { - if (hash_leaves.count(n) > 0) { - automatic_elements.insert(n); - } - } - REQUIRE(by_hand_elements == automatic_elements); - } - - THEN("We can correctly find the blocks adjacent to an edge") { - auto base_loc = LogicalLocation(2, 2, 3, 3); - - auto possible_neighbors = - base_loc.GetPossibleBlocksSurroundingTopologicalElement(1, -1, 0); - std::unordered_set by_hand_elements, automatic_elements; - // There should be five total neighboring blocks of this edge since one neighbor is - // refined - by_hand_elements.insert(LogicalLocation(2, 2, 2, 3)); - by_hand_elements.insert(LogicalLocation(2, 3, 2, 3)); - by_hand_elements.insert(LogicalLocation(2, 2, 3, 3)); - by_hand_elements.insert(LogicalLocation(3, 6, 6, 6)); - by_hand_elements.insert(LogicalLocation(3, 6, 6, 7)); - - for (auto &n : possible_neighbors) { - if (hash_leaves.count(n) > 0) { - automatic_elements.insert(n); - } - } - REQUIRE(by_hand_elements == automatic_elements); - } - - THEN("We can correctly find the blocks adjacent to a node") { - auto base_loc = LogicalLocation(2, 2, 3, 3); - - auto possible_neighbors = - base_loc.GetPossibleBlocksSurroundingTopologicalElement(1, -1, -1); - std::unordered_set by_hand_elements, automatic_elements; - // There should be eight total neighboring blocks for this node - by_hand_elements.insert(LogicalLocation(2, 2, 2, 3)); - by_hand_elements.insert(LogicalLocation(2, 2, 2, 2)); - by_hand_elements.insert(LogicalLocation(2, 3, 2, 3)); - by_hand_elements.insert(LogicalLocation(2, 3, 2, 2)); - by_hand_elements.insert(LogicalLocation(2, 2, 3, 3)); - by_hand_elements.insert(LogicalLocation(2, 2, 3, 2)); - by_hand_elements.insert(LogicalLocation(3, 6, 6, 6)); - by_hand_elements.insert(LogicalLocation(2, 3, 3, 2)); - - for (auto &n : possible_neighbors) { - if (hash_leaves.count(n) > 0) { - automatic_elements.insert(n); - } - } - REQUIRE(by_hand_elements == automatic_elements); - } - THEN("We can find the ownership array of a block") { LogicalLocation base_loc(2, 2, 3, 3); auto owns = DetermineOwnershipForest(base_loc, neighbor_locs); From 933f48dc282cc06195c7ae42bf418594a0f30d3d Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:28:46 -0600 Subject: [PATCH 17/39] Remove dead code --- src/mesh/mesh.cpp | 63 ++++------------------------------------------- 1 file changed, 5 insertions(+), 58 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index b6ea2ba53e81..8bc1d38f6213 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1183,31 +1183,10 @@ std::shared_ptr Mesh::FindMeshBlock(int tgid) const { bool Mesh::SetBlockSizeAndBoundaries(LogicalLocation loc, RegionSize &block_size, BoundaryFlag *block_bcs) { bool valid_region = true; - block_size = GetBlockSize(loc); - if (loc.tree() >= 0) { - auto bcs = forest.GetBlockBCs(loc); - for (int i = 0; i < BOUNDARY_NFACES; ++i) - block_bcs[i] = bcs[i]; - return valid_region; - } - - for (auto &dir : {X1DIR, X2DIR, X3DIR}) { - if (!block_size.symmetry(dir)) { - std::int64_t nrbx_ll = nrbx[dir - 1] << (loc.level() - root_level); - if (loc.level() < root_level) { - std::int64_t fac = 1 << (root_level - loc.level()); - nrbx_ll = nrbx[dir - 1] / fac + (nrbx[dir - 1] % fac != 0); - } - block_bcs[GetInnerBoundaryFace(dir)] = - loc.l(dir - 1) == 0 ? mesh_bcs[GetInnerBoundaryFace(dir)] : BoundaryFlag::block; - block_bcs[GetOuterBoundaryFace(dir)] = loc.l(dir - 1) == nrbx_ll - 1 - ? mesh_bcs[GetOuterBoundaryFace(dir)] - : BoundaryFlag::block; - } else { - block_bcs[GetInnerBoundaryFace(dir)] = mesh_bcs[GetInnerBoundaryFace(dir)]; - block_bcs[GetOuterBoundaryFace(dir)] = mesh_bcs[GetOuterBoundaryFace(dir)]; - } - } + block_size = forest.GetBlockDomain(loc); + auto bcs = forest.GetBlockBCs(loc); + for (int i = 0; i < BOUNDARY_NFACES; ++i) + block_bcs[i] = bcs[i]; return valid_region; } @@ -1218,39 +1197,7 @@ bool Mesh::SetBlockSizeAndBoundaries(LogicalLocation loc, RegionSize &block_size RegionSize Mesh::GetBlockSize(const LogicalLocation &loc) const { // TODO(LFR): Update this - if (loc.tree() >= 0) { - // Implies this is a location in a forest, not in the old Athena tree - return forest.GetBlockDomain(loc); - } - RegionSize block_size = GetBlockSize(); - for (auto &dir : {X1DIR, X2DIR, X3DIR}) { - block_size.xrat(dir) = mesh_size.xrat(dir); - block_size.symmetry(dir) = mesh_size.symmetry(dir); - if (!block_size.symmetry(dir)) { - std::int64_t nrbx_ll = nrbx[dir - 1] << (loc.level() - root_level); - if (loc.level() < root_level) { - std::int64_t fac = 1 << (root_level - loc.level()); - nrbx_ll = nrbx[dir - 1] / fac + (nrbx[dir - 1] % fac != 0); - } - block_size.xmin(dir) = GetMeshCoordinate(dir, BlockLocation::Left, loc); - block_size.xmax(dir) = GetMeshCoordinate(dir, BlockLocation::Right, loc); - // Correct for possible overshooting, since the root grid may not cover the - // entire logical level zero block of the mesh - if (block_size.xmax(dir) > mesh_size.xmax(dir) || loc.level() < 0) { - // Need integer reduction factor, so transform location back to root level - PARTHENON_REQUIRE(loc.level() < root_level, "Something is messed up."); - std::int64_t loc_low = loc.l(dir - 1) << (root_level - loc.level()); - std::int64_t loc_hi = (loc.l(dir - 1) + 1) << (root_level - loc.level()); - block_size.nx(dir) = - block_size.nx(dir) * (nrbx[dir - 1] - loc_low) / (loc_hi - loc_low); - block_size.xmax(dir) = mesh_size.xmax(dir); - } - } else { - block_size.xmin(dir) = mesh_size.xmin(dir); - block_size.xmax(dir) = mesh_size.xmax(dir); - } - } - return block_size; + return forest.GetBlockDomain(loc); } std::int64_t Mesh::GetTotalCells() { From 8966e48ccb46d4a070ee2c314d9c05703757b00e Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:32:44 -0600 Subject: [PATCH 18/39] Remove forest naming --- src/mesh/forest/block_ownership.cpp | 4 ++-- src/mesh/forest/block_ownership.hpp | 2 +- src/mesh/forest/logical_location.cpp | 6 +++--- src/mesh/forest/logical_location.hpp | 8 +++----- src/mesh/forest/tree.cpp | 4 ++-- src/mesh/mesh-gmg.cpp | 4 ++-- tst/unit/test_logical_location.cpp | 8 ++++---- 7 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/mesh/forest/block_ownership.cpp b/src/mesh/forest/block_ownership.cpp index 047496593b9b..f61093bfaaf2 100644 --- a/src/mesh/forest/block_ownership.cpp +++ b/src/mesh/forest/block_ownership.cpp @@ -38,7 +38,7 @@ namespace parthenon { block_ownership_t -DetermineOwnershipForest(const LogicalLocation &main_block, +DetermineOwnership(const LogicalLocation &main_block, const std::vector &allowed_neighbors, const std::unordered_set &newly_refined) { block_ownership_t main_owns; @@ -68,7 +68,7 @@ DetermineOwnershipForest(const LogicalLocation &main_block, main_owns(ox1, ox2, ox3) = true; for (const auto &n : allowed_neighbors) { if (ownership_less_than(main_block, n.global_loc) && - main_block.IsNeighborOfTEForest(n.origin_loc, {ox1, ox2, ox3})) { + main_block.IsNeighborOfTE(n.origin_loc, {ox1, ox2, ox3})) { main_owns(ox1, ox2, ox3) = false; break; } diff --git a/src/mesh/forest/block_ownership.hpp b/src/mesh/forest/block_ownership.hpp index d03fa86c6d41..4ed9b0308c10 100644 --- a/src/mesh/forest/block_ownership.hpp +++ b/src/mesh/forest/block_ownership.hpp @@ -83,7 +83,7 @@ struct block_ownership_t { }; block_ownership_t -DetermineOwnershipForest(const LogicalLocation &main_block, +DetermineOwnership(const LogicalLocation &main_block, const std::vector &allowed_neighbors, const std::unordered_set &newly_refined = {}); diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index c33936fd99d3..7119f26321b0 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -54,7 +54,7 @@ bool LogicalLocation::Contains(const LogicalLocation &containee) const { } std::array -LogicalLocation::GetSameLevelOffsetsForest(const LogicalLocation &neighbor) const { +LogicalLocation::GetSameLevelOffsets(const LogicalLocation &neighbor) const { std::array offsets; const int level_shift_neigh = std::max(neighbor.level() - level(), 0); const int level_shift_me = std::max(level() - neighbor.level(), 0); @@ -65,7 +65,7 @@ LogicalLocation::GetSameLevelOffsetsForest(const LogicalLocation &neighbor) cons return offsets; } -bool LogicalLocation::IsNeighborForest(const LogicalLocation &in) const { +bool LogicalLocation::IsNeighbor(const LogicalLocation &in) const { PARTHENON_REQUIRE(tree() == in.tree(), "Trying to compare locations not in the same octree."); const int max_level = std::max(in.level(), level()); @@ -86,7 +86,7 @@ bool LogicalLocation::IsNeighborForest(const LogicalLocation &in) const { return neighbors; } -bool LogicalLocation::IsNeighborOfTEForest(const LogicalLocation &in, +bool LogicalLocation::IsNeighborOfTE(const LogicalLocation &in, const std::array &te_offset) const { PARTHENON_REQUIRE(tree() == in.tree(), "Trying to compare locations not in the same octree."); diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 65fde3d5de02..083690a3a858 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -130,14 +130,12 @@ class LogicalLocation { // aggregate and POD type bool Contains(const LogicalLocation &containee) const; - // TODO(LFR): Remove the corresponding non-forest routine once GMG is working - std::array GetSameLevelOffsetsForest(const LogicalLocation &neighbor) const; + std::array GetSameLevelOffsets(const LogicalLocation &neighbor) const; // Being a neighbor implies that you share a face, edge, or node and don't share a // volume - bool IsNeighborForest(const LogicalLocation &in) const; - // TODO(LFR): Remove the corresponding non-forest routine once GMG is working - bool IsNeighborOfTEForest(const LogicalLocation &in, + bool IsNeighbor(const LogicalLocation &in) const; + bool IsNeighborOfTE(const LogicalLocation &in, const std::array &te_offset) const; LogicalLocation GetSameLevelNeighbor(int ox1, int ox2, int ox3) const { diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index 284b54e2a2f8..ae3eb0ac867e 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -208,7 +208,7 @@ void Tree::FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int o if (include_fine) { auto daughters = tneigh.GetDaughters(neighbor_tree->ndim); for (auto &n : daughters) { - if (tloc.IsNeighborForest(n)) + if (tloc.IsNeighbor(n)) neighbor_locs->push_back({n, orientation.TransformBack(n, GetId())}); } } else if (include_internal) { @@ -220,7 +220,7 @@ void Tree::FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int o // because our communication algorithm packs this extra data by hand, we do not wish // to duplicate coarser blocks in the neighbor list. Therefore, we only include the // coarse block in one offset position - auto sl_offset = loc.GetSameLevelOffsetsForest(neighp); + auto sl_offset = loc.GetSameLevelOffsets(neighp); if (sl_offset[0] == ox1 && sl_offset[1] == ox2 && sl_offset[2] == ox3) neighbor_locs->push_back({tneigh.GetParent(), neighp}); } diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 815d19a7e704..01c581a8dbdc 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -60,7 +60,7 @@ void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, int buf_id = 0; for (const auto &nloc : neighbors) { auto gid = forest.GetGid(nloc.global_loc); - auto offsets = loc.GetSameLevelOffsetsForest(nloc.origin_loc); + auto offsets = loc.GetSameLevelOffsets(nloc.origin_loc); auto f = loc.GetAthenaXXFaceOffsets(nloc.origin_loc, offsets[0], offsets[1], offsets[2]); int bid = buffer_id.GetID(offsets[0], offsets[1], offsets[2], f[0], f[1]); @@ -80,7 +80,7 @@ void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc); nb.ownership = - DetermineOwnershipForest(nloc.global_loc, neighbor_neighbors, newly_refined); + DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined); nb.ownership.initialized = true; } diff --git a/tst/unit/test_logical_location.cpp b/tst/unit/test_logical_location.cpp index 045719898fc1..665496d9012c 100644 --- a/tst/unit/test_logical_location.cpp +++ b/tst/unit/test_logical_location.cpp @@ -188,7 +188,7 @@ TEST_CASE("Logical Location", "[Logical Location]") { THEN("We can find the ownership array of a block") { LogicalLocation base_loc(2, 2, 3, 3); - auto owns = DetermineOwnershipForest(base_loc, neighbor_locs); + auto owns = DetermineOwnership(base_loc, neighbor_locs); // Determined by drawing and inspecting diagram block_ownership_t by_hand; @@ -218,7 +218,7 @@ TEST_CASE("Logical Location", "[Logical Location]") { THEN("We can find the ownership array of another block") { LogicalLocation base_loc(2, 1, 1, 1); - auto owns = DetermineOwnershipForest(base_loc, neighbor_locs); + auto owns = DetermineOwnership(base_loc, neighbor_locs); // Determined by drawing and inspecting diagram block_ownership_t by_hand; @@ -240,7 +240,7 @@ TEST_CASE("Logical Location", "[Logical Location]") { THEN("We can find the ownership array of yet another block") { LogicalLocation base_loc(2, 0, 0, 0); - auto owns = DetermineOwnershipForest(base_loc, neighbor_locs); + auto owns = DetermineOwnership(base_loc, neighbor_locs); // Determined by drawing and inspecting diagram, this should be the // ownership structure for every block in a uniform grid @@ -269,7 +269,7 @@ TEST_CASE("Logical Location", "[Logical Location]") { THEN("We can find the ownership array of yet another block") { LogicalLocation base_loc(3, 7, 7, 7); - auto owns = DetermineOwnershipForest(base_loc, neighbor_locs); + auto owns = DetermineOwnership(base_loc, neighbor_locs); // Determined by drawing and inspecting diagram, this is // the upper rightmost block in the grid on the finest refinement From d751b569323fda6f5d1bb0634f21751b2d84676b Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:35:05 -0600 Subject: [PATCH 19/39] rename for clarity --- src/mesh/mesh-amr_loadbalance.cpp | 2 +- src/mesh/mesh-gmg.cpp | 4 ++-- src/mesh/mesh.cpp | 2 +- src/mesh/mesh.hpp | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index 9e5776977927..f6d32b84b755 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -792,7 +792,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput BlockList_t new_block_list(nbe - nbs + 1); { // AMR Construct new MeshBlockList region PARTHENON_INSTRUMENT - RegionSize block_size = GetBlockSize(); + RegionSize block_size = GetDefaultBlockSize(); for (int n = nbs; n <= nbe; n++) { int on = newtoold[n]; diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 01c581a8dbdc..a13fa80c031c 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -102,7 +102,7 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app // See how many times we can go below logical level zero based on the // number of times a blocks zones can be reduced by 2^D int gmg_level_offset = std::numeric_limits::max(); - auto block_size_default = GetBlockSize(); + auto block_size_default = GetDefaultBlockSize(); for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (!mesh_size.symmetry(dir)) { int dir_allowed_levels = @@ -138,7 +138,7 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app // determines which process internal blocks live on auto loc = pmb->loc.GetParent(); while (loc.level() >= gmg_min_level && loc.morton() == pmb->loc.morton()) { - RegionSize block_size = GetBlockSize(); + RegionSize block_size = GetDefaultBlockSize(); BoundaryFlag block_bcs[6]; SetBlockSizeAndBoundaries(loc, block_size, block_bcs); gmg_block_lists[loc.level()].push_back( diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 8bc1d38f6213..f5d340429518 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1209,7 +1209,7 @@ std::int64_t Mesh::GetTotalCells() { int Mesh::GetNumberOfMeshBlockCells() const { return block_list.front()->GetNumberOfMeshBlockCells(); } -const RegionSize &Mesh::GetBlockSize() const { return base_block_size; } +const RegionSize &Mesh::GetDefaultBlockSize() const { return base_block_size; } const IndexShape &Mesh::GetLeafBlockCellBounds(CellLevel level) const { // TODO(JMM): Luke this is for your Metadata::fine stuff. diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 8e5225c16479..4026dccb064c 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -93,7 +93,7 @@ class Mesh { std::int64_t GetTotalCells(); // TODO(JMM): Move block_size into mesh. int GetNumberOfMeshBlockCells() const; - const RegionSize &GetBlockSize() const; + const RegionSize &GetDefaultBlockSize() const; RegionSize GetBlockSize(const LogicalLocation &loc) const; const IndexShape &GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const; From 9f2057439bc954c2249f058c45dea754521b4448 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:47:05 -0600 Subject: [PATCH 20/39] Split function into two --- src/mesh/mesh-amr_loadbalance.cpp | 7 +++++-- src/mesh/mesh-gmg.cpp | 13 +++++++++++-- src/mesh/mesh.cpp | 6 ++++-- src/mesh/mesh.hpp | 3 ++- 4 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index f6d32b84b755..a3d32f9eb0c9 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -927,7 +927,9 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput loclist = std::move(newloc); ranklist = std::move(newrank); costlist = std::move(newcost); - + + BuildGMGBlockLists(pin, app_in); + // Make sure all old sends/receives are done before we reconfigure the mesh #ifdef MPI_PARALLEL if (send_reqs.size() != 0) @@ -954,6 +956,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput // Thus we rebuild and synchronize the mesh now, but using a unique // neighbor precedence favoring the "old" fine blocks over "new" ones SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist, newly_refined); + SetGMGNeighbors(); BuildTagMapAndBoundaryBuffers(); std::string noncc = "mesh_internal_noncc"; for (int i = 0; i < DefaultNumPartitions(); ++i) { @@ -973,7 +976,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput // Rebuild just the ownership model, this time weighting the "new" fine blocks just // like any other blocks at their level. SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); - BuildGMGHierarchy(nbs, pin, app_in); + SetGMGNeighbors(); // Ownership does not impact anything about the buffers, so we don't need to // rebuild them if they were built above if (noncc_names.size() == 0) BuildTagMapAndBoundaryBuffers(); diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index a13fa80c031c..cb3a10c7a680 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -96,7 +96,7 @@ void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, } } -void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app_in) { +void Mesh::BuildGMGBlockLists(ParameterInput *pin, ApplicationInput *app_in) { if (!multigrid) return; // See how many times we can go below logical level zero based on the @@ -148,12 +148,20 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } } - // Sort the gmg block lists by gid and find neighbors + // Sort the gmg block lists by gid for (auto &[level, bl] : gmg_block_lists) { std::sort(bl.begin(), bl.end(), [](auto &a, auto &b){ return a->gid < b->gid;}); + } +} +void Mesh::SetGMGNeighbors() { + if (!multigrid) return; + const int gmg_min_level = GetGMGMinLevel(); + // Sort the gmg block lists by gid and find neighbors + for (auto &[level, bl] : gmg_block_lists) { for (auto &pmb : bl) { // Coarser neighbor + pmb->gmg_coarser_neighbors.clear(); if (pmb->loc.level() > gmg_min_level) { auto ploc = pmb->loc.GetParent(); int gid = forest.GetGid(ploc); @@ -165,6 +173,7 @@ void Mesh::BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app } // Finer neighbor(s) + pmb->gmg_finer_neighbors.clear(); if (pmb->loc.level() < current_level) { auto dlocs = pmb->loc.GetDaughters(ndim); for (auto &d : dlocs) { diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index f5d340429518..1c453a090a3a 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -426,8 +426,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages, MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, packages, resolved_packages, gflag); } + BuildGMGBlockLists(pin, app_in); SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); - BuildGMGHierarchy(nbs, pin, app_in); + SetGMGNeighbors(); ResetLoadBalanceVariables(); } @@ -683,8 +684,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr, MeshBlock::Make(i, i - nbs, loclist[i], block_size, block_bcs, this, pin, app_in, packages, resolved_packages, gflag, costlist[i]); } + BuildGMGBlockLists(pin, app_in); SetMeshBlockNeighbors(GridIdentifier::leaf(), block_list, ranklist); - BuildGMGHierarchy(nbs, pin, app_in); + SetGMGNeighbors(); ResetLoadBalanceVariables(); } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 4026dccb064c..3eb8575d424e 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -315,7 +315,8 @@ class Mesh { bool GatherCostListAndCheckBalance(); void RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput *app_in, int ntot); - void BuildGMGHierarchy(int nbs, ParameterInput *pin, ApplicationInput *app_in); + void BuildGMGBlockLists(ParameterInput *pin, ApplicationInput *app_in); + void SetGMGNeighbors(); void SetMeshBlockNeighbors(GridIdentifier grid_id, BlockList_t &block_list, const std::vector &ranklist, From 932966092bfbc44b7ee62e76b37079478b749c5d Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:52:11 -0600 Subject: [PATCH 21/39] More cleanup --- src/mesh/mesh.cpp | 14 ++------------ src/mesh/mesh.hpp | 4 ++-- 2 files changed, 4 insertions(+), 14 deletions(-) diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp index 1c453a090a3a..56079af8f700 100644 --- a/src/mesh/mesh.cpp +++ b/src/mesh/mesh.cpp @@ -1064,7 +1064,7 @@ void Mesh::FillDerived() { //---------------------------------------------------------------------------------------- // \!fn void Mesh::Initialize(bool init_problem, ParameterInput *pin) -// \brief initialization before the main loop as well as during remeshing +// \brief initialization before the main loop void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput *app_in) { PARTHENON_INSTRUMENT @@ -1192,26 +1192,16 @@ bool Mesh::SetBlockSizeAndBoundaries(LogicalLocation loc, RegionSize &block_size return valid_region; } -//---------------------------------------------------------------------------------------- -// \!fn void Mesh::GetBlockSize(const LogicalLocation &loc) const -// \brief Find the (hyper-)rectangular region of the grid covered by the block at -// logical location loc - -RegionSize Mesh::GetBlockSize(const LogicalLocation &loc) const { - // TODO(LFR): Update this - return forest.GetBlockDomain(loc); -} - std::int64_t Mesh::GetTotalCells() { auto &pmb = block_list.front(); return static_cast(nbtotal) * pmb->block_size.nx(X1DIR) * pmb->block_size.nx(X2DIR) * pmb->block_size.nx(X3DIR); } + // TODO(JMM): Move block_size into mesh. int Mesh::GetNumberOfMeshBlockCells() const { return block_list.front()->GetNumberOfMeshBlockCells(); } -const RegionSize &Mesh::GetDefaultBlockSize() const { return base_block_size; } const IndexShape &Mesh::GetLeafBlockCellBounds(CellLevel level) const { // TODO(JMM): Luke this is for your Metadata::fine stuff. diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 3eb8575d424e..8617e1eb1254 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -93,8 +93,8 @@ class Mesh { std::int64_t GetTotalCells(); // TODO(JMM): Move block_size into mesh. int GetNumberOfMeshBlockCells() const; - const RegionSize &GetDefaultBlockSize() const; - RegionSize GetBlockSize(const LogicalLocation &loc) const; + const RegionSize &GetDefaultBlockSize() const { return base_block_size; } + RegionSize GetBlockSize(const LogicalLocation &loc) const { return forest.GetBlockDomain(loc); } const IndexShape &GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const; const forest::Forest &Forest() const { return forest; } From 364d71a1285b17cda00366181f707b0e2ad7c8d4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:56:27 -0600 Subject: [PATCH 22/39] Remove RootGridInfo --- src/mesh/forest/logical_location.hpp | 14 -------------- src/mesh/mesh.hpp | 13 +------------ 2 files changed, 1 insertion(+), 26 deletions(-) diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 083690a3a858..e024f1360375 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -48,20 +48,6 @@ struct std::hash { }; namespace parthenon { - -// TODO(LFR): This can go away once MG is fixed for forests, and probably any routine that -// depends on it. -struct RootGridInfo { - int level; - std::array n; - std::array periodic; - // Defaults to root grid of single block at the - // coarsest level - RootGridInfo() : level(0), n{1, 1, 1}, periodic{false, false, false} {} - RootGridInfo(int level, int nx1, int nx2, int nx3, bool p1, bool p2, bool p3) - : level(level), n{nx1, nx2, nx3}, periodic{p1, p2, p3} {} -}; - //-------------------------------------------------------------------------------------- //! \struct LogicalLocation // \brief stores logical location and level of MeshBlock diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 8617e1eb1254..2cbd3a837c25 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -191,13 +191,6 @@ class Mesh { return forest.root_level + forest.forest_level; } - RootGridInfo GetRootGridInfo() const noexcept { - return RootGridInfo( - root_level, nrbx[0], nrbx[1], nrbx[2], - mesh_bcs[BoundaryFace::inner_x1] == BoundaryFlag::periodic && ndim > 0, - mesh_bcs[BoundaryFace::inner_x2] == BoundaryFlag::periodic && ndim > 1, - mesh_bcs[BoundaryFace::inner_x3] == BoundaryFlag::periodic && ndim > 2); - } int GetMaxLevel() const noexcept { return max_level; } int GetCurrentLevel() const noexcept { return current_level; } std::vector GetNbList() const noexcept { return nblist; } @@ -321,11 +314,7 @@ class Mesh { SetMeshBlockNeighbors(GridIdentifier grid_id, BlockList_t &block_list, const std::vector &ranklist, const std::unordered_set &newly_refined = {}); - void - SetSameLevelNeighbors(BlockList_t &block_list, const LogicalLocMap_t &loc_map, - RootGridInfo root_grid, int nbs, bool gmg_neighbors, - int composite_logical_level = 0, - const std::unordered_set &newly_refined = {}); + // defined in either the prob file or default_pgen.cpp in ../pgen/ static void InitUserMeshDataDefault(Mesh *mesh, ParameterInput *pin); std::function InitUserMeshData = From a4758a61ac06d54ec486d703221c679300b81f59 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 11:57:56 -0600 Subject: [PATCH 23/39] format --- src/basic_types.hpp | 6 ++- src/mesh/forest/block_ownership.cpp | 4 +- src/mesh/forest/block_ownership.hpp | 4 +- src/mesh/forest/forest.hpp | 4 +- src/mesh/forest/logical_location.cpp | 3 +- src/mesh/forest/logical_location.hpp | 2 +- src/mesh/forest/tree.cpp | 34 +++++++------- src/mesh/forest/tree.hpp | 6 ++- src/mesh/mesh-amr_loadbalance.cpp | 4 +- src/mesh/mesh-gmg.cpp | 69 ++++++++++++++-------------- src/mesh/mesh.hpp | 4 +- 11 files changed, 75 insertions(+), 65 deletions(-) diff --git a/src/basic_types.hpp b/src/basic_types.hpp index 1557c07f981b..567d04485827 100644 --- a/src/basic_types.hpp +++ b/src/basic_types.hpp @@ -71,8 +71,10 @@ struct GridIdentifier { GridType type = GridType::none; int logical_level = 0; - static GridIdentifier leaf() {return GridIdentifier{GridType::leaf, 0};} - static GridIdentifier two_level_composite(int level) {return GridIdentifier{GridType::two_level_composite, level};} + static GridIdentifier leaf() { return GridIdentifier{GridType::leaf, 0}; } + static GridIdentifier two_level_composite(int level) { + return GridIdentifier{GridType::two_level_composite, level}; + } }; constexpr bool IsSender(BoundaryType btype) { diff --git a/src/mesh/forest/block_ownership.cpp b/src/mesh/forest/block_ownership.cpp index f61093bfaaf2..70f6e27ef941 100644 --- a/src/mesh/forest/block_ownership.cpp +++ b/src/mesh/forest/block_ownership.cpp @@ -39,8 +39,8 @@ namespace parthenon { block_ownership_t DetermineOwnership(const LogicalLocation &main_block, - const std::vector &allowed_neighbors, - const std::unordered_set &newly_refined) { + const std::vector &allowed_neighbors, + const std::unordered_set &newly_refined) { block_ownership_t main_owns; auto ownership_level = [&](const LogicalLocation &a) { diff --git a/src/mesh/forest/block_ownership.hpp b/src/mesh/forest/block_ownership.hpp index 4ed9b0308c10..314db4438787 100644 --- a/src/mesh/forest/block_ownership.hpp +++ b/src/mesh/forest/block_ownership.hpp @@ -84,8 +84,8 @@ struct block_ownership_t { block_ownership_t DetermineOwnership(const LogicalLocation &main_block, - const std::vector &allowed_neighbors, - const std::unordered_set &newly_refined = {}); + const std::vector &allowed_neighbors, + const std::unordered_set &newly_refined = {}); // Given a topological element, ownership array of the sending block, and offset indices // defining the location of an index region within the block (i.e. the ghost zones passed diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index 52a5f4763d57..dbb239b0d716 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -84,7 +84,9 @@ class Forest { return trees.at(loc.tree())->FindNeighbors(loc, ox1, ox2, ox3); } - std::vector FindNeighbors(const LogicalLocation &loc, GridIdentifier grid_id = GridIdentifier::leaf()) const { + std::vector + FindNeighbors(const LogicalLocation &loc, + GridIdentifier grid_id = GridIdentifier::leaf()) const { return trees.at(loc.tree())->FindNeighbors(loc, grid_id); } std::size_t CountMeshBlock() const { diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index 7119f26321b0..bffa1110fa90 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -87,7 +87,7 @@ bool LogicalLocation::IsNeighbor(const LogicalLocation &in) const { } bool LogicalLocation::IsNeighborOfTE(const LogicalLocation &in, - const std::array &te_offset) const { + const std::array &te_offset) const { PARTHENON_REQUIRE(tree() == in.tree(), "Trying to compare locations not in the same octree."); const int max_level = std::max(in.level(), level()); @@ -135,5 +135,4 @@ std::vector LogicalLocation::GetDaughters(int ndim) const { return daughters; } - } // namespace parthenon diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index e024f1360375..78657b69c354 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -122,7 +122,7 @@ class LogicalLocation { // aggregate and POD type // volume bool IsNeighbor(const LogicalLocation &in) const; bool IsNeighborOfTE(const LogicalLocation &in, - const std::array &te_offset) const; + const std::array &te_offset) const; LogicalLocation GetSameLevelNeighbor(int ox1, int ox2, int ox3) const { return LogicalLocation(tree(), level(), lx1() + ox1, lx2() + ox2, lx3() + ox3); diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index ae3eb0ac867e..b1088ec88d1b 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -56,7 +56,7 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz } // Build in negative levels - for (int l = -20; l < 0; ++l) { + for (int l = -20; l < 0; ++l) { internal_nodes.emplace(LocMapEntry(LogicalLocation(my_id, l, 0, 0, 0), -1, -1)); } @@ -142,7 +142,8 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { return nadded; } -std::vector Tree::FindNeighbors(const LogicalLocation &loc, GridIdentifier grid_id) const { +std::vector Tree::FindNeighbors(const LogicalLocation &loc, + GridIdentifier grid_id) const { const Indexer3D offsets({ndim > 0 ? -1 : 0, ndim > 0 ? 1 : 0}, {ndim > 1 ? -1 : 0, ndim > 1 ? 1 : 0}, {ndim > 2 ? -1 : 0, ndim > 2 ? 1 : 0}); @@ -166,31 +167,32 @@ std::vector Tree::FindNeighbors(const LogicalLocation &loc, in } void Tree::FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int ox3, - std::vector *neighbor_locs, + std::vector *neighbor_locs, GridIdentifier grid_id) const { PARTHENON_REQUIRE( loc.tree() == my_id, "Trying to find neighbors in a tree with a LogicalLocation on a different tree."); - PARTHENON_REQUIRE((leaves.count(loc) == 1 || internal_nodes.count(loc) == 1), "Location must be in the tree to find neighbors."); + PARTHENON_REQUIRE((leaves.count(loc) == 1 || internal_nodes.count(loc) == 1), + "Location must be in the tree to find neighbors."); auto neigh = loc.GetSameLevelNeighbor(ox1, ox2, ox3); int n_idx = neigh.NeighborTreeIndex(); - bool include_same, include_fine, include_internal, include_coarse; + bool include_same, include_fine, include_internal, include_coarse; if (grid_id.type == GridType::leaf) { include_same = true; - include_fine = true; + include_fine = true; include_internal = false; include_coarse = true; - } else if (grid_id.type == GridType::two_level_composite) { - if (loc.level() == grid_id.logical_level) { + } else if (grid_id.type == GridType::two_level_composite) { + if (loc.level() == grid_id.logical_level) { include_same = true; - include_fine = false; + include_fine = false; include_internal = true; - include_coarse = true; - } else if (loc.level() == grid_id.logical_level - 1) { + include_coarse = true; + } else if (loc.level() == grid_id.logical_level - 1) { include_same = false; - include_fine = true; - include_internal = false; + include_fine = true; + include_internal = false; include_coarse = false; } else { PARTHENON_FAIL("Logic is wrong somewhere."); @@ -205,14 +207,14 @@ void Tree::FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int o if (neighbor_tree->leaves.count(tneigh) && include_same) { neighbor_locs->push_back({tneigh, orientation.TransformBack(tneigh, GetId())}); } else if (neighbor_tree->internal_nodes.count(tneigh)) { - if (include_fine) { + if (include_fine) { auto daughters = tneigh.GetDaughters(neighbor_tree->ndim); for (auto &n : daughters) { if (tloc.IsNeighbor(n)) neighbor_locs->push_back({n, orientation.TransformBack(n, GetId())}); } - } else if (include_internal) { - neighbor_locs->push_back({tneigh, orientation.TransformBack(tneigh, GetId())}); + } else if (include_internal) { + neighbor_locs->push_back({tneigh, orientation.TransformBack(tneigh, GetId())}); } } else if (neighbor_tree->leaves.count(tneigh.GetParent()) && include_coarse) { auto neighp = orientation.TransformBack(tneigh.GetParent(), GetId()); diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index ea5b3c286cb2..595b9af7b933 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -69,7 +69,9 @@ class Tree : public std::enable_shared_from_this { std::vector GetSortedInternalNodeList() const; RegionSize GetBlockDomain(const LogicalLocation &loc) const; std::array GetBlockBCs(const LogicalLocation &loc) const; - std::vector FindNeighbors(const LogicalLocation &loc, GridIdentifier grid_id = GridIdentifier::leaf()) const; + std::vector + FindNeighbors(const LogicalLocation &loc, + GridIdentifier grid_id = GridIdentifier::leaf()) const; std::vector FindNeighbors(const LogicalLocation &loc, int ox1, int ox2, int ox3) const; @@ -133,7 +135,7 @@ class Tree : public std::enable_shared_from_this { private: void FindNeighborsImpl(const LogicalLocation &loc, int ox1, int ox2, int ox3, - std::vector *neighbor_locs, + std::vector *neighbor_locs, GridIdentifier grid_type) const; int ndim; diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp index a3d32f9eb0c9..4afd80ffd85a 100644 --- a/src/mesh/mesh-amr_loadbalance.cpp +++ b/src/mesh/mesh-amr_loadbalance.cpp @@ -927,9 +927,9 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput loclist = std::move(newloc); ranklist = std::move(newrank); costlist = std::move(newcost); - + BuildGMGBlockLists(pin, app_in); - + // Make sure all old sends/receives are done before we reconfigure the mesh #ifdef MPI_PARALLEL if (send_reqs.size() != 0) diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index cb3a10c7a680..fc0c6da14b74 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -43,8 +43,8 @@ namespace parthenon { -void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, - BlockList_t &block_list, const std::vector &ranklist, +void Mesh::SetMeshBlockNeighbors( + GridIdentifier grid_id, BlockList_t &block_list, const std::vector &ranklist, const std::unordered_set &newly_refined) { Indexer3D offsets({ndim > 0 ? -1 : 0, ndim > 0 ? 1 : 0}, {ndim > 1 ? -1 : 0, ndim > 1 ? 1 : 0}, @@ -83,14 +83,14 @@ void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined); nb.ownership.initialized = true; } - + if (grid_id.type == GridType::leaf) { pmb->neighbors = all_neighbors; - } else if (grid_id.type == GridType::two_level_composite - && pmb->loc.level() == grid_id.logical_level) { + } else if (grid_id.type == GridType::two_level_composite && + pmb->loc.level() == grid_id.logical_level) { pmb->gmg_same_neighbors = all_neighbors; - } else if (grid_id.type == GridType::two_level_composite - && pmb->loc.level() == grid_id.logical_level - 1) { + } else if (grid_id.type == GridType::two_level_composite && + pmb->loc.level() == grid_id.logical_level - 1) { pmb->gmg_composite_finer_neighbors = all_neighbors; } } @@ -98,15 +98,14 @@ void Mesh::SetMeshBlockNeighbors(GridIdentifier grid_id, void Mesh::BuildGMGBlockLists(ParameterInput *pin, ApplicationInput *app_in) { if (!multigrid) return; - - // See how many times we can go below logical level zero based on the + + // See how many times we can go below logical level zero based on the // number of times a blocks zones can be reduced by 2^D int gmg_level_offset = std::numeric_limits::max(); auto block_size_default = GetDefaultBlockSize(); for (auto dir : {X1DIR, X2DIR, X3DIR}) { if (!mesh_size.symmetry(dir)) { - int dir_allowed_levels = - NumberOfBinaryTrailingZeros(block_size_default.nx(dir)); + int dir_allowed_levels = NumberOfBinaryTrailingZeros(block_size_default.nx(dir)); gmg_level_offset = std::min(dir_allowed_levels, gmg_level_offset); } } @@ -121,20 +120,20 @@ void Mesh::BuildGMGBlockLists(ParameterInput *pin, ApplicationInput *app_in) { // Create MeshData objects for GMG for (auto &[l, mdc] : gmg_mesh_data) mdc.SetMeshPointer(this); - - // Fill/create gmg block lists based on this ranks block list - for (auto &pmb : block_list) { - const int level = pmb->loc.level(); + + // Fill/create gmg block lists based on this ranks block list + for (auto &pmb : block_list) { + const int level = pmb->loc.level(); // Add the leaf block to its level gmg_block_lists[level].push_back(pmb); - + // Add the leaf block to the next finer level if required - if (level < current_level) { + if (level < current_level) { gmg_block_lists[level + 1].push_back(pmb); } - // Create internal blocks that share a Morton number with this block - // and add them to gmg two-level composite grid block lists. This + // Create internal blocks that share a Morton number with this block + // and add them to gmg two-level composite grid block lists. This // determines which process internal blocks live on auto loc = pmb->loc.GetParent(); while (loc.level() >= gmg_min_level && loc.morton() == pmb->loc.morton()) { @@ -142,15 +141,15 @@ void Mesh::BuildGMGBlockLists(ParameterInput *pin, ApplicationInput *app_in) { BoundaryFlag block_bcs[6]; SetBlockSizeAndBoundaries(loc, block_size, block_bcs); gmg_block_lists[loc.level()].push_back( - MeshBlock::Make(forest.GetGid(loc), -1, loc, block_size, block_bcs, this, pin, - app_in, packages, resolved_packages, gflag)); + MeshBlock::Make(forest.GetGid(loc), -1, loc, block_size, block_bcs, this, pin, + app_in, packages, resolved_packages, gflag)); loc = loc.GetParent(); } } // Sort the gmg block lists by gid - for (auto &[level, bl] : gmg_block_lists) { - std::sort(bl.begin(), bl.end(), [](auto &a, auto &b){ return a->gid < b->gid;}); + for (auto &[level, bl] : gmg_block_lists) { + std::sort(bl.begin(), bl.end(), [](auto &a, auto &b) { return a->gid < b->gid; }); } } @@ -158,34 +157,36 @@ void Mesh::SetGMGNeighbors() { if (!multigrid) return; const int gmg_min_level = GetGMGMinLevel(); // Sort the gmg block lists by gid and find neighbors - for (auto &[level, bl] : gmg_block_lists) { + for (auto &[level, bl] : gmg_block_lists) { for (auto &pmb : bl) { - // Coarser neighbor + // Coarser neighbor pmb->gmg_coarser_neighbors.clear(); if (pmb->loc.level() > gmg_min_level) { auto ploc = pmb->loc.GetParent(); int gid = forest.GetGid(ploc); - if (gid >= 0) { + if (gid >= 0) { int leaf_gid = forest.GetLeafGid(ploc); - pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, ploc, ranklist[leaf_gid], gid, - std::array{0, 0, 0}, 0, 0, 0, 0); + pmb->gmg_coarser_neighbors.emplace_back(pmb->pmy_mesh, ploc, ranklist[leaf_gid], + gid, std::array{0, 0, 0}, 0, 0, + 0, 0); } } - + // Finer neighbor(s) pmb->gmg_finer_neighbors.clear(); if (pmb->loc.level() < current_level) { auto dlocs = pmb->loc.GetDaughters(ndim); - for (auto &d : dlocs) { + for (auto &d : dlocs) { int gid = forest.GetGid(d); - if (gid >= 0) { + if (gid >= 0) { int leaf_gid = forest.GetLeafGid(d); - pmb->gmg_finer_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], gid, - std::array{0, 0, 0}, 0, 0, 0, 0); + pmb->gmg_finer_neighbors.emplace_back(pmb->pmy_mesh, d, ranklist[leaf_gid], + gid, std::array{0, 0, 0}, 0, 0, + 0, 0); } } } - + // Same level neighbors SetMeshBlockNeighbors(GridIdentifier::two_level_composite(level), bl, ranklist); } diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp index 2cbd3a837c25..5e7647ecb166 100644 --- a/src/mesh/mesh.hpp +++ b/src/mesh/mesh.hpp @@ -94,7 +94,9 @@ class Mesh { // TODO(JMM): Move block_size into mesh. int GetNumberOfMeshBlockCells() const; const RegionSize &GetDefaultBlockSize() const { return base_block_size; } - RegionSize GetBlockSize(const LogicalLocation &loc) const { return forest.GetBlockDomain(loc); } + RegionSize GetBlockSize(const LogicalLocation &loc) const { + return forest.GetBlockDomain(loc); + } const IndexShape &GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const; const forest::Forest &Forest() const { return forest; } From a5da7c92264ef395df3f7d9cf0b0c67b1ed5d0e8 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 12:05:21 -0600 Subject: [PATCH 24/39] remove unecessary location lists --- src/mesh/forest/tree.cpp | 20 +------------------- src/mesh/forest/tree.hpp | 3 --- 2 files changed, 1 insertion(+), 22 deletions(-) diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index b1088ec88d1b..b219dd05a81c 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -37,7 +37,7 @@ namespace forest { Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSize domain, std::array bcs) - : my_id(id), ndim(ndim), domain(domain), boundary_conditions(bcs), gmg_tlc_grids(50) { + : my_id(id), ndim(ndim), domain(domain), boundary_conditions(bcs) { // Add internal and leaf nodes of the initial tree for (int l = 0; l <= root_level; ++l) { for (int k = 0; k < (ndim > 2 ? (1LL << l) : 1); ++k) { @@ -49,7 +49,6 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz } else { internal_nodes.emplace(LocMapEntry(loc, -1, -1)); } - gmg_tlc_grids[l].emplace(LocMapEntry(loc, -1, -1)); } } } @@ -59,9 +58,6 @@ Tree::Tree(Tree::private_t, std::int64_t id, int ndim, int root_level, RegionSiz for (int l = -20; l < 0; ++l) { internal_nodes.emplace(LocMapEntry(LogicalLocation(my_id, l, 0, 0, 0), -1, -1)); } - - // Pre-populate the next finest two-level composite grid with coarser blocks - gmg_tlc_grids[root_level + 1] = gmg_tlc_grids[root_level]; } int Tree::AddMeshBlock(const LogicalLocation &loc, bool enforce_proper_nesting) { @@ -105,17 +101,6 @@ int Tree::Refine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) { } int nadded = daughters.size() - 1; - // Update multigrid levels - auto &gmg_grid_1 = gmg_tlc_grids[ref_loc.level() + 1]; - auto &gmg_grid_2 = gmg_tlc_grids[ref_loc.level() + 2]; - if (gmg_grid_1.count(ref_loc)) gmg_grid_1.erase(ref_loc); - for (auto &d : daughters) { - // Insert as fine leaf blocks on coarser two-level grid - gmg_grid_1.insert(LocMapEntry(d, gid_parent, -1)); - // Insert as coarse leaf blocks on finer two-level grid - gmg_grid_2.insert(LocMapEntry(d, gid_parent, -1)); - } - if (enforce_proper_nesting) { LogicalLocation parent = ref_loc.GetParent(); int ox1 = ref_loc.lx1() - (parent.lx1() << 1); @@ -268,16 +253,13 @@ int Tree::Derefine(const LogicalLocation &ref_loc, bool enforce_proper_nesting) } // Derefinement is ok - auto &gmg_grid = gmg_tlc_grids[ref_loc.level() + 1]; std::int64_t dgid = std::numeric_limits::max(); for (auto &d : daughters) { - gmg_grid.erase(d); auto node = leaves.extract(d); dgid = std::min(dgid, node.mapped().first); } internal_nodes.erase(ref_loc); leaves.insert(LocMapEntry(ref_loc, dgid, -1)); - gmg_grid.insert(LocMapEntry(ref_loc, dgid, -1)); return daughters.size() - 1; } diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 595b9af7b933..ba53c935d59d 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -151,9 +151,6 @@ class Tree : public std::enable_shared_from_this { LocMap_t leaves; LocMap_t internal_nodes; - // Two-level composite grids for geometric multigrid - std::vector gmg_tlc_grids; - // This contains all of the neighbor information for this tree, for each of the // 3^3 possible neighbor connections. Since an edge or node connection can have // multiple neighbors generally, we keep a map at each neighbor location from From c29147a12da209bba51e1fbad232c536940a46c1 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 12:17:33 -0600 Subject: [PATCH 25/39] update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f02eb296972c..2667057f04bd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ - [[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 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 From b2541875dbdf1189bbc7f300ef3fc5a77df49acd Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 12:23:43 -0600 Subject: [PATCH 26/39] update copyrights --- src/driver/driver.cpp | 2 +- src/interface/data_collection.cpp | 2 +- src/interface/mesh_data.hpp | 2 +- src/mesh/forest/block_ownership.cpp | 4 ++-- src/mesh/forest/cell_center_offsets.hpp | 2 +- src/mesh/forest/forest.hpp | 2 +- src/mesh/forest/logical_location.cpp | 4 ++-- src/mesh/forest/logical_location.hpp | 4 ++-- src/mesh/forest/relative_orientation.hpp | 2 +- src/mesh/forest/tree.hpp | 2 +- src/mesh/mesh-gmg.cpp | 4 ++-- tst/unit/test_logical_location.cpp | 2 +- 12 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp index 372ba30f0964..3b0e129ad2bd 100644 --- a/src/driver/driver.cpp +++ b/src/driver/driver.cpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/src/interface/data_collection.cpp b/src/interface/data_collection.cpp index 589c2c7128a6..275be3c7ff8f 100644 --- a/src/interface/data_collection.cpp +++ b/src/interface/data_collection.cpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/src/interface/mesh_data.hpp b/src/interface/mesh_data.hpp index 4df20bf26d88..e44e7def467e 100644 --- a/src/interface/mesh_data.hpp +++ b/src/interface/mesh_data.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/src/mesh/forest/block_ownership.cpp b/src/mesh/forest/block_ownership.cpp index 70f6e27ef941..a3a2f3f46c0b 100644 --- a/src/mesh/forest/block_ownership.cpp +++ b/src/mesh/forest/block_ownership.cpp @@ -4,10 +4,10 @@ // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== // Parthenon performance portable AMR framework -// Copyright(C) 2020-2023 The Parthenon collaboration +// Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/src/mesh/forest/cell_center_offsets.hpp b/src/mesh/forest/cell_center_offsets.hpp index a1237a514201..8641fe2f73cf 100644 --- a/src/mesh/forest/cell_center_offsets.hpp +++ b/src/mesh/forest/cell_center_offsets.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2023-2024. 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 diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index dbb239b0d716..4a9ebd5bf7ee 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2023-2024. 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 diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index bffa1110fa90..2db9c5856690 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -4,10 +4,10 @@ // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== // Parthenon performance portable AMR framework -// Copyright(C) 2020-2023 The Parthenon collaboration +// Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 78657b69c354..3a57f55dd3d5 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -4,10 +4,10 @@ // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== // Parthenon performance portable AMR framework -// Copyright(C) 2020-2023 The Parthenon collaboration +// Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/src/mesh/forest/relative_orientation.hpp b/src/mesh/forest/relative_orientation.hpp index 7bb711da47bc..f91aa3f5f10a 100644 --- a/src/mesh/forest/relative_orientation.hpp +++ b/src/mesh/forest/relative_orientation.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2024. 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 diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index ba53c935d59d..72cf337bd28b 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -1,5 +1,5 @@ //======================================================================================== -// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2024. 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 diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index fc0c6da14b74..c72c05b9f7f0 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -1,9 +1,9 @@ //======================================================================================== // Parthenon performance portable AMR framework -// Copyright(C) 2020-2023 The Parthenon collaboration +// Copyright(C) 2020-2024 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2020-2024. 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 diff --git a/tst/unit/test_logical_location.cpp b/tst/unit/test_logical_location.cpp index 665496d9012c..ddc6a18ece56 100644 --- a/tst/unit/test_logical_location.cpp +++ b/tst/unit/test_logical_location.cpp @@ -3,7 +3,7 @@ // Copyright(C) 2023 The Parthenon collaboration // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== -// (C) (or copyright) 2021-2023. Triad National Security, LLC. All rights reserved. +// (C) (or copyright) 2021-2024. 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 From 59bfdd15e36d601a7f1145d379526cc20b34c6fe Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 12:46:54 -0600 Subject: [PATCH 27/39] Move definitions out of header --- src/mesh/forest/forest.hpp | 3 ++ src/mesh/forest/tree.cpp | 41 ++++++++++++++++++++++++++ src/mesh/forest/tree.hpp | 59 +++++++------------------------------- 3 files changed, 54 insertions(+), 49 deletions(-) diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index 4a9ebd5bf7ee..5c90be03086e 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -135,6 +135,9 @@ class Forest { PARTHENON_REQUIRE(gids_resolved, "Asking for GID in invalid state."); return trees.at(loc.tree())->GetGid(loc); } + + // Get the gid of the leaf block with the same Morton number + // as loc (on the same tree) std::int64_t GetLeafGid(const LogicalLocation &loc) const { PARTHENON_REQUIRE(gids_resolved, "Asking for GID in invalid state."); return trees.at(loc.tree())->GetLeafGid(loc); diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index b219dd05a81c..cc7e9946cb43 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -330,5 +330,46 @@ void Tree::AddNeighborTree(CellCentOffsets offset, std::shared_ptr neighbo if (fidx >= 0) boundary_conditions[fidx] = BoundaryFlag::block; } +void Tree::InsertGid(const LogicalLocation &loc, std::int64_t gid) { + if (leaves.count(loc)) { + leaves[loc].second = leaves[loc].first; + leaves[loc].first = gid; + } else if (internal_nodes.count(loc)) { + internal_nodes[loc].second = internal_nodes[loc].first; + internal_nodes[loc].first = gid; + } else { + PARTHENON_FAIL("Tried to assign gid to non-existent block."); + } +} + +std::int64_t Tree::GetGid(const LogicalLocation &loc) const { + if (leaves.count(loc)) { + return leaves.at(loc).first; + } else if (internal_nodes.count(loc)) { + return internal_nodes.at(loc).first; + } + return -1; +} + +// Get the gid of the leaf block with the same Morton number +// as loc +std::int64_t Tree::GetLeafGid(const LogicalLocation &loc) const { + if (leaves.count(loc)) { + return leaves.at(loc).first; + } else if (internal_nodes.count(loc)) { + return GetLeafGid(loc.GetDaughter(0, 0, 0)); + } + return -1; +} + +std::int64_t Tree::GetOldGid(const LogicalLocation &loc) const { + if (leaves.count(loc)) { + return leaves.at(loc).second; + } else if (internal_nodes.count(loc)) { + return internal_nodes.at(loc).second; + } + return -1; +} + } // namespace forest } // namespace parthenon diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 72cf337bd28b..40130af7105f 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -74,62 +74,23 @@ class Tree : public std::enable_shared_from_this { GridIdentifier grid_id = GridIdentifier::leaf()) const; std::vector FindNeighbors(const LogicalLocation &loc, int ox1, int ox2, int ox3) const; - - std::vector - GetLocalLocationsFromNeighborLocation(const LogicalLocation &loc); - std::size_t CountMeshBlock() const { return leaves.size(); } + // Gid related methods + void InsertGid(const LogicalLocation &loc, std::int64_t gid); + std::int64_t GetGid(const LogicalLocation &loc) const; + std::int64_t GetOldGid(const LogicalLocation &loc) const; + // Get the gid of the leaf block with the same Morton number + // as loc + std::int64_t GetLeafGid(const LogicalLocation &loc) const; + // Methods for building tree connectivity void AddNeighborTree(CellCentOffsets offset, std::shared_ptr neighbor_tree, RelativeOrientation orient); - + + // Global id of the tree std::uint64_t GetId() const { return my_id; } - const std::unordered_map> & - GetLeaves() const { - return leaves; - } - - void InsertGid(const LogicalLocation &loc, std::int64_t gid) { - if (leaves.count(loc)) { - leaves[loc].second = leaves[loc].first; - leaves[loc].first = gid; - } else if (internal_nodes.count(loc)) { - internal_nodes[loc].second = internal_nodes[loc].first; - internal_nodes[loc].first = gid; - } else { - PARTHENON_FAIL("Tried to assign gid to non-existent block."); - } - } - - std::int64_t GetGid(const LogicalLocation &loc) const { - if (leaves.count(loc)) { - return leaves.at(loc).first; - } else if (internal_nodes.count(loc)) { - return internal_nodes.at(loc).first; - } - return -1; - } - - std::int64_t GetLeafGid(const LogicalLocation &loc) const { - if (leaves.count(loc)) { - return leaves.at(loc).first; - } else if (internal_nodes.count(loc)) { - return GetLeafGid(loc.GetDaughter(0, 0, 0)); - } - return -1; - } - - std::int64_t GetOldGid(const LogicalLocation &loc) const { - if (leaves.count(loc)) { - return leaves.at(loc).second; - } else if (internal_nodes.count(loc)) { - return internal_nodes.at(loc).second; - } - return -1; - } - // TODO(LFR): Eventually remove this. LogicalLocation athena_forest_loc; From a5c56d9ccae7b02c7753f6d4995b655bf0bc0038 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 27 Mar 2024 12:48:35 -0600 Subject: [PATCH 28/39] format --- src/mesh/forest/forest.hpp | 2 +- src/mesh/forest/tree.cpp | 2 +- src/mesh/forest/tree.hpp | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index 5c90be03086e..1fc6f65e54ee 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -136,7 +136,7 @@ class Forest { return trees.at(loc.tree())->GetGid(loc); } - // Get the gid of the leaf block with the same Morton number + // Get the gid of the leaf block with the same Morton number // as loc (on the same tree) std::int64_t GetLeafGid(const LogicalLocation &loc) const { PARTHENON_REQUIRE(gids_resolved, "Asking for GID in invalid state."); diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp index cc7e9946cb43..302216818af7 100644 --- a/src/mesh/forest/tree.cpp +++ b/src/mesh/forest/tree.cpp @@ -351,7 +351,7 @@ std::int64_t Tree::GetGid(const LogicalLocation &loc) const { return -1; } -// Get the gid of the leaf block with the same Morton number +// Get the gid of the leaf block with the same Morton number // as loc std::int64_t Tree::GetLeafGid(const LogicalLocation &loc) const { if (leaves.count(loc)) { diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index 40130af7105f..6c61d24c1ba5 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -76,18 +76,18 @@ class Tree : public std::enable_shared_from_this { int ox2, int ox3) const; std::size_t CountMeshBlock() const { return leaves.size(); } - // Gid related methods + // Gid related methods void InsertGid(const LogicalLocation &loc, std::int64_t gid); std::int64_t GetGid(const LogicalLocation &loc) const; std::int64_t GetOldGid(const LogicalLocation &loc) const; - // Get the gid of the leaf block with the same Morton number + // Get the gid of the leaf block with the same Morton number // as loc std::int64_t GetLeafGid(const LogicalLocation &loc) const; // Methods for building tree connectivity void AddNeighborTree(CellCentOffsets offset, std::shared_ptr neighbor_tree, RelativeOrientation orient); - + // Global id of the tree std::uint64_t GetId() const { return my_id; } From 51bf560d31b83fcf11c2e3e3252d01fa95263c68 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 28 Mar 2024 13:33:16 -0600 Subject: [PATCH 29/39] fix neighbor tree index finding --- src/mesh/forest/logical_location.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 3a57f55dd3d5..ecfe86ca3e6a 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -98,9 +98,9 @@ class LogicalLocation { // aggregate and POD type int NeighborTreeIndex() const { int up = 1LL << std::max(level(), 0); - int i1 = (l_[0] >= 0) - (l_[0] < (1LL << level())) + 1; - int i2 = (l_[1] >= 0) - (l_[1] < (1LL << level())) + 1; - int i3 = (l_[2] >= 0) - (l_[2] < (1LL << level())) + 1; + int i1 = (l_[0] >= 0) - (l_[0] < up) + 1; + int i2 = (l_[1] >= 0) - (l_[1] < up) + 1; + int i3 = (l_[2] >= 0) - (l_[2] < up) + 1; return i1 + 3 * i2 + 9 * i3; } From df545d226103507f885941c48de78a2efbe7b916 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 28 Mar 2024 14:51:05 -0600 Subject: [PATCH 30/39] cleanup --- src/mesh/forest/logical_location.cpp | 43 +++++++++++++++++ src/mesh/forest/logical_location.hpp | 70 +++++++--------------------- 2 files changed, 59 insertions(+), 54 deletions(-) diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index 2db9c5856690..7a5a61181b0a 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -35,6 +35,34 @@ #include "utils/morton_number.hpp" namespace parthenon { +std::string LogicalLocation::label() const { + return "([" + std::to_string(tree_idx_) + "] " + std::to_string(level_) + ": " + + std::to_string(l_[0]) + ", " + std::to_string(l_[1]) + ", " + + std::to_string(l_[2]) + ")"; +} + +bool LogicalLocation::IsInTree(int nghost) const { + const int low = -nghost; + const int up = (1LL << std::max(level(), 0)) + nghost; + return (l_[0] >= low) && (l_[0] < up) && (l_[1] >= low) && (l_[1] < up) && + (l_[2] >= low) && (l_[2] < up); +} + +int LogicalLocation::NeighborTreeIndex() const { + auto up = 1LL << std::max(level(), 0); + int i1 = (l_[0] >= 0) - (l_[0] < (1LL << up) + 1; + int i2 = (l_[1] >= 0) - (l_[1] < (1LL << up) + 1; + int i3 = (l_[2] >= 0) - (l_[2] < (1LL << up) + 1; + int idx = i1 + 3 * i2 + 9 * i3; + PARTHENON_REQUIRE(idx >= 0 && idx < 27, "Bad index."); + return idx; +} + +Real LogicalLocation::LLCoord(CoordinateDirection dir, BlockLocation bloc) const { + auto nblocks_tot = 1 << std::max(level(), 0); + return (static_cast(l(dir - 1)) + 0.5 * static_cast(bloc)) / + static_cast(nblocks_tot); +} bool LogicalLocation::IsContainedIn(const LogicalLocation &container) const { if (container.level() > level()) return false; @@ -135,4 +163,19 @@ std::vector LogicalLocation::GetDaughters(int ndim) const { return daughters; } +std::array +LogicalLocation::GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, int ox2, + int ox3) const { + // The neighbor block struct should only use the first two, but we have three to allow + // for this being a parent of neighbor, this should be checked for elsewhere + std::array f{0, 0, 0}; + if (neighbor.level() == level() + 1) { + int idx = 0; + if (ox1 == 0) f[idx++] = neighbor.lx1() % 2; + if (ox2 == 0) f[idx++] = neighbor.lx2() % 2; + if (ox3 == 0) f[idx++] = neighbor.lx3() % 2; + } + return f; +} + } // namespace parthenon diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index ecfe86ca3e6a..5fe4944e1355 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -36,17 +36,6 @@ #include "utils/error_checking.hpp" #include "utils/morton_number.hpp" -namespace parthenon { -class LogicalLocation; -} - -// This must be declared before an unordered_set of LogicalLocation is used -// below, but must be *implemented* after the class definition -template <> -struct std::hash { - std::size_t operator()(const parthenon::LogicalLocation &key) const noexcept; -}; - namespace parthenon { //-------------------------------------------------------------------------------------- //! \struct LogicalLocation @@ -71,11 +60,7 @@ class LogicalLocation { // aggregate and POD type : l_{l1, l2, l3}, level_{lev}, tree_idx_{tree}, morton_(lev, l1, l2, l3) {} LogicalLocation() : LogicalLocation(0, 0, 0, 0) {} - std::string label() const { - return "([" + std::to_string(tree_idx_) + "] " + std::to_string(level_) + ": " + - std::to_string(l_[0]) + ", " + std::to_string(l_[1]) + ", " + - std::to_string(l_[2]) + ")"; - } + std::string label() const; const auto &l(int i) const { return l_[i]; } const auto &lx1() const { return l_[0]; } const auto &lx2() const { return l_[1]; } @@ -86,31 +71,16 @@ class LogicalLocation { // aggregate and POD type // Check if this logical location is actually in the domain of the tree, // possibly including a ghost block halo around the tree - bool IsInTree(int nghost = 0) const { - const int low = -nghost; - const int up = (1LL << std::max(level(), 0)) + nghost; - return (l_[0] >= low) && (l_[0] < up) && (l_[1] >= low) && (l_[1] < up) && - (l_[2] >= low) && (l_[2] < up); - } + bool IsInTree(int nghost = 0) const; // Check if a LL is in the ghost block halo of the tree it is associated with bool IsInHalo(int nghost) const { return IsInTree(nghost) && !IsInTree(0); } - int NeighborTreeIndex() const { - int up = 1LL << std::max(level(), 0); - int i1 = (l_[0] >= 0) - (l_[0] < up) + 1; - int i2 = (l_[1] >= 0) - (l_[1] < up) + 1; - int i3 = (l_[2] >= 0) - (l_[2] < up) + 1; - return i1 + 3 * i2 + 9 * i3; - } + int NeighborTreeIndex() const; // Returns the coordinate in the range [0, 1] of the left side of // a logical location in a given direction on refinement level level - Real LLCoord(CoordinateDirection dir, BlockLocation bloc = BlockLocation::Left) const { - auto nblocks_tot = 1 << std::max(level(), 0); - return (static_cast(l(dir - 1)) + 0.5 * static_cast(bloc)) / - static_cast(nblocks_tot); - } + Real LLCoord(CoordinateDirection dir, BlockLocation bloc = BlockLocation::Left) const; bool IsContainedIn(const LogicalLocation &container) const; @@ -146,19 +116,8 @@ class LogicalLocation { // aggregate and POD type // 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 // correction is combined with boundary communication. - auto GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, int ox2, - int ox3) const { - // The neighbor block struct should only use the first two, but we have three to allow - // for this being a parent of neighbor, this should be checked for elsewhere - std::array f{0, 0, 0}; - if (neighbor.level() == level() + 1) { - int idx = 0; - if (ox1 == 0) f[idx++] = neighbor.lx1() % 2; - if (ox2 == 0) f[idx++] = neighbor.lx2() % 2; - if (ox3 == 0) f[idx++] = neighbor.lx3() % 2; - } - return f; - } + std::array GetAthenaXXFaceOffsets(const LogicalLocation &neighbor, int ox1, + int ox2, int ox3) const; }; inline bool operator<(const LogicalLocation &lhs, const LogicalLocation &rhs) { @@ -193,12 +152,15 @@ struct NeighborLocation { } // namespace parthenon -inline std::size_t std::hash::operator()( - const parthenon::LogicalLocation &key) const noexcept { - // TODO(LFR): Think more carefully about what the best choice for this key is, - // probably the least significant sizeof(size_t) * 8 bits of the morton number - // with 3 * (level - 21) trailing bits removed. - return key.morton().bits[0]; -} +// Inject hash function for LogicalLocation into the std namespace +template <> +struct std::hash { + std::size_t operator()(const parthenon::LogicalLocation &key) const noexcept { + // TODO(LFR): Think more carefully about what the best choice for this key is, + // probably the least significant sizeof(size_t) * 8 bits of the morton number + // with 3 * (level - 21) trailing bits removed. + return key.morton().bits[0]; + } +}; #endif // MESH_FOREST_LOGICAL_LOCATION_HPP_ From f1aba7a09734aab4b0b6fc86ddb685ea06f37c99 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 28 Mar 2024 15:02:12 -0600 Subject: [PATCH 31/39] not actually working --- src/mesh/forest/logical_location.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index 7a5a61181b0a..b2db34c9aeb1 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -50,9 +50,9 @@ bool LogicalLocation::IsInTree(int nghost) const { int LogicalLocation::NeighborTreeIndex() const { auto up = 1LL << std::max(level(), 0); - int i1 = (l_[0] >= 0) - (l_[0] < (1LL << up) + 1; - int i2 = (l_[1] >= 0) - (l_[1] < (1LL << up) + 1; - int i3 = (l_[2] >= 0) - (l_[2] < (1LL << up) + 1; + int i1 = (l_[0] >= 0) - (l_[0] < (1LL << up)) + 1; + int i2 = (l_[1] >= 0) - (l_[1] < (1LL << up)) + 1; + int i3 = (l_[2] >= 0) - (l_[2] < (1LL << up)) + 1; int idx = i1 + 3 * i2 + 9 * i3; PARTHENON_REQUIRE(idx >= 0 && idx < 27, "Bad index."); return idx; From 45b8b3d0a5e468caeff8056538f9dfa195079dc4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 12:50:43 -0600 Subject: [PATCH 32/39] format --- src/mesh/forest/tree.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp index e30c1e7feca8..7906215260fc 100644 --- a/src/mesh/forest/tree.hpp +++ b/src/mesh/forest/tree.hpp @@ -25,7 +25,6 @@ #include "basic_types.hpp" #include "defs.hpp" -#include "utils/cell_center_offsets.hpp" #include "mesh/forest/logical_location.hpp" #include "mesh/forest/relative_orientation.hpp" #include "utils/bit_hacks.hpp" From d4daeebe94f6b2a36e90dfd68f35753a30f0f0ea Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 12:58:05 -0600 Subject: [PATCH 33/39] remove dead code --- src/bvals/neighbor_block.cpp | 50 ------------------------------------ src/bvals/neighbor_block.hpp | 3 --- 2 files changed, 53 deletions(-) diff --git a/src/bvals/neighbor_block.cpp b/src/bvals/neighbor_block.cpp index 2d0c0b13a597..71f0958777b4 100644 --- a/src/bvals/neighbor_block.cpp +++ b/src/bvals/neighbor_block.cpp @@ -37,56 +37,6 @@ namespace parthenon { -//---------------------------------------------------------------------------------------- -// \!fn void NeighborBlock::SetNeighbor(int irank, int ilevel, int igid, int ilid, -// int iox1, int iox2, int iox3, NeighborConnect itype, -// int ibid, int itargetid, int ifi1=0, int ifi2=0) -// \brief Set neighbor information - -void NeighborBlock::SetNeighbor(LogicalLocation inloc, int irank, int ilevel, int igid, - int ilid, int iox1, int iox2, int iox3, - NeighborConnect itype, int ibid, int itargetid, - int ifi1, // =0 - int ifi2 // =0 -) { - snb.rank = irank; - snb.level = ilevel; - snb.gid = igid; - snb.lid = ilid; - ni.ox1 = iox1; - ni.ox2 = iox2; - ni.ox3 = iox3; - ni.type = itype; - ni.fi1 = ifi1; - ni.fi2 = ifi2; - bufid = ibid; - targetid = itargetid; - loc = inloc; - if (ni.type == NeighborConnect::face) { - if (ni.ox1 == -1) - fid = BoundaryFace::inner_x1; - else if (ni.ox1 == 1) - fid = BoundaryFace::outer_x1; - else if (ni.ox2 == -1) - fid = BoundaryFace::inner_x2; - else if (ni.ox2 == 1) - fid = BoundaryFace::outer_x2; - else if (ni.ox3 == -1) - fid = BoundaryFace::inner_x3; - else if (ni.ox3 == 1) - fid = BoundaryFace::outer_x3; - } - if (ni.type == NeighborConnect::edge) { - if (ni.ox3 == 0) - eid = ((((ni.ox1 + 1) >> 1) | ((ni.ox2 + 1) & 2))); - else if (ni.ox2 == 0) - eid = (4 + (((ni.ox1 + 1) >> 1) | ((ni.ox3 + 1) & 2))); - else if (ni.ox1 == 0) - eid = (8 + (((ni.ox2 + 1) >> 1) | ((ni.ox3 + 1) & 2))); - } - return; -} - NeighborConnect NCFromOffsets(const std::array offsets) { int connect_indicator = std::abs(offsets[0]) + std::abs(offsets[1]) + std::abs(offsets[2]); diff --git a/src/bvals/neighbor_block.hpp b/src/bvals/neighbor_block.hpp index 99b849985142..7e2b3cf1262b 100644 --- a/src/bvals/neighbor_block.hpp +++ b/src/bvals/neighbor_block.hpp @@ -146,9 +146,6 @@ struct NeighborBlock { block_ownership_t ownership; RegionSize block_size; - void SetNeighbor(LogicalLocation inloc, int irank, int ilevel, int igid, int ilid, - int iox1, int iox2, int iox3, NeighborConnect itype, int ibid, - int itargetid, int ifi1 = 0, int ifi2 = 0); NeighborBlock() = default; NeighborBlock(Mesh *mesh, LogicalLocation loc, int rank, int gid, int lid, std::array offsets, NeighborConnect type, int bid, int target_id, From c0624e9ffee7c2f0c1f6d9b343d44062dccbd718 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 16:02:43 -0600 Subject: [PATCH 34/39] fix merge bug --- src/mesh/forest/logical_location.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index b2db34c9aeb1..baa13dd9c3b3 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -50,9 +50,9 @@ bool LogicalLocation::IsInTree(int nghost) const { int LogicalLocation::NeighborTreeIndex() const { auto up = 1LL << std::max(level(), 0); - int i1 = (l_[0] >= 0) - (l_[0] < (1LL << up)) + 1; - int i2 = (l_[1] >= 0) - (l_[1] < (1LL << up)) + 1; - int i3 = (l_[2] >= 0) - (l_[2] < (1LL << up)) + 1; + int i1 = (l_[0] >= 0) - (l_[0] < up) + 1; + int i2 = (l_[1] >= 0) - (l_[1] < up) + 1; + int i3 = (l_[2] >= 0) - (l_[2] < up) + 1; int idx = i1 + 3 * i2 + 9 * i3; PARTHENON_REQUIRE(idx >= 0 && idx < 27, "Bad index."); return idx; From 42fc5cf440eb03c99a8729fb4640bf6058d595c6 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 16:04:23 -0600 Subject: [PATCH 35/39] correct morton numbers for negative levels --- src/mesh/forest/logical_location.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 5fe4944e1355..4ccabbe89123 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -54,10 +54,10 @@ class LogicalLocation { // aggregate and POD type // No check is provided that the requested LogicalLocation is in the allowed // range of logical location in the requested level. LogicalLocation(int lev, std::int64_t l1, std::int64_t l2, std::int64_t l3) - : l_{l1, l2, l3}, level_{lev}, tree_idx_{-1}, morton_(lev, l1, l2, l3) {} + : l_{l1, l2, l3}, level_{lev}, tree_idx_{-1}, morton_(std::max(lev, 0), l1, l2, l3) {} LogicalLocation(std::int64_t tree, int lev, std::int64_t l1, std::int64_t l2, std::int64_t l3) - : l_{l1, l2, l3}, level_{lev}, tree_idx_{tree}, morton_(lev, l1, l2, l3) {} + : l_{l1, l2, l3}, level_{lev}, tree_idx_{tree}, morton_(std::max(lev, 0), l1, l2, l3) {} LogicalLocation() : LogicalLocation(0, 0, 0, 0) {} std::string label() const; From a5515c26052ef7e4e5d7db6e255741696d773c08 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 16:29:02 -0600 Subject: [PATCH 36/39] More carefully deal with negative levels --- src/mesh/forest/logical_location.cpp | 21 +++++++++++++++++++++ src/mesh/forest/logical_location.hpp | 8 ++------ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index baa13dd9c3b3..9af59fcb913c 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -143,6 +143,27 @@ bool LogicalLocation::IsNeighborOfTE(const LogicalLocation &in, return neighbors; } +LogicalLocation LogicalLocation::GetParent(int nlevel) const { + // Shift all locations so that ghost tree locations have positive values, + // then shift back. Looks overly complicated to deal with ghosts and negative + // levels, but basically boils down to + // return LogicalLocation(tree(), level() - nlevel, lx1() >> nlevel, lx2() >> nlevel, lx3() >> nlevel); + // for most cases + const int norig = 1LL << std::max(level(), 0); + const int nparent = 1LL << std::max(level() - nlevel, 0); + constexpr int nmax_tree_offset = 5; + std::array lparent; + for (int dir = 0; dir < 3; ++dir) { + PARTHENON_DEBUG_REQUIRE((l(dir) >= -norig * nmax_tree_offset && l(dir) < (1 + nmax_tree_offset) * norig), "More than maximum number of tree offset."); + const int offset_l = l(dir) + nmax_tree_offset * norig; + lparent[dir] = (offset_l % norig) >> nlevel; + lparent[dir] += (offset_l / norig - nmax_tree_offset) * nparent; + } + + return LogicalLocation(tree(), level() - nlevel, lparent[0], lparent[1], lparent[2]); + } + + std::vector LogicalLocation::GetDaughters(int ndim) const { std::vector daughters; if (level() < 0) { diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 4ccabbe89123..13927351049c 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -98,16 +98,12 @@ class LogicalLocation { // aggregate and POD type return LogicalLocation(tree(), level(), lx1() + ox1, lx2() + ox2, lx3() + ox3); } - LogicalLocation GetParent(int nlevel = 1) const { - if (level() - nlevel < 0) return LogicalLocation(tree(), level() - nlevel, 0, 0, 0); - return LogicalLocation(tree(), level() - nlevel, lx1() >> nlevel, lx2() >> nlevel, - lx3() >> nlevel); - } + LogicalLocation GetParent(int nlevel = 1) const; std::vector GetDaughters(int ndim = 3) const; LogicalLocation GetDaughter(int ox1, int ox2, int ox3) const { - if (level() < 0) return LogicalLocation(tree(), level() + 1, 0, 0, 0); + if (level() < 0) return LogicalLocation(tree(), level() + 1, lx1(), lx2(), lx3()); return LogicalLocation(tree(), level() + 1, (lx1() << 1) + ox1, (lx2() << 1) + ox2, (lx3() << 1) + ox3); } From b767ba243aa46a75b24dedade65b7821aa9c8839 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 16:29:37 -0600 Subject: [PATCH 37/39] format --- src/mesh/forest/logical_location.cpp | 35 ++++++++++++++-------------- src/mesh/forest/logical_location.hpp | 6 +++-- 2 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/mesh/forest/logical_location.cpp b/src/mesh/forest/logical_location.cpp index 9af59fcb913c..7eb7382c2ef1 100644 --- a/src/mesh/forest/logical_location.cpp +++ b/src/mesh/forest/logical_location.cpp @@ -144,25 +144,26 @@ bool LogicalLocation::IsNeighborOfTE(const LogicalLocation &in, } LogicalLocation LogicalLocation::GetParent(int nlevel) const { - // Shift all locations so that ghost tree locations have positive values, - // then shift back. Looks overly complicated to deal with ghosts and negative - // levels, but basically boils down to - // return LogicalLocation(tree(), level() - nlevel, lx1() >> nlevel, lx2() >> nlevel, lx3() >> nlevel); - // for most cases - const int norig = 1LL << std::max(level(), 0); - const int nparent = 1LL << std::max(level() - nlevel, 0); - constexpr int nmax_tree_offset = 5; - std::array lparent; - for (int dir = 0; dir < 3; ++dir) { - PARTHENON_DEBUG_REQUIRE((l(dir) >= -norig * nmax_tree_offset && l(dir) < (1 + nmax_tree_offset) * norig), "More than maximum number of tree offset."); - const int offset_l = l(dir) + nmax_tree_offset * norig; - lparent[dir] = (offset_l % norig) >> nlevel; - lparent[dir] += (offset_l / norig - nmax_tree_offset) * nparent; - } - - return LogicalLocation(tree(), level() - nlevel, lparent[0], lparent[1], lparent[2]); + // Shift all locations so that ghost tree locations have positive values, + // then shift back. Looks overly complicated to deal with ghosts and negative + // levels, but basically boils down to + // return LogicalLocation(tree(), level() - nlevel, lx1() >> nlevel, lx2() >> nlevel, + // lx3() >> nlevel); for most cases + const int norig = 1LL << std::max(level(), 0); + const int nparent = 1LL << std::max(level() - nlevel, 0); + constexpr int nmax_tree_offset = 5; + std::array lparent; + for (int dir = 0; dir < 3; ++dir) { + PARTHENON_DEBUG_REQUIRE( + (l(dir) >= -norig * nmax_tree_offset && l(dir) < (1 + nmax_tree_offset) * norig), + "More than maximum number of tree offset."); + const int offset_l = l(dir) + nmax_tree_offset * norig; + lparent[dir] = (offset_l % norig) >> nlevel; + lparent[dir] += (offset_l / norig - nmax_tree_offset) * nparent; } + return LogicalLocation(tree(), level() - nlevel, lparent[0], lparent[1], lparent[2]); +} std::vector LogicalLocation::GetDaughters(int ndim) const { std::vector daughters; diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index 13927351049c..1ced6b035cb2 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -54,10 +54,12 @@ class LogicalLocation { // aggregate and POD type // No check is provided that the requested LogicalLocation is in the allowed // range of logical location in the requested level. LogicalLocation(int lev, std::int64_t l1, std::int64_t l2, std::int64_t l3) - : l_{l1, l2, l3}, level_{lev}, tree_idx_{-1}, morton_(std::max(lev, 0), l1, l2, l3) {} + : l_{l1, l2, l3}, level_{lev}, tree_idx_{-1}, + morton_(std::max(lev, 0), l1, l2, l3) {} LogicalLocation(std::int64_t tree, int lev, std::int64_t l1, std::int64_t l2, std::int64_t l3) - : l_{l1, l2, l3}, level_{lev}, tree_idx_{tree}, morton_(std::max(lev, 0), l1, l2, l3) {} + : l_{l1, l2, l3}, level_{lev}, tree_idx_{tree}, + morton_(std::max(lev, 0), l1, l2, l3) {} LogicalLocation() : LogicalLocation(0, 0, 0, 0) {} std::string label() const; From a9060c0aeb8b34cecbef4e6f84ece6f3d03bdb5d Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 2 Apr 2024 11:40:47 -0600 Subject: [PATCH 38/39] add back gmg tests --- .github/workflows/ci-extended.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci-extended.yml b/.github/workflows/ci-extended.yml index 44b2ad9d8b48..3d216ba5c716 100644 --- a/.github/workflows/ci-extended.yml +++ b/.github/workflows/ci-extended.yml @@ -60,7 +60,7 @@ jobs: cd build # Pick GPU with most available memory export CUDA_VISIBLE_DEVICES=$(nvidia-smi --query-gpu=memory.free,index --format=csv,nounits,noheader | sort -nr | head -1 | awk '{ print $NF }') - ctest -L performance -LE perf-reg -E gmg + ctest -L performance -LE perf-reg # run regression tests - name: Regression tests @@ -68,7 +68,7 @@ jobs: cd build # Pick GPU with most available memory export CUDA_VISIBLE_DEVICES=$(nvidia-smi --query-gpu=memory.free,index --format=csv,nounits,noheader | sort -nr | head -1 | awk '{ print $NF }') - ctest -L regression -L ${{ matrix.parallel }} -LE perf-reg -E gmg --timeout 3600 + ctest -L regression -L ${{ matrix.parallel }} -LE perf-reg --timeout 3600 # Test Ascent integration (only most complex setup with MPI and on device) - name: Ascent tests From b9ab262ff4dfefb38dc641089981a4911ec48bca Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 10 Apr 2024 13:28:31 -0600 Subject: [PATCH 39/39] fix possible gmg ownership bug --- src/mesh/mesh-gmg.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index c72c05b9f7f0..e5ec8a6c1ca4 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -77,7 +77,7 @@ void Mesh::SetMeshBlockNeighbors( // Set neighbor block ownership auto &nb = all_neighbors.back(); - auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc); + auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc, grid_id); nb.ownership = DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined);