From 8aaa28856a9302531955f5085c8949632ce62dd4 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 10:31:20 -0600 Subject: [PATCH 01/20] store simple boundary information in nodes --- src/mesh/forest/forest.hpp | 2 ++ src/mesh/forest/forest_node.hpp | 4 +++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp index db9de2e0025b..7f4d93931ae9 100644 --- a/src/mesh/forest/forest.hpp +++ b/src/mesh/forest/forest.hpp @@ -68,6 +68,8 @@ class ForestDefinition { PARTHENON_REQUIRE(periodic_connection, "Must specify another edge for periodic boundary conditions."); bc_edges.emplace_back(ForestBC{edge, bf, periodic_connection}); + for (auto &pnode : edge.nodes) + pnode->on_physical_boundary = true; } void AddInitialRefinement(const LogicalLocation &loc) { diff --git a/src/mesh/forest/forest_node.hpp b/src/mesh/forest/forest_node.hpp index 57249a89e7ea..9fbfaf3e4deb 100644 --- a/src/mesh/forest/forest_node.hpp +++ b/src/mesh/forest/forest_node.hpp @@ -32,7 +32,8 @@ namespace forest { class Face; class Node { public: - Node(int id_in, std::array pos) : id(id_in), x(pos), associated_faces{} {} + Node(int id_in, std::array pos) + : id(id_in), x(pos), associated_faces{}, on_physical_boundary{false} {} static std::shared_ptr create(int id, std::array pos) { return std::make_shared(id, pos); @@ -41,6 +42,7 @@ class Node { std::uint32_t id; std::array x; std::unordered_set> associated_faces; + bool on_physical_boundary; }; } // namespace forest } // namespace parthenon From 32ece06d309293ab72b5a621a224b1173016f5a3 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 10:32:55 -0600 Subject: [PATCH 02/20] add tree block selector --- src/pack/block_selector.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/pack/block_selector.hpp b/src/pack/block_selector.hpp index 1283ecffd573..8bcfd2e97fc4 100644 --- a/src/pack/block_selector.hpp +++ b/src/pack/block_selector.hpp @@ -74,6 +74,12 @@ inline block_selector_func_t OnPhysicalBoundary() { }; } +inline block_selector_func_t OnTree(std::size_t tree_id) { + return [tree_id](MeshBlockData *pmbd) { + return pmbd->GetBlockPointer()->loc.tree() == tree_id; + }; +} + inline block_selector_func_t FineOnCompositeGrid(const MeshData *pmd) { if (pmd->grid.type == GridType::two_level_composite) { const int fine_level = pmd->grid.logical_level; From 21949ed9a84b5da1553f09e693a2d0e4c883c2bc Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 10:33:30 -0600 Subject: [PATCH 03/20] Add transfinite interpolation to boundary comm example --- .../boundary_exchange/boundary_exchange.cpp | 97 ++++++++++++++++++- .../boundary_exchange/boundary_exchange.hpp | 9 ++ .../boundary_exchange_driver.cpp | 73 ++++++++++---- 3 files changed, 158 insertions(+), 21 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp index 14f966485ca5..1c943030e7d1 100644 --- a/example/boundary_exchange/boundary_exchange.cpp +++ b/example/boundary_exchange/boundary_exchange.cpp @@ -12,6 +12,7 @@ //======================================================================================== // Standard Includes +#include #include #include #include @@ -101,6 +102,97 @@ TaskStatus SetBlockValues(MeshData *md) { return TaskStatus::complete; } +struct ParameterizedLine { + const Real x1, y1; + const Real x2, y2; + + using node = parthenon::forest::Node; + ParameterizedLine(std::shared_ptr start, std::shared_ptr end) + : x1{start->x[0]}, y1{start->x[1]}, x2{end->x[0]}, y2{end->x[1]} {} + + KOKKOS_INLINE_FUNCTION + Real GetX(Real u) const { return x1 * (1.0 - u) + x2 * u; } + + KOKKOS_INLINE_FUNCTION + Real GetY(Real u) const { return y1 * (1.0 - u) + y2 * u; } +}; + +TaskStatus SetCoordinates(MeshData *md) { + using TE = parthenon::TopologicalElement; + auto pmesh = md->GetMeshPointer(); + auto desc = parthenon::MakePackDescriptor(md); + + for (auto &ptree : pmesh->forest.GetTrees()) { + auto tree_id = ptree->GetId(); + auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId())); + printf("Tree: %i\n", ptree->GetId()); + int i{0}; + Real posx[4], posy[4]; + for (auto &pnode : ptree->forest_nodes) { + posx[i] = pnode->x[0]; + posy[i] = pnode->x[1]; + printf(" Node %i:(%e, %e)\n", i++, pnode->x[0], pnode->x[1]); + } + + auto &pnodes = ptree->forest_nodes; + ParameterizedLine c1(pnodes[0], pnodes[1]); + ParameterizedLine c3(pnodes[2], pnodes[3]); + ParameterizedLine c2(pnodes[0], pnodes[2]); + ParameterizedLine c4(pnodes[1], pnodes[3]); + + IndexRange ib = md->GetBoundsI(IndexDomain::interior, TE::NN); + IndexRange jb = md->GetBoundsJ(IndexDomain::interior, TE::NN); + IndexRange kb = md->GetBoundsK(IndexDomain::interior, TE::NN); + + parthenon::par_for( + parthenon::loop_pattern_mdrange_tag, "SetPosition", DevExecSpace(), 0, + pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + auto &coords = pack.GetCoordinates(b); + const Real u = coords.X(k, j, i); + const Real v = coords.X(k, j, i); + const Real mu = 1.0 - u; + const Real mv = 1.0 - v; + + Real x = mv * c1.GetX(u) + v * c3.GetX(u) + mu * c2.GetX(v) + u * c4.GetX(v) - + (mu * mv * c1.GetX(0.0) + u * mv * c1.GetX(1.0) + + mu * v * c2.GetX(1.0) + u * v * c3.GetX(1.0)); + Real y = mv * c1.GetY(u) + v * c3.GetY(u) + mu * c2.GetY(v) + u * c4.GetY(v) - + (mu * mv * c1.GetY(0.0) + u * mv * c1.GetY(1.0) + + mu * v * c2.GetY(1.0) + u * v * c3.GetY(1.0)); + pack(b, TE::NN, position(0), k, j, i) = x; // x-position + pack(b, TE::NN, position(1), k, j, i) = y; // y-position + }); + } + return TaskStatus::complete; +} + +TaskStatus FixTrivalentNodes(MeshData *md) { + using TE = parthenon::TopologicalElement; + auto pmesh = md->GetMeshPointer(); + auto desc = parthenon::MakePackDescriptor(md); + + IndexRange ib_in = md->GetBoundsI(IndexDomain::interior, TE::NN); + IndexRange jb_in = md->GetBoundsJ(IndexDomain::interior, TE::NN); + IndexRange kb_in = md->GetBoundsK(IndexDomain::interior, TE::NN); + + for (auto &ptree : pmesh->forest.GetTrees()) { + auto tree_id = ptree->GetId(); + auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId())); + printf("Tree: %i\n", ptree->GetId()); + int pos{0}; + for (auto &pnode : ptree->forest_nodes) { + bool trivalent = + (pnode->associated_faces.size() == 3) && !pnode->on_physical_boundary; + if (trivalent) { + printf(" Node %i is trivalent.\n"); + } + pos++; + } + } + return TaskStatus::complete; +} + std::shared_ptr Initialize(ParameterInput *pin) { auto package = std::make_shared("boundary_exchange"); Params ¶ms = package->AllParams(); @@ -109,8 +201,11 @@ std::shared_ptr Initialize(ParameterInput *pin) { std::vector{8}); m.RegisterRefinementOps(); - package->AddField(neighbor_info::name(), m); + package->AddField(m); + Metadata m_node({Metadata::Node, Metadata::Independent, Metadata::FillGhost}, + std::vector{2}); + package->AddField(m_node); return package; } diff --git a/example/boundary_exchange/boundary_exchange.hpp b/example/boundary_exchange/boundary_exchange.hpp index 7b0cbb87bd88..bb83391a61c7 100644 --- a/example/boundary_exchange/boundary_exchange.hpp +++ b/example/boundary_exchange/boundary_exchange.hpp @@ -34,7 +34,16 @@ struct neighbor_info : public parthenon::variable_names::base_t { static std::string name() { return "neighbor_info"; } }; +struct position : public parthenon::variable_names::base_t { + template + KOKKOS_INLINE_FUNCTION position(Ts &&...args) + : parthenon::variable_names::base_t(std::forward(args)...) {} + static std::string name() { return "position"; } +}; + TaskStatus SetBlockValues(MeshData *rc); +TaskStatus SetCoordinates(MeshData *rc); +TaskStatus FixTrivalentNodes(MeshData *rc); std::shared_ptr Initialize(ParameterInput *pin); } // namespace boundary_exchange diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp index 640705ca77ef..f3bfd98525df 100644 --- a/example/boundary_exchange/boundary_exchange_driver.cpp +++ b/example/boundary_exchange/boundary_exchange_driver.cpp @@ -34,24 +34,7 @@ using boundary_exchange::BoundaryExchangeDriver; Packages_t ProcessPackages(std::unique_ptr &pin); -int main(int argc, char *argv[]) { - ParthenonManager pman; - - pman.app_input->ProcessPackages = ProcessPackages; - - // This is called on each mesh block whenever the mesh changes. - // pman.app_input->InitMeshBlockUserData = &calculate_pi::SetInOrOutBlock; - - auto manager_status = pman.ParthenonInitEnv(argc, argv); - if (manager_status == ParthenonStatus::complete) { - pman.ParthenonFinalize(); - return 0; - } - if (manager_status == ParthenonStatus::error) { - pman.ParthenonFinalize(); - return 1; - } - +parthenon::forest::ForestDefinition MakeFourTreeForestWithOneTreeRotated() { // Create the nodes for the forest, the x-y positions are only used for // visualizing the forest configuration and *do not* determine the global // coordinates of the trees. @@ -95,7 +78,55 @@ int main(int argc, char *argv[]) { forest_def.AddBC(edge_t({n[7], n[8]})); forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, 1, 0, 0, 0)); - pman.ParthenonInitPackagesAndMesh(forest_def); + + return forest_def; +} + +parthenon::forest::ForestDefinition n_blocks(int nblocks) { + using namespace parthenon::forest; + std::unordered_map> nodes; + ForestDefinition forest_def; + parthenon::Real xoffset = 0.0; + + for (int point = 0; point < 2 * nblocks; ++point) { + nodes[point] = Node::create(point, {std::sin(point * M_PI / nblocks) + xoffset, + std::cos(point * M_PI / nblocks)}); + } + + using edge_t = parthenon::forest::Edge; + for (int point = 0; point < 2 * nblocks; ++point) { + forest_def.AddBC( + edge_t({nodes[point % (2 * nblocks)], nodes[(point + 1) % (2 * nblocks)]})); + } + + nodes[2 * nblocks] = Node::create(2 * nblocks, {0.0 + xoffset, 0.0}); + auto &n = nodes; + for (int t = 0; t < nblocks; ++t) + forest_def.AddFace( + t, {n[2 * t + 1], n[2 * t], n[(2 * t + 2) % (2 * nblocks)], n[2 * nblocks]}); + + return forest_def; +} + +int main(int argc, char *argv[]) { + ParthenonManager pman; + + pman.app_input->ProcessPackages = ProcessPackages; + + // This is called on each mesh block whenever the mesh changes. + // pman.app_input->InitMeshBlockUserData = &calculate_pi::SetInOrOutBlock; + + auto manager_status = pman.ParthenonInitEnv(argc, argv); + if (manager_status == ParthenonStatus::complete) { + pman.ParthenonFinalize(); + return 0; + } + if (manager_status == ParthenonStatus::error) { + pman.ParthenonFinalize(); + return 1; + } + + pman.ParthenonInitPackagesAndMesh(n_blocks(3)); // This needs to be scoped so that the driver object is destructed before Finalize { @@ -150,7 +181,9 @@ TaskCollection BoundaryExchangeDriver::MakeTaskCollection(T &blocks) { auto &md = pmesh->mesh_data.GetOrAdd("base", i); TaskID none(0); auto fill = tl.AddTask(none, SetBlockValues, md.get()); - auto bound = AddBoundaryExchangeTasks(fill, tl, md, true); + auto set_coords = tl.AddTask(none, SetCoordinates, md.get()); + auto bound = AddBoundaryExchangeTasks(fill | set_coords, tl, md, true); + auto fix = tl.AddTask(bound, FixTrivalentNodes, md.get()); } } From ca89c0f079b4c43d20596b694ce6782d4903dcdf Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:25:05 -0600 Subject: [PATCH 04/20] convenience --- src/utils/cell_center_offsets.hpp | 5 ++--- src/utils/indexer.hpp | 9 +++++++++ 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/utils/cell_center_offsets.hpp b/src/utils/cell_center_offsets.hpp index 2bc41e323f7e..41d95148fb20 100644 --- a/src/utils/cell_center_offsets.hpp +++ b/src/utils/cell_center_offsets.hpp @@ -26,7 +26,6 @@ #include "basic_types.hpp" #include "defs.hpp" #include "utils/bit_hacks.hpp" -#include "utils/indexer.hpp" namespace parthenon { // CellCentOffsets defines the position of a topological element @@ -39,7 +38,7 @@ namespace parthenon { // TODO(LFR): Consider switching this to C-style enum within a namespace to avoid // static_cast -enum class Offset : int { Low = -1, Middle = 0, Up = 1 }; +enum Offset : int { Low = -1, Middle = 0, Up = 1 }; inline int operator+(Offset a, int b) { return static_cast(a) + b; } inline int operator+(int b, Offset a) { return static_cast(a) + b; } inline Offset operator-(Offset in) { return static_cast(-static_cast(in)); } @@ -58,7 +57,7 @@ struct CellCentOffsets { Offset &operator[](int idx) { return u[idx]; } const Offset &operator[](int idx) const { return u[idx]; } - int operator()(CoordinateDirection dir) const { return static_cast(u[dir - 1]); } + Offset operator()(CoordinateDirection dir) const { return u[dir - 1]; } operator std::array() const { return {static_cast(u[0]), static_cast(u[1]), static_cast(u[2])}; diff --git a/src/utils/indexer.hpp b/src/utils/indexer.hpp index f996b9c89281..a452c0beafe4 100644 --- a/src/utils/indexer.hpp +++ b/src/utils/indexer.hpp @@ -19,6 +19,7 @@ #include #include +#include "utils/cell_center_offsets.hpp" #include "utils/concepts_lite.hpp" #include "utils/utils.hpp" @@ -34,6 +35,14 @@ struct block_ownership_t { bool &operator()(int ox1, int ox2, int ox3) { return ownership[ox1 + 1][ox2 + 1][ox3 + 1]; } + KOKKOS_FORCEINLINE_FUNCTION + bool &operator()(const CellCentOffsets &offset) { + return ownership[offset[0]+1][offset[1]+1][offset[2]+1]; + } + KOKKOS_FORCEINLINE_FUNCTION + const bool &operator()(const CellCentOffsets &offset) const { + return ownership[offset[0]+1][offset[1]+1][offset[2]+1]; + } KOKKOS_FORCEINLINE_FUNCTION block_ownership_t() : block_ownership_t(false) {} From 93239a091426fcca821ac3bdb077aaa73e248b8e Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:25:24 -0600 Subject: [PATCH 05/20] more convenience --- src/mesh/forest/logical_location.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/mesh/forest/logical_location.hpp b/src/mesh/forest/logical_location.hpp index e1417d5d49ac..7a07fce55ecb 100644 --- a/src/mesh/forest/logical_location.hpp +++ b/src/mesh/forest/logical_location.hpp @@ -34,6 +34,7 @@ #include "basic_types.hpp" #include "defs.hpp" +#include "utils/cell_center_offsets.hpp" #include "utils/error_checking.hpp" #include "utils/morton_number.hpp" @@ -123,6 +124,10 @@ class LogicalLocation { // aggregate and POD type return {(lx1() & 1LL) == 1LL, (lx2() & 1LL) == 1LL, (lx3() & 1LL) == 1LL}; } + bool IsOnTreeBoundary(CellCentOffsets offset) const { + return IsOnTreeBoundary(offset(X1DIR), offset(X2DIR), offset(X3DIR)); + } + bool IsOnTreeBoundary(int ox1, int ox2, int ox3) const { const int nup = (1 << std::max(level(), 0)) - 1; const bool bound1 = From 68969697aba56a09de75620f5ec54d4286f21965 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:25:39 -0600 Subject: [PATCH 06/20] Add some more transformation types --- .../logical_coordinate_transformation.cpp | 35 +++++++++++++++++++ .../logical_coordinate_transformation.hpp | 12 +++++-- 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/src/mesh/forest/logical_coordinate_transformation.cpp b/src/mesh/forest/logical_coordinate_transformation.cpp index 34d4bec27a6f..c6d97bdeaba1 100644 --- a/src/mesh/forest/logical_coordinate_transformation.cpp +++ b/src/mesh/forest/logical_coordinate_transformation.cpp @@ -100,6 +100,41 @@ CellCentOffsets LogicalCoordinateTransformation::Transform(CellCentOffsets in) c return out; } +CellCentOffsets LogicalCoordinateTransformation::InverseTransform(CellCentOffsets in) const { + CellCentOffsets out; + for (int dir = 0; dir < 3; ++dir) { + const int indir = abs(dir_connection[dir]); + out.u[dir] = dir_flip[dir] ? -in.u[indir] : in.u[indir]; + } + return out; +} + +block_ownership_t LogicalCoordinateTransformation::Transform(const block_ownership_t &in) const { + block_ownership_t out; + for (int ox3 : {-1 , 0, 1}) { + for (int ox2 : {-1 , 0, 1}) { + for (int ox1 : {-1 , 0, 1}) { + auto offset = CellCentOffsets(ox1, ox2, ox3); + out(Transform(offset)) = in(offset); + } + } + } + return out; +} + +block_ownership_t LogicalCoordinateTransformation::InverseTransform(const block_ownership_t &in) const { + block_ownership_t out; + for (int ox3 : {-1 , 0, 1}) { + for (int ox2 : {-1 , 0, 1}) { + for (int ox1 : {-1 , 0, 1}) { + auto offset = CellCentOffsets(ox1, ox2, ox3); + out(offset) = in(InverseTransform(offset)); + } + } + } + return out; +} + LogicalCoordinateTransformation ComposeTransformations(const LogicalCoordinateTransformation &first, const LogicalCoordinateTransformation &second) { diff --git a/src/mesh/forest/logical_coordinate_transformation.hpp b/src/mesh/forest/logical_coordinate_transformation.hpp index 0e0d0bc293cf..8faf16ab1dcd 100644 --- a/src/mesh/forest/logical_coordinate_transformation.hpp +++ b/src/mesh/forest/logical_coordinate_transformation.hpp @@ -51,6 +51,9 @@ struct LogicalCoordinateTransformation { LogicalLocation InverseTransform(const LogicalLocation &loc_in, std::int64_t origin) const; CellCentOffsets Transform(CellCentOffsets in) const; + CellCentOffsets InverseTransform(CellCentOffsets in) const; + block_ownership_t Transform(const block_ownership_t &in) const; + block_ownership_t InverseTransform(const block_ownership_t &in) const; KOKKOS_INLINE_FUNCTION std::tuple Transform(TopologicalElement el) const { @@ -112,9 +115,12 @@ struct NeighborLocation { NeighborLocation(const LogicalLocation &g, const LogicalLocation &o, const LogicalCoordinateTransformation &lcoord_trans) : global_loc(g), origin_loc(o), lcoord_trans(lcoord_trans) {} - LogicalLocation global_loc; // Global location of neighboring block - LogicalLocation - origin_loc; // Logical location of neighboring block in index space of origin block + // Global location of neighboring block + LogicalLocation global_loc; + // Logical location of neighboring block in index space of origin block + LogicalLocation origin_loc; + // Coordinate transform to coords of neighbor tree from origin tree + // (i.e. global_loc = lcoord_trans.Transform(origin_loc)) LogicalCoordinateTransformation lcoord_trans; }; From f652ed548d40182845c03ab50f1d313a61ae2941 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:25:57 -0600 Subject: [PATCH 07/20] transform ownership correctly --- 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 791449aa7acd..33e65c77588d 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -81,7 +81,7 @@ void Mesh::SetMeshBlockNeighbors( auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc, grid_id); nb.ownership = - DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined); + nloc.lcoord_trans.InverseTransform(DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined)); nb.ownership.initialized = true; // Set logical coordinate transformation from this block to the neighbor From 26096368186daafb1a553fef7a9aef435940c9cd Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:26:31 -0600 Subject: [PATCH 08/20] Allow tree selector to include offsets --- src/pack/block_selector.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/pack/block_selector.hpp b/src/pack/block_selector.hpp index 8bcfd2e97fc4..71b25f1cf59a 100644 --- a/src/pack/block_selector.hpp +++ b/src/pack/block_selector.hpp @@ -32,6 +32,7 @@ #include "interface/variable.hpp" #include "pack/pack_utils.hpp" #include "pack/sparse_pack_base.hpp" +#include "utils/cell_center_offsets.hpp" #include "utils/concepts_lite.hpp" #include "utils/type_list.hpp" #include "utils/utils.hpp" @@ -74,9 +75,11 @@ inline block_selector_func_t OnPhysicalBoundary() { }; } -inline block_selector_func_t OnTree(std::size_t tree_id) { - return [tree_id](MeshBlockData *pmbd) { - return pmbd->GetBlockPointer()->loc.tree() == tree_id; +inline block_selector_func_t OnTree(std::size_t tree_id, + CellCentOffsets offset = CellCentOffsets(0, 0, 0)) { + return [tree_id, offset](MeshBlockData *pmbd) { + const auto &loc = pmbd->GetBlockPointer()->loc; + return (loc.tree() == tree_id) && loc.IsOnTreeBoundary(offset); }; } From cb88d7aa7f638f265cf686e6db62ece792715c20 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:26:48 -0600 Subject: [PATCH 09/20] work with new overloads --- src/bvals/comms/boundary_communication.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp index 78121cd3fac9..4376bf007cd8 100644 --- a/src/bvals/comms/boundary_communication.cpp +++ b/src/bvals/comms/boundary_communication.cpp @@ -303,7 +303,7 @@ TaskStatus SetBounds(std::shared_ptr> &md) { Kokkos::parallel_for( Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) { const auto [il, jl, kl] = - lcoord_trans.InverseTransform({ii + m, jj, kk}); + lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); if (idxer.IsActive(kl, jl, il)) var(iel, tt, uu, vv, kl, jl, il) = fac * buf[m]; }); @@ -323,7 +323,7 @@ TaskStatus SetBounds(std::shared_ptr> &md) { Kokkos::parallel_for( Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) { const auto [il, jl, kl] = - lcoord_trans.InverseTransform({ii + m, jj, kk}); + lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); if (idxer.IsActive(kl, jl, il)) var(iel, tt, uu, vv, kl, jl, il) = default_val; }); From 3e774e5cb02fceedfdc9838d4239148e601258a0 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Wed, 7 May 2025 16:26:57 -0600 Subject: [PATCH 10/20] working version --- .../boundary_exchange/boundary_exchange.cpp | 26 ++++++++++++++++--- .../boundary_exchange/boundary_exchange.hpp | 2 +- .../boundary_exchange_driver.cpp | 15 ++++++++--- 3 files changed, 34 insertions(+), 9 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp index 1c943030e7d1..f6e4e681a807 100644 --- a/example/boundary_exchange/boundary_exchange.cpp +++ b/example/boundary_exchange/boundary_exchange.cpp @@ -167,25 +167,43 @@ TaskStatus SetCoordinates(MeshData *md) { return TaskStatus::complete; } -TaskStatus FixTrivalentNodes(MeshData *md) { +TaskStatus FixTrivalentNodes2D(MeshData *md) { using TE = parthenon::TopologicalElement; auto pmesh = md->GetMeshPointer(); auto desc = parthenon::MakePackDescriptor(md); IndexRange ib_in = md->GetBoundsI(IndexDomain::interior, TE::NN); IndexRange jb_in = md->GetBoundsJ(IndexDomain::interior, TE::NN); - IndexRange kb_in = md->GetBoundsK(IndexDomain::interior, TE::NN); + IndexRange kb = md->GetBoundsK(IndexDomain::interior, TE::NN); for (auto &ptree : pmesh->forest.GetTrees()) { auto tree_id = ptree->GetId(); - auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId())); + printf("Tree: %i\n", ptree->GetId()); int pos{0}; for (auto &pnode : ptree->forest_nodes) { bool trivalent = (pnode->associated_faces.size() == 3) && !pnode->on_physical_boundary; if (trivalent) { - printf(" Node %i is trivalent.\n"); + parthenon::CellCentOffsets offset(2 * (pos % 2) - 1, 2 * (pos / 2) - 1, 0); + auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId(), offset)); + printf(" Node %i is trivalent and we picked out %i blocks.\n", pos, pack.GetNBlocks()); + + const int icorner = (offset(X1DIR) == parthenon::Offset::Low) ? ib_in.s : ib_in.e; + const int jcorner = (offset(X2DIR) == parthenon::Offset::Low) ? jb_in.s : jb_in.e; + IndexRange ib = offset(X1DIR) == parthenon::Offset::Low ? IndexRange{0, ib_in.s - 1} : IndexRange{ib_in.e + 1, ib_in.e + parthenon::Globals::nghost}; + IndexRange jb = offset(X2DIR) == parthenon::Offset::Low ? IndexRange{0, jb_in.s - 1} : IndexRange{jb_in.e + 1, ib_in.e + parthenon::Globals::nghost}; + parthenon::par_for( + parthenon::loop_pattern_mdrange_tag, "SetPosition", DevExecSpace(), 0, + pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, + KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { + const int enclosed_area = std::abs((i - icorner) * (j - jcorner)); + int ioff{1}; + while (ioff * ioff < enclosed_area) ioff++; + ioff *= offset(X1DIR); + pack(b, TE::NN, position(0), k, j, i) = pack(b, TE::NN, position(0), k, jcorner, icorner + ioff); // x-position + pack(b, TE::NN, position(1), k, j, i) = pack(b, TE::NN, position(1), k, jcorner, icorner + ioff); // y-position + }); } pos++; } diff --git a/example/boundary_exchange/boundary_exchange.hpp b/example/boundary_exchange/boundary_exchange.hpp index bb83391a61c7..f8853e00be63 100644 --- a/example/boundary_exchange/boundary_exchange.hpp +++ b/example/boundary_exchange/boundary_exchange.hpp @@ -43,7 +43,7 @@ struct position : public parthenon::variable_names::base_t { TaskStatus SetBlockValues(MeshData *rc); TaskStatus SetCoordinates(MeshData *rc); -TaskStatus FixTrivalentNodes(MeshData *rc); +TaskStatus FixTrivalentNodes2D(MeshData *rc); std::shared_ptr Initialize(ParameterInput *pin); } // namespace boundary_exchange diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp index f3bfd98525df..d924abcedb2e 100644 --- a/example/boundary_exchange/boundary_exchange_driver.cpp +++ b/example/boundary_exchange/boundary_exchange_driver.cpp @@ -86,11 +86,12 @@ parthenon::forest::ForestDefinition n_blocks(int nblocks) { using namespace parthenon::forest; std::unordered_map> nodes; ForestDefinition forest_def; - parthenon::Real xoffset = 0.0; + parthenon::Real xoffset = 0.1; + parthenon::Real yoffset = 0.1; for (int point = 0; point < 2 * nblocks; ++point) { nodes[point] = Node::create(point, {std::sin(point * M_PI / nblocks) + xoffset, - std::cos(point * M_PI / nblocks)}); + std::cos(point * M_PI / nblocks) + yoffset}); } using edge_t = parthenon::forest::Edge; @@ -99,11 +100,17 @@ parthenon::forest::ForestDefinition n_blocks(int nblocks) { edge_t({nodes[point % (2 * nblocks)], nodes[(point + 1) % (2 * nblocks)]})); } - nodes[2 * nblocks] = Node::create(2 * nblocks, {0.0 + xoffset, 0.0}); + nodes[2 * nblocks] = Node::create(2 * nblocks, {xoffset, yoffset}); auto &n = nodes; for (int t = 0; t < nblocks; ++t) forest_def.AddFace( t, {n[2 * t + 1], n[2 * t], n[(2 * t + 2) % (2 * nblocks)], n[2 * nblocks]}); + + int level = 1; + printf("loc: %i\n", (1 << level) - 1); + forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, level, (1 << level) - 1, (1 << level) - 1, 0)); + forest_def.AddInitialRefinement(parthenon::LogicalLocation(1, level, (1 << level) - 1, (1 << level) - 1, 0)); + forest_def.AddInitialRefinement(parthenon::LogicalLocation(2, level, (1 << level) - 1, (1 << level) - 1, 0)); return forest_def; } @@ -183,7 +190,7 @@ TaskCollection BoundaryExchangeDriver::MakeTaskCollection(T &blocks) { auto fill = tl.AddTask(none, SetBlockValues, md.get()); auto set_coords = tl.AddTask(none, SetCoordinates, md.get()); auto bound = AddBoundaryExchangeTasks(fill | set_coords, tl, md, true); - auto fix = tl.AddTask(bound, FixTrivalentNodes, md.get()); + auto fix = tl.AddTask(bound, FixTrivalentNodes2D, md.get()); } } From 39b9a80a3ba228f77c9099c7a72b738354a4c429 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:15:17 -0600 Subject: [PATCH 11/20] fix masking to be before transformation --- src/bvals/comms/boundary_communication.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp index 4376bf007cd8..8a6d43d8f06a 100644 --- a/src/bvals/comms/boundary_communication.cpp +++ b/src/bvals/comms/boundary_communication.cpp @@ -302,10 +302,11 @@ TaskStatus SetBounds(std::shared_ptr> &md) { const int ii = i; Kokkos::parallel_for( Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) { - const auto [il, jl, kl] = - lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); - if (idxer.IsActive(kl, jl, il)) + if (idxer.IsActive(kk, jj, ii + m)) { + const auto [il, jl, kl] = + lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); var(iel, tt, uu, vv, kl, jl, il) = fac * buf[m]; + } }); }); } else if (bnd_info(b).allocated && bound_type != BoundaryType::flxcor_recv) { @@ -322,10 +323,11 @@ TaskStatus SetBounds(std::shared_ptr> &md) { const int ii = i; Kokkos::parallel_for( Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) { - const auto [il, jl, kl] = - lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); - if (idxer.IsActive(kl, jl, il)) + if (idxer.IsActive(kk, jj, ii + m)) { + const auto [il, jl, kl] = + lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); var(iel, tt, uu, vv, kl, jl, il) = default_val; + } }); }); } From f77cb3f1ffdd025441e5fcc8c79472293a061cc6 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:49:31 -0600 Subject: [PATCH 12/20] add equality operator --- src/utils/cell_center_offsets.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/utils/cell_center_offsets.hpp b/src/utils/cell_center_offsets.hpp index 41d95148fb20..b41b625d85b3 100644 --- a/src/utils/cell_center_offsets.hpp +++ b/src/utils/cell_center_offsets.hpp @@ -73,6 +73,10 @@ struct CellCentOffsets { // (in cyclic order, XY, YZ, ZX, XYZ) along with the offset of the // element in that direction from the cell center. std::vector> GetNormals() const; + + bool operator==(const CellCentOffsets& other) const { + return u == other.u; + } bool IsNode() const { return 3 == abs(static_cast(u[0])) + abs(static_cast(u[1])) + From 4c932269288f42939ff79abc6dd779511d1d46cd Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:49:51 -0600 Subject: [PATCH 13/20] add cubed circle definition --- .../boundary_exchange_driver.cpp | 36 +++++++++++++++++-- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp index d924abcedb2e..17325c157c56 100644 --- a/example/boundary_exchange/boundary_exchange_driver.cpp +++ b/example/boundary_exchange/boundary_exchange_driver.cpp @@ -90,8 +90,8 @@ parthenon::forest::ForestDefinition n_blocks(int nblocks) { parthenon::Real yoffset = 0.1; for (int point = 0; point < 2 * nblocks; ++point) { - nodes[point] = Node::create(point, {std::sin(point * M_PI / nblocks) + xoffset, - std::cos(point * M_PI / nblocks) + yoffset}); + nodes[point] = Node::create(point, {-std::sin(point * M_PI / nblocks) + xoffset, + -std::cos(point * M_PI / nblocks) + yoffset}); } using edge_t = parthenon::forest::Edge; @@ -115,6 +115,35 @@ parthenon::forest::ForestDefinition n_blocks(int nblocks) { return forest_def; } +parthenon::forest::ForestDefinition cubed_circle(std::vector radii) { + using namespace parthenon::forest; + std::vector> nodes_x, nodes_y, nodes_c; + std::size_t node_idx{0}; + std::shared_ptr origin = Node::create(node_idx++, {0.0, 0.0}); + double fac = 1.0 / sqrt(2.0); + for (auto rad : radii) { + nodes_x.push_back(Node::create(node_idx++, {rad, 0.0})); + nodes_c.push_back(Node::create(node_idx++, {rad * fac, rad * fac})); + nodes_y.push_back(Node::create(node_idx++, {0.0, rad})); + } + + ForestDefinition forest_def; + std::size_t face_idx{0}; + forest_def.AddFace(face_idx++, {origin, nodes_x[0], nodes_y[0], nodes_c[0]}); + forest_def.AddBC(Edge({origin, nodes_x[0]})); + forest_def.AddBC(Edge({origin, nodes_y[0]})); + for (int i = 0; i < nodes_x.size() - 1; ++i) { + forest_def.AddFace(face_idx++, {nodes_x[i], nodes_x[i + 1], nodes_c[i], nodes_c[i + 1]}); + forest_def.AddFace(face_idx++, {nodes_c[i], nodes_c[i + 1], nodes_y[i], nodes_y[i + 1]}); + forest_def.AddBC(Edge({nodes_x[i], nodes_x[i + 1]})); + forest_def.AddBC(Edge({nodes_y[i], nodes_y[i + 1]})); + } + forest_def.AddBC(Edge({nodes_x.back(), nodes_c.back()})); + forest_def.AddBC(Edge({nodes_c.back(), nodes_y.back()})); + + return forest_def; +} + int main(int argc, char *argv[]) { ParthenonManager pman; @@ -133,7 +162,8 @@ int main(int argc, char *argv[]) { return 1; } - pman.ParthenonInitPackagesAndMesh(n_blocks(3)); + pman.ParthenonInitPackagesAndMesh(cubed_circle({1.0, 2.0, 3.0})); + //pman.ParthenonInitPackagesAndMesh(n_blocks(3)); // This needs to be scoped so that the driver object is destructed before Finalize { From 72f55de5a6c90d91b4c03384b8b6d079ecc51059 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:50:30 -0600 Subject: [PATCH 14/20] add janky arc for transfinite interpolation --- .../boundary_exchange/boundary_exchange.cpp | 33 ++++++++++++++++--- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp index f6e4e681a807..9305b7531dae 100644 --- a/example/boundary_exchange/boundary_exchange.cpp +++ b/example/boundary_exchange/boundary_exchange.cpp @@ -103,18 +103,43 @@ TaskStatus SetBlockValues(MeshData *md) { } struct ParameterizedLine { + enum class ltype {straight, arc}; + const Real x1, y1; const Real x2, y2; - + Real r, delta, phi; + ltype type; + using node = parthenon::forest::Node; ParameterizedLine(std::shared_ptr start, std::shared_ptr end) - : x1{start->x[0]}, y1{start->x[1]}, x2{end->x[0]}, y2{end->x[1]} {} + : x1{start->x[0]}, y1{start->x[1]}, x2{end->x[0]}, y2{end->x[1]}, type{ltype::straight} { + Real d1 = std::sqrt(x1 * x1 + y1 * y1); + Real d2 = std::sqrt(x2 * x2 + y2 * y2); + if (std::abs(d1 - d2) < 1.e-8 && d1 > 1.1) { + type = ltype::arc; + r = d1; + delta = M_PI / 4.0; + phi = 0.0; + if (y1 > 1.e-5) phi = delta; + } + } KOKKOS_INLINE_FUNCTION - Real GetX(Real u) const { return x1 * (1.0 - u) + x2 * u; } + Real GetX(Real u) const { + if (type == ltype::straight) + return x1 * (1.0 - u) + x2 * u; + else + return r * cos(delta * u + phi); + } + KOKKOS_INLINE_FUNCTION - Real GetY(Real u) const { return y1 * (1.0 - u) + y2 * u; } + Real GetY(Real u) const { + if (type == ltype::straight) + return y1 * (1.0 - u) + y2 * u; + else + return r * sin(delta * u + phi); + } }; TaskStatus SetCoordinates(MeshData *md) { From 3729fc524028d85bf61eb7633674c6a45a7f722f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:51:24 -0600 Subject: [PATCH 15/20] get all of the frames right for neighbor ownership --- src/bvals/comms/bnd_info.cpp | 6 +++++- src/mesh/forest/logical_coordinate_transformation.cpp | 2 +- src/mesh/mesh-gmg.cpp | 10 +++++----- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/bvals/comms/bnd_info.cpp b/src/bvals/comms/bnd_info.cpp index 2b4948278c38..4467927ae849 100644 --- a/src/bvals/comms/bnd_info.cpp +++ b/src/bvals/comms/bnd_info.cpp @@ -244,7 +244,11 @@ CalcIndices(const NeighborBlock &nb, MeshBlock *pmb, if (sox2 == 0) sox2 = loc.l(1) % 2 == 1 ? 1 : -1; if (sox3 == 0) sox3 = loc.l(2) % 2 == 1 ? 1 : -1; } - owns = GetIndexRangeMaskFromOwnership(el, nb.ownership, sox1, sox2, sox3); + // We are working in the frame of the origin block (i.e. the receiving block), + // so we use the ownership array in that frame to get the masking initially. + // This mask must then be transformed to the frame of the neighbor block (i.e. + // sending block) since the index range we are masking is defined in that frame. + owns = lcoord_trans.Transform(GetIndexRangeMaskFromOwnership(el, nb.origin_ownership, sox1, sox2, sox3)); } return SpatiallyMaskedIndexer6D(owns, {0, tensor_shape[0] - 1}, {0, tensor_shape[1] - 1}, {0, tensor_shape[2] - 1}, diff --git a/src/mesh/forest/logical_coordinate_transformation.cpp b/src/mesh/forest/logical_coordinate_transformation.cpp index c6d97bdeaba1..7f2f075f0987 100644 --- a/src/mesh/forest/logical_coordinate_transformation.cpp +++ b/src/mesh/forest/logical_coordinate_transformation.cpp @@ -128,7 +128,7 @@ block_ownership_t LogicalCoordinateTransformation::InverseTransform(const block_ for (int ox2 : {-1 , 0, 1}) { for (int ox1 : {-1 , 0, 1}) { auto offset = CellCentOffsets(ox1, ox2, ox3); - out(offset) = in(InverseTransform(offset)); + out(offset) = in(Transform(offset)); } } } diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 33e65c77588d..965fa9ca61dd 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -76,14 +76,14 @@ void Mesh::SetMeshBlockNeighbors( all_neighbors.emplace_back(pmb->pmy_mesh, nloc.global_loc, nloc.origin_loc, ranklist[lgid], gid, offsets, bid, tid, f[0], f[1]); - // Set neighbor block ownership + // Set neighbor block ownership (note that this is ownership in coordinates of the neighbor block) auto &nb = all_neighbors.back(); auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc, grid_id); - - nb.ownership = - nloc.lcoord_trans.InverseTransform(DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined)); + nb.ownership = DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined); nb.ownership.initialized = true; - + nb.origin_ownership = nloc.lcoord_trans.InverseTransform(nb.ownership); + nb.origin_ownership.initialized = true; + // Set logical coordinate transformation from this block to the neighbor nb.lcoord_trans = nloc.lcoord_trans; } From 7922d12225deda4fc1c88c90b53176f6ce96a7f7 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:53:03 -0600 Subject: [PATCH 16/20] add origin ownership to neighbor block structure --- src/bvals/neighbor_block.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/bvals/neighbor_block.hpp b/src/bvals/neighbor_block.hpp index 78f954b04386..821373e4571d 100644 --- a/src/bvals/neighbor_block.hpp +++ b/src/bvals/neighbor_block.hpp @@ -63,7 +63,10 @@ struct NeighborBlock { // Offset of the neighbor block relative to origin block CellCentOffsets offsets; // Ownership of neighbor block of different topological elements + // (in the frame of the neighbor block) block_ownership_t ownership; + // (in the frame of the origin block) + block_ownership_t origin_ownership; // Logical coordinate transformation from the main block to this neighbor forest::LogicalCoordinateTransformation lcoord_trans; From 5d4edd47f95d10654b4bf67416cfdda5f4e337c0 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:53:20 -0600 Subject: [PATCH 17/20] add trivalent node correction --- .../boundary_exchange/boundary_exchange.cpp | 51 ++++++++++++++----- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp index 9305b7531dae..d9c09f6aa543 100644 --- a/example/boundary_exchange/boundary_exchange.cpp +++ b/example/boundary_exchange/boundary_exchange.cpp @@ -150,13 +150,11 @@ TaskStatus SetCoordinates(MeshData *md) { for (auto &ptree : pmesh->forest.GetTrees()) { auto tree_id = ptree->GetId(); auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId())); - printf("Tree: %i\n", ptree->GetId()); int i{0}; Real posx[4], posy[4]; for (auto &pnode : ptree->forest_nodes) { posx[i] = pnode->x[0]; posy[i] = pnode->x[1]; - printf(" Node %i:(%e, %e)\n", i++, pnode->x[0], pnode->x[1]); } auto &pnodes = ptree->forest_nodes; @@ -203,8 +201,6 @@ TaskStatus FixTrivalentNodes2D(MeshData *md) { for (auto &ptree : pmesh->forest.GetTrees()) { auto tree_id = ptree->GetId(); - - printf("Tree: %i\n", ptree->GetId()); int pos{0}; for (auto &pnode : ptree->forest_nodes) { bool trivalent = @@ -212,22 +208,53 @@ TaskStatus FixTrivalentNodes2D(MeshData *md) { if (trivalent) { parthenon::CellCentOffsets offset(2 * (pos % 2) - 1, 2 * (pos / 2) - 1, 0); auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId(), offset)); - printf(" Node %i is trivalent and we picked out %i blocks.\n", pos, pack.GetNBlocks()); + if (pack.GetNBlocks() == 0) return TaskStatus::complete; + // Copy shared boundary data + // This pack has to have only one block + auto gid = pack.GetGIDHost(0); + parthenon::MeshBlock *pmb; + for (auto &pmbd : md->GetAllBlockData()) + if (pmbd->GetParentPointer()->gid == gid) pmb = pmbd->GetParentPointer(); + + // Now check which row of corner shared elements has been set by an owning block + parthenon::CellCentOffsets offsetX1(2 * (pos % 2) - 1, 0, 0); + parthenon::CellCentOffsets offsetX2(0, 2 * (pos / 2) - 1, 0); + int dir = X2DIR; + for (auto &neighbor : pmb->neighbors) { + if (neighbor.offsets == offsetX1 && neighbor.origin_ownership(offsetX2)) + dir = X1DIR; + if (neighbor.offsets == offsetX2 && neighbor.origin_ownership(offsetX1)) + dir = X2DIR; + } + + // Select the reference location on the node shared by all three blocks const int icorner = (offset(X1DIR) == parthenon::Offset::Low) ? ib_in.s : ib_in.e; const int jcorner = (offset(X2DIR) == parthenon::Offset::Low) ? jb_in.s : jb_in.e; - IndexRange ib = offset(X1DIR) == parthenon::Offset::Low ? IndexRange{0, ib_in.s - 1} : IndexRange{ib_in.e + 1, ib_in.e + parthenon::Globals::nghost}; - IndexRange jb = offset(X2DIR) == parthenon::Offset::Low ? IndexRange{0, jb_in.s - 1} : IndexRange{jb_in.e + 1, ib_in.e + parthenon::Globals::nghost}; + // Select the index space of elements that need to get overwritten + IndexRange ib = offset(X1DIR) == parthenon::Offset::Low ? IndexRange{0, ib_in.s - (dir == X2DIR)} : IndexRange{ib_in.e + (dir == X2DIR), ib_in.e + parthenon::Globals::nghost}; + IndexRange jb = offset(X2DIR) == parthenon::Offset::Low ? IndexRange{0, jb_in.s - (dir == X1DIR)} : IndexRange{jb_in.e + (dir == X1DIR), ib_in.e + parthenon::Globals::nghost}; parthenon::par_for( parthenon::loop_pattern_mdrange_tag, "SetPosition", DevExecSpace(), 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { const int enclosed_area = std::abs((i - icorner) * (j - jcorner)); - int ioff{1}; - while (ioff * ioff < enclosed_area) ioff++; - ioff *= offset(X1DIR); - pack(b, TE::NN, position(0), k, j, i) = pack(b, TE::NN, position(0), k, jcorner, icorner + ioff); // x-position - pack(b, TE::NN, position(1), k, j, i) = pack(b, TE::NN, position(1), k, jcorner, icorner + ioff); // y-position + int ioff{0}; + int joff{0}; + // This is an element shared by both of the neighbors, but only one communicated it, so copy it to the other side of the corner + if (enclosed_area == 0) { + joff = std::abs(icorner - i) * offset(X2DIR); + ioff = std::abs(jcorner - j) * offset(X1DIR); + } + // These are non-existent elements that we fill with data copied from their nearest neighbors + else { + int off = 1; + while (off * off < enclosed_area) off++; + ioff = offset(X1DIR) * off * (dir == X1DIR); + joff = offset(X2DIR) * off * (dir == X2DIR); + } + pack(b, TE::NN, position(0), k, j, i) = pack(b, TE::NN, position(0), k, jcorner + joff, icorner + ioff); // x-position + pack(b, TE::NN, position(1), k, j, i) = pack(b, TE::NN, position(1), k, jcorner + joff, icorner + ioff); // y-position }); } pos++; From 2ce72e3833b0ae56981a8e5a1b71bf037db3c1e0 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Thu, 8 May 2025 15:55:25 -0600 Subject: [PATCH 18/20] format and lint --- .../boundary_exchange/boundary_exchange.cpp | 81 +++++++++++-------- .../boundary_exchange_driver.cpp | 36 +++++---- src/bvals/comms/bnd_info.cpp | 5 +- src/bvals/comms/boundary_communication.cpp | 8 +- .../logical_coordinate_transformation.cpp | 25 +++--- .../logical_coordinate_transformation.hpp | 6 +- src/mesh/mesh-gmg.cpp | 8 +- src/utils/cell_center_offsets.hpp | 6 +- src/utils/indexer.hpp | 10 +-- 9 files changed, 103 insertions(+), 82 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp index d9c09f6aa543..7f98004dfb3a 100644 --- a/example/boundary_exchange/boundary_exchange.cpp +++ b/example/boundary_exchange/boundary_exchange.cpp @@ -103,41 +103,41 @@ TaskStatus SetBlockValues(MeshData *md) { } struct ParameterizedLine { - enum class ltype {straight, arc}; + enum class ltype { straight, arc }; const Real x1, y1; const Real x2, y2; Real r, delta, phi; ltype type; - + using node = parthenon::forest::Node; ParameterizedLine(std::shared_ptr start, std::shared_ptr end) - : x1{start->x[0]}, y1{start->x[1]}, x2{end->x[0]}, y2{end->x[1]}, type{ltype::straight} { + : x1{start->x[0]}, y1{start->x[1]}, x2{end->x[0]}, y2{end->x[1]}, + type{ltype::straight} { Real d1 = std::sqrt(x1 * x1 + y1 * y1); - Real d2 = std::sqrt(x2 * x2 + y2 * y2); - if (std::abs(d1 - d2) < 1.e-8 && d1 > 1.1) { + Real d2 = std::sqrt(x2 * x2 + y2 * y2); + if (std::abs(d1 - d2) < 1.e-8 && d1 > 1.1) { type = ltype::arc; - r = d1; - delta = M_PI / 4.0; + r = d1; + delta = M_PI / 4.0; phi = 0.0; if (y1 > 1.e-5) phi = delta; } } KOKKOS_INLINE_FUNCTION - Real GetX(Real u) const { + Real GetX(Real u) const { if (type == ltype::straight) - return x1 * (1.0 - u) + x2 * u; - else + return x1 * (1.0 - u) + x2 * u; + else return r * cos(delta * u + phi); } - KOKKOS_INLINE_FUNCTION - Real GetY(Real u) const { + Real GetY(Real u) const { if (type == ltype::straight) - return y1 * (1.0 - u) + y2 * u; - else + return y1 * (1.0 - u) + y2 * u; + else return r * sin(delta * u + phi); } }; @@ -207,16 +207,17 @@ TaskStatus FixTrivalentNodes2D(MeshData *md) { (pnode->associated_faces.size() == 3) && !pnode->on_physical_boundary; if (trivalent) { parthenon::CellCentOffsets offset(2 * (pos % 2) - 1, 2 * (pos / 2) - 1, 0); - auto pack = desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId(), offset)); + auto pack = + desc.GetPack(md, parthenon::GetBlockSelector::OnTree(ptree->GetId(), offset)); if (pack.GetNBlocks() == 0) return TaskStatus::complete; - // Copy shared boundary data - - // This pack has to have only one block + // Copy shared boundary data + + // This pack has to have only one block auto gid = pack.GetGIDHost(0); - parthenon::MeshBlock *pmb; + parthenon::MeshBlock *pmb; for (auto &pmbd : md->GetAllBlockData()) - if (pmbd->GetParentPointer()->gid == gid) pmb = pmbd->GetParentPointer(); - + if (pmbd->GetParentPointer()->gid == gid) pmb = pmbd->GetParentPointer(); + // Now check which row of corner shared elements has been set by an owning block parthenon::CellCentOffsets offsetX1(2 * (pos % 2) - 1, 0, 0); parthenon::CellCentOffsets offsetX2(0, 2 * (pos / 2) - 1, 0); @@ -232,29 +233,41 @@ TaskStatus FixTrivalentNodes2D(MeshData *md) { const int icorner = (offset(X1DIR) == parthenon::Offset::Low) ? ib_in.s : ib_in.e; const int jcorner = (offset(X2DIR) == parthenon::Offset::Low) ? jb_in.s : jb_in.e; // Select the index space of elements that need to get overwritten - IndexRange ib = offset(X1DIR) == parthenon::Offset::Low ? IndexRange{0, ib_in.s - (dir == X2DIR)} : IndexRange{ib_in.e + (dir == X2DIR), ib_in.e + parthenon::Globals::nghost}; - IndexRange jb = offset(X2DIR) == parthenon::Offset::Low ? IndexRange{0, jb_in.s - (dir == X1DIR)} : IndexRange{jb_in.e + (dir == X1DIR), ib_in.e + parthenon::Globals::nghost}; + IndexRange ib = offset(X1DIR) == parthenon::Offset::Low + ? IndexRange{0, ib_in.s - (dir == X2DIR)} + : IndexRange{ib_in.e + (dir == X2DIR), + ib_in.e + parthenon::Globals::nghost}; + IndexRange jb = offset(X2DIR) == parthenon::Offset::Low + ? IndexRange{0, jb_in.s - (dir == X1DIR)} + : IndexRange{jb_in.e + (dir == X1DIR), + ib_in.e + parthenon::Globals::nghost}; parthenon::par_for( parthenon::loop_pattern_mdrange_tag, "SetPosition", DevExecSpace(), 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) { - const int enclosed_area = std::abs((i - icorner) * (j - jcorner)); + const int enclosed_area = std::abs((i - icorner) * (j - jcorner)); int ioff{0}; int joff{0}; - // This is an element shared by both of the neighbors, but only one communicated it, so copy it to the other side of the corner - if (enclosed_area == 0) { - joff = std::abs(icorner - i) * offset(X2DIR); + // This is an element shared by both of the neighbors, but only one + // communicated it, so copy it to the other side of the corner + if (enclosed_area == 0) { + joff = std::abs(icorner - i) * offset(X2DIR); ioff = std::abs(jcorner - j) * offset(X1DIR); - } - // These are non-existent elements that we fill with data copied from their nearest neighbors - else { - int off = 1; - while (off * off < enclosed_area) off++; + } else { + // These are non-existent elements that we fill with data copied from + // their nearest neighbors + int off = 1; + while (off * off < enclosed_area) + off++; ioff = offset(X1DIR) * off * (dir == X1DIR); joff = offset(X2DIR) * off * (dir == X2DIR); } - pack(b, TE::NN, position(0), k, j, i) = pack(b, TE::NN, position(0), k, jcorner + joff, icorner + ioff); // x-position - pack(b, TE::NN, position(1), k, j, i) = pack(b, TE::NN, position(1), k, jcorner + joff, icorner + ioff); // y-position + pack(b, TE::NN, position(0), k, j, i) = + pack(b, TE::NN, position(0), k, jcorner + joff, + icorner + ioff); // x-position + pack(b, TE::NN, position(1), k, j, i) = + pack(b, TE::NN, position(1), k, jcorner + joff, + icorner + ioff); // y-position }); } pos++; diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp index 17325c157c56..c523d5fe1f1a 100644 --- a/example/boundary_exchange/boundary_exchange_driver.cpp +++ b/example/boundary_exchange/boundary_exchange_driver.cpp @@ -105,23 +105,25 @@ parthenon::forest::ForestDefinition n_blocks(int nblocks) { for (int t = 0; t < nblocks; ++t) forest_def.AddFace( t, {n[2 * t + 1], n[2 * t], n[(2 * t + 2) % (2 * nblocks)], n[2 * nblocks]}); - + int level = 1; - printf("loc: %i\n", (1 << level) - 1); - forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, level, (1 << level) - 1, (1 << level) - 1, 0)); - forest_def.AddInitialRefinement(parthenon::LogicalLocation(1, level, (1 << level) - 1, (1 << level) - 1, 0)); - forest_def.AddInitialRefinement(parthenon::LogicalLocation(2, level, (1 << level) - 1, (1 << level) - 1, 0)); + forest_def.AddInitialRefinement( + parthenon::LogicalLocation(0, level, (1 << level) - 1, (1 << level) - 1, 0)); + forest_def.AddInitialRefinement( + parthenon::LogicalLocation(1, level, (1 << level) - 1, (1 << level) - 1, 0)); + forest_def.AddInitialRefinement( + parthenon::LogicalLocation(2, level, (1 << level) - 1, (1 << level) - 1, 0)); return forest_def; } -parthenon::forest::ForestDefinition cubed_circle(std::vector radii) { +parthenon::forest::ForestDefinition cubed_circle(std::vector radii) { using namespace parthenon::forest; - std::vector> nodes_x, nodes_y, nodes_c; + std::vector> nodes_x, nodes_y, nodes_c; std::size_t node_idx{0}; std::shared_ptr origin = Node::create(node_idx++, {0.0, 0.0}); double fac = 1.0 / sqrt(2.0); - for (auto rad : radii) { + for (auto rad : radii) { nodes_x.push_back(Node::create(node_idx++, {rad, 0.0})); nodes_c.push_back(Node::create(node_idx++, {rad * fac, rad * fac})); nodes_y.push_back(Node::create(node_idx++, {0.0, rad})); @@ -133,15 +135,17 @@ parthenon::forest::ForestDefinition cubed_circle(std::vector radii) { forest_def.AddBC(Edge({origin, nodes_x[0]})); forest_def.AddBC(Edge({origin, nodes_y[0]})); for (int i = 0; i < nodes_x.size() - 1; ++i) { - forest_def.AddFace(face_idx++, {nodes_x[i], nodes_x[i + 1], nodes_c[i], nodes_c[i + 1]}); - forest_def.AddFace(face_idx++, {nodes_c[i], nodes_c[i + 1], nodes_y[i], nodes_y[i + 1]}); - forest_def.AddBC(Edge({nodes_x[i], nodes_x[i + 1]})); - forest_def.AddBC(Edge({nodes_y[i], nodes_y[i + 1]})); - } + forest_def.AddFace(face_idx++, + {nodes_x[i], nodes_x[i + 1], nodes_c[i], nodes_c[i + 1]}); + forest_def.AddFace(face_idx++, + {nodes_c[i], nodes_c[i + 1], nodes_y[i], nodes_y[i + 1]}); + forest_def.AddBC(Edge({nodes_x[i], nodes_x[i + 1]})); + forest_def.AddBC(Edge({nodes_y[i], nodes_y[i + 1]})); + } forest_def.AddBC(Edge({nodes_x.back(), nodes_c.back()})); forest_def.AddBC(Edge({nodes_c.back(), nodes_y.back()})); - - return forest_def; + + return forest_def; } int main(int argc, char *argv[]) { @@ -163,7 +167,7 @@ int main(int argc, char *argv[]) { } pman.ParthenonInitPackagesAndMesh(cubed_circle({1.0, 2.0, 3.0})); - //pman.ParthenonInitPackagesAndMesh(n_blocks(3)); + // pman.ParthenonInitPackagesAndMesh(n_blocks(3)); // This needs to be scoped so that the driver object is destructed before Finalize { diff --git a/src/bvals/comms/bnd_info.cpp b/src/bvals/comms/bnd_info.cpp index 4467927ae849..370dcac9962b 100644 --- a/src/bvals/comms/bnd_info.cpp +++ b/src/bvals/comms/bnd_info.cpp @@ -246,9 +246,10 @@ CalcIndices(const NeighborBlock &nb, MeshBlock *pmb, } // We are working in the frame of the origin block (i.e. the receiving block), // so we use the ownership array in that frame to get the masking initially. - // This mask must then be transformed to the frame of the neighbor block (i.e. + // This mask must then be transformed to the frame of the neighbor block (i.e. // sending block) since the index range we are masking is defined in that frame. - owns = lcoord_trans.Transform(GetIndexRangeMaskFromOwnership(el, nb.origin_ownership, sox1, sox2, sox3)); + owns = lcoord_trans.Transform( + GetIndexRangeMaskFromOwnership(el, nb.origin_ownership, sox1, sox2, sox3)); } return SpatiallyMaskedIndexer6D(owns, {0, tensor_shape[0] - 1}, {0, tensor_shape[1] - 1}, {0, tensor_shape[2] - 1}, diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp index 8a6d43d8f06a..6d4ca985923e 100644 --- a/src/bvals/comms/boundary_communication.cpp +++ b/src/bvals/comms/boundary_communication.cpp @@ -303,8 +303,8 @@ TaskStatus SetBounds(std::shared_ptr> &md) { Kokkos::parallel_for( Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) { if (idxer.IsActive(kk, jj, ii + m)) { - const auto [il, jl, kl] = - lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); + const auto [il, jl, kl] = lcoord_trans.InverseTransform( + std::array{ii + m, jj, kk}); var(iel, tt, uu, vv, kl, jl, il) = fac * buf[m]; } }); @@ -324,8 +324,8 @@ TaskStatus SetBounds(std::shared_ptr> &md) { Kokkos::parallel_for( Kokkos::ThreadVectorRange<>(team_member, Ni), [&](int m) { if (idxer.IsActive(kk, jj, ii + m)) { - const auto [il, jl, kl] = - lcoord_trans.InverseTransform(std::array{ii + m, jj, kk}); + const auto [il, jl, kl] = lcoord_trans.InverseTransform( + std::array{ii + m, jj, kk}); var(iel, tt, uu, vv, kl, jl, il) = default_val; } }); diff --git a/src/mesh/forest/logical_coordinate_transformation.cpp b/src/mesh/forest/logical_coordinate_transformation.cpp index 7f2f075f0987..55242844c867 100644 --- a/src/mesh/forest/logical_coordinate_transformation.cpp +++ b/src/mesh/forest/logical_coordinate_transformation.cpp @@ -100,7 +100,8 @@ CellCentOffsets LogicalCoordinateTransformation::Transform(CellCentOffsets in) c return out; } -CellCentOffsets LogicalCoordinateTransformation::InverseTransform(CellCentOffsets in) const { +CellCentOffsets +LogicalCoordinateTransformation::InverseTransform(CellCentOffsets in) const { CellCentOffsets out; for (int dir = 0; dir < 3; ++dir) { const int indir = abs(dir_connection[dir]); @@ -109,11 +110,12 @@ CellCentOffsets LogicalCoordinateTransformation::InverseTransform(CellCentOffset return out; } -block_ownership_t LogicalCoordinateTransformation::Transform(const block_ownership_t &in) const { - block_ownership_t out; - for (int ox3 : {-1 , 0, 1}) { - for (int ox2 : {-1 , 0, 1}) { - for (int ox1 : {-1 , 0, 1}) { +block_ownership_t +LogicalCoordinateTransformation::Transform(const block_ownership_t &in) const { + block_ownership_t out; + for (int ox3 : {-1, 0, 1}) { + for (int ox2 : {-1, 0, 1}) { + for (int ox1 : {-1, 0, 1}) { auto offset = CellCentOffsets(ox1, ox2, ox3); out(Transform(offset)) = in(offset); } @@ -122,11 +124,12 @@ block_ownership_t LogicalCoordinateTransformation::Transform(const block_ownersh return out; } -block_ownership_t LogicalCoordinateTransformation::InverseTransform(const block_ownership_t &in) const { - block_ownership_t out; - for (int ox3 : {-1 , 0, 1}) { - for (int ox2 : {-1 , 0, 1}) { - for (int ox1 : {-1 , 0, 1}) { +block_ownership_t +LogicalCoordinateTransformation::InverseTransform(const block_ownership_t &in) const { + block_ownership_t out; + for (int ox3 : {-1, 0, 1}) { + for (int ox2 : {-1, 0, 1}) { + for (int ox1 : {-1, 0, 1}) { auto offset = CellCentOffsets(ox1, ox2, ox3); out(offset) = in(Transform(offset)); } diff --git a/src/mesh/forest/logical_coordinate_transformation.hpp b/src/mesh/forest/logical_coordinate_transformation.hpp index 8faf16ab1dcd..57aec7db9691 100644 --- a/src/mesh/forest/logical_coordinate_transformation.hpp +++ b/src/mesh/forest/logical_coordinate_transformation.hpp @@ -117,10 +117,10 @@ struct NeighborLocation { : global_loc(g), origin_loc(o), lcoord_trans(lcoord_trans) {} // Global location of neighboring block LogicalLocation global_loc; - // Logical location of neighboring block in index space of origin block + // Logical location of neighboring block in index space of origin block LogicalLocation origin_loc; - // Coordinate transform to coords of neighbor tree from origin tree - // (i.e. global_loc = lcoord_trans.Transform(origin_loc)) + // Coordinate transform to coords of neighbor tree from origin tree + // (i.e. global_loc = lcoord_trans.Transform(origin_loc)) LogicalCoordinateTransformation lcoord_trans; }; diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp index 965fa9ca61dd..0910050af08d 100644 --- a/src/mesh/mesh-gmg.cpp +++ b/src/mesh/mesh-gmg.cpp @@ -76,14 +76,16 @@ void Mesh::SetMeshBlockNeighbors( all_neighbors.emplace_back(pmb->pmy_mesh, nloc.global_loc, nloc.origin_loc, ranklist[lgid], gid, offsets, bid, tid, f[0], f[1]); - // Set neighbor block ownership (note that this is ownership in coordinates of the neighbor block) + // Set neighbor block ownership (note that this is ownership in coordinates of the + // neighbor block) auto &nb = all_neighbors.back(); auto neighbor_neighbors = forest.FindNeighbors(nloc.global_loc, grid_id); - nb.ownership = DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined); + nb.ownership = + DetermineOwnership(nloc.global_loc, neighbor_neighbors, newly_refined); nb.ownership.initialized = true; nb.origin_ownership = nloc.lcoord_trans.InverseTransform(nb.ownership); nb.origin_ownership.initialized = true; - + // Set logical coordinate transformation from this block to the neighbor nb.lcoord_trans = nloc.lcoord_trans; } diff --git a/src/utils/cell_center_offsets.hpp b/src/utils/cell_center_offsets.hpp index b41b625d85b3..86cbc01fc3a0 100644 --- a/src/utils/cell_center_offsets.hpp +++ b/src/utils/cell_center_offsets.hpp @@ -73,10 +73,8 @@ struct CellCentOffsets { // (in cyclic order, XY, YZ, ZX, XYZ) along with the offset of the // element in that direction from the cell center. std::vector> GetNormals() const; - - bool operator==(const CellCentOffsets& other) const { - return u == other.u; - } + + bool operator==(const CellCentOffsets &other) const { return u == other.u; } bool IsNode() const { return 3 == abs(static_cast(u[0])) + abs(static_cast(u[1])) + diff --git a/src/utils/indexer.hpp b/src/utils/indexer.hpp index a452c0beafe4..9ffb384b9d62 100644 --- a/src/utils/indexer.hpp +++ b/src/utils/indexer.hpp @@ -35,13 +35,13 @@ struct block_ownership_t { bool &operator()(int ox1, int ox2, int ox3) { return ownership[ox1 + 1][ox2 + 1][ox3 + 1]; } - KOKKOS_FORCEINLINE_FUNCTION - bool &operator()(const CellCentOffsets &offset) { - return ownership[offset[0]+1][offset[1]+1][offset[2]+1]; + KOKKOS_FORCEINLINE_FUNCTION + bool &operator()(const CellCentOffsets &offset) { + return ownership[offset[0] + 1][offset[1] + 1][offset[2] + 1]; } KOKKOS_FORCEINLINE_FUNCTION - const bool &operator()(const CellCentOffsets &offset) const { - return ownership[offset[0]+1][offset[1]+1][offset[2]+1]; + const bool &operator()(const CellCentOffsets &offset) const { + return ownership[offset[0] + 1][offset[1] + 1][offset[2] + 1]; } KOKKOS_FORCEINLINE_FUNCTION From 9d3b6c047100b9cbf826eb115db18d29dfd0849f Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Fri, 9 May 2025 19:32:11 -0600 Subject: [PATCH 19/20] split out prolongation --- src/bvals/comms/boundary_communication.cpp | 90 +++++++++++++++++++--- src/bvals/comms/bvals_in_one.hpp | 5 +- 2 files changed, 84 insertions(+), 11 deletions(-) diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp index 6d4ca985923e..222c1f5c58fc 100644 --- a/src/bvals/comms/boundary_communication.cpp +++ b/src/bvals/comms/boundary_communication.cpp @@ -403,10 +403,77 @@ ProlongateBounds(std::shared_ptr> &); template TaskStatus ProlongateBounds(std::shared_ptr> &); +template +TaskStatus ProlongateInternalBounds(std::shared_ptr> &md) { + PARTHENON_INSTRUMENT + + Mesh *pmesh = md->GetMeshPointer(); + auto &cache = md->GetBvarsCache().GetSubCache(bound_type, false); + + auto [rebuild, nbound] = CheckReceiveBufferCacheForRebuild(md); + + if (rebuild) { + if constexpr (bound_type == BoundaryType::gmg_prolongate_recv) { + RebuildBufferCache(md, nbound, BndInfo::GetSetBndInfo, + ProResInfo::GetInteriorProlongate); + } else if constexpr (bound_type == BoundaryType::gmg_restrict_recv) { + RebuildBufferCache(md, nbound, BndInfo::GetSetBndInfo, + ProResInfo::GetNull); + } else { + RebuildBufferCache(md, nbound, BndInfo::GetSetBndInfo, + ProResInfo::GetSet); + } + } + + if (nbound > 0 && pmesh->multilevel && md->NumBlocks() > 0) { + auto pmb = md->GetBlockData(0)->GetBlockPointer(); + StateDescriptor *resolved_packages = pmb->resolved_packages.get(); + + // Prolongate from coarse buffer + refinement::ProlongateInternal(resolved_packages, cache.prores_cache, pmb->cellbounds, + pmb->c_cellbounds); + } + return TaskStatus::complete; +} + +template +TaskStatus ProlongateSharedBounds(std::shared_ptr> &md) { + PARTHENON_INSTRUMENT + + Mesh *pmesh = md->GetMeshPointer(); + auto &cache = md->GetBvarsCache().GetSubCache(bound_type, false); + + auto [rebuild, nbound] = CheckReceiveBufferCacheForRebuild(md); + + if (rebuild) { + if constexpr (bound_type == BoundaryType::gmg_prolongate_recv) { + RebuildBufferCache(md, nbound, BndInfo::GetSetBndInfo, + ProResInfo::GetInteriorProlongate); + } else if constexpr (bound_type == BoundaryType::gmg_restrict_recv) { + RebuildBufferCache(md, nbound, BndInfo::GetSetBndInfo, + ProResInfo::GetNull); + } else { + RebuildBufferCache(md, nbound, BndInfo::GetSetBndInfo, + ProResInfo::GetSet); + } + } + + if (nbound > 0 && pmesh->multilevel && md->NumBlocks() > 0) { + auto pmb = md->GetBlockData(0)->GetBlockPointer(); + StateDescriptor *resolved_packages = pmb->resolved_packages.get(); + + // Prolongate from coarse buffer + refinement::ProlongateShared(resolved_packages, cache.prores_cache, pmb->cellbounds, + pmb->c_cellbounds); + } + return TaskStatus::complete; +} + // Adds all relevant boundary communication to a single task list template TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl, - std::shared_ptr> &md, bool multilevel) { + std::shared_ptr> &md, bool multilevel, + bound_task_adder_t extra_bounds) { // TODO(LFR): Splitting up the boundary tasks while doing prolongation can cause some // possible issues for sparse fields. In particular, the order in which // fields are allocated and then set could potentially result in different @@ -439,19 +506,22 @@ TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl, auto pro = set; if (md->GetMeshPointer()->multilevel) { auto cbound = tl.AddTask(set, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, true); + if (extra_bounds) cbound = extra_bounds(cbound, &tl, md, true); pro = tl.AddTask(cbound, TF(ProlongateBounds), md); + auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false); + if (extra_bounds) fbound = extra_bounds(fbound, &tl, md, false); + return tl.AddTask(fbound, TF(ProlongateInternalBounds), md); + } else { + auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false); + if (extra_bounds) fbound = extra_bounds(fbound, &tl, md, false); + return fbound; } - auto fbound = tl.AddTask(pro, TF(ApplyBoundaryConditionsOnCoarseOrFineMD), md, false); - - return fbound; } -template TaskID -AddBoundaryExchangeTasks(TaskID, TaskList &, - std::shared_ptr> &, bool); +template TaskID AddBoundaryExchangeTasks( + TaskID, TaskList &, std::shared_ptr> &, bool, bound_task_adder_t); -template TaskID -AddBoundaryExchangeTasks(TaskID, TaskList &, - std::shared_ptr> &, bool); +template TaskID AddBoundaryExchangeTasks( + TaskID, TaskList &, std::shared_ptr> &, bool, bound_task_adder_t); TaskID AddFluxCorrectionTasks(TaskID dependency, TaskList &tl, std::shared_ptr> &md, bool multilevel) { diff --git a/src/bvals/comms/bvals_in_one.hpp b/src/bvals/comms/bvals_in_one.hpp index 0180bd6a8305..a38952ccc58e 100644 --- a/src/bvals/comms/bvals_in_one.hpp +++ b/src/bvals/comms/bvals_in_one.hpp @@ -80,9 +80,12 @@ static TaskStatus SetFluxCorrections(std::shared_ptr> &md) { } // Adds all relevant boundary communication to a single task list +using bound_task_adder_t = + std::function>, bool)>; template TaskID AddBoundaryExchangeTasks(TaskID dependency, TaskList &tl, - std::shared_ptr> &md, bool multilevel); + std::shared_ptr> &md, bool multilevel, + bound_task_adder_t extra_bounds = bound_task_adder_t()); // Adds all relevant flux correction tasks to a single task list TaskID AddFluxCorrectionTasks(TaskID dependency, TaskList &tl, From b24eaba361f69668e8ad5c121b85abd8779fd7c7 Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Fri, 9 May 2025 19:32:47 -0600 Subject: [PATCH 20/20] make trivalent fixup work with AMR --- example/boundary_exchange/boundary_exchange.cpp | 16 +++++++++++----- example/boundary_exchange/boundary_exchange.hpp | 2 +- .../boundary_exchange_driver.cpp | 13 +++++++++++-- 3 files changed, 23 insertions(+), 8 deletions(-) diff --git a/example/boundary_exchange/boundary_exchange.cpp b/example/boundary_exchange/boundary_exchange.cpp index 7f98004dfb3a..6fdc42f572a1 100644 --- a/example/boundary_exchange/boundary_exchange.cpp +++ b/example/boundary_exchange/boundary_exchange.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -190,14 +191,19 @@ TaskStatus SetCoordinates(MeshData *md) { return TaskStatus::complete; } -TaskStatus FixTrivalentNodes2D(MeshData *md) { +TaskStatus FixTrivalentNodes2D(MeshData *md, bool coarse) { using TE = parthenon::TopologicalElement; auto pmesh = md->GetMeshPointer(); - auto desc = parthenon::MakePackDescriptor(md); + std::set opts{}; + if (coarse) opts.insert(parthenon::PDOpt::Coarse); + + auto desc = parthenon::MakePackDescriptor(md, {}, opts); - IndexRange ib_in = md->GetBoundsI(IndexDomain::interior, TE::NN); - IndexRange jb_in = md->GetBoundsJ(IndexDomain::interior, TE::NN); - IndexRange kb = md->GetBoundsK(IndexDomain::interior, TE::NN); + const parthenon::CellLevel clevel = + coarse ? parthenon::CellLevel::coarse : parthenon::CellLevel::same; + IndexRange ib_in = md->GetBoundsI(clevel, IndexDomain::interior, TE::NN); + IndexRange jb_in = md->GetBoundsJ(clevel, IndexDomain::interior, TE::NN); + IndexRange kb = md->GetBoundsK(clevel, IndexDomain::interior, TE::NN); for (auto &ptree : pmesh->forest.GetTrees()) { auto tree_id = ptree->GetId(); diff --git a/example/boundary_exchange/boundary_exchange.hpp b/example/boundary_exchange/boundary_exchange.hpp index f8853e00be63..bd8a60cbfcdc 100644 --- a/example/boundary_exchange/boundary_exchange.hpp +++ b/example/boundary_exchange/boundary_exchange.hpp @@ -43,7 +43,7 @@ struct position : public parthenon::variable_names::base_t { TaskStatus SetBlockValues(MeshData *rc); TaskStatus SetCoordinates(MeshData *rc); -TaskStatus FixTrivalentNodes2D(MeshData *rc); +TaskStatus FixTrivalentNodes2D(MeshData *rc, bool coarse); std::shared_ptr Initialize(ParameterInput *pin); } // namespace boundary_exchange diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp index c523d5fe1f1a..0b0fc408a7d4 100644 --- a/example/boundary_exchange/boundary_exchange_driver.cpp +++ b/example/boundary_exchange/boundary_exchange_driver.cpp @@ -145,6 +145,10 @@ parthenon::forest::ForestDefinition cubed_circle(std::vector radii) { forest_def.AddBC(Edge({nodes_x.back(), nodes_c.back()})); forest_def.AddBC(Edge({nodes_c.back(), nodes_y.back()})); + int level = 2; + forest_def.AddInitialRefinement( + parthenon::LogicalLocation(1, level, 0, (1 << level) - 1, 0)); + return forest_def; } @@ -223,8 +227,13 @@ TaskCollection BoundaryExchangeDriver::MakeTaskCollection(T &blocks) { TaskID none(0); auto fill = tl.AddTask(none, SetBlockValues, md.get()); auto set_coords = tl.AddTask(none, SetCoordinates, md.get()); - auto bound = AddBoundaryExchangeTasks(fill | set_coords, tl, md, true); - auto fix = tl.AddTask(bound, FixTrivalentNodes2D, md.get()); + auto bound = AddBoundaryExchangeTasks( + fill | set_coords, tl, md, true, + [](TaskID depend, TaskList *ptl, std::shared_ptr> md, + bool coarse) { + return ptl->AddTask(depend, FixTrivalentNodes2D, md.get(), coarse); + }); + // auto fix = tl.AddTask(bound, FixTrivalentNodes2D, md.get(), false); } }