From 1dafe264a0f43c4f5171963482417b7bb5a0ddb4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 20 Mar 2024 15:37:17 -0600 Subject: [PATCH 1/7] 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 2/7] 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 3/7] 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 4/7] 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 5/7] 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 45b8b3d0a5e468caeff8056538f9dfa195079dc4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Mon, 1 Apr 2024 12:50:43 -0600 Subject: [PATCH 6/7] 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 7/7] 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,