From 40dd59647dee23a7a950f9ce9fd6b287231b467e Mon Sep 17 00:00:00 2001 From: Luke Roberts Date: Tue, 20 Jun 2023 15:57:45 -0600 Subject: [PATCH] Add constant derivative boundary condition and expose `GenericBC` (#884) * Change CellV* to V* * Make boundary communication fail if non-cell centered variables have FillGhost or WithFluxes * Remove cc * Rename namespace cell_centered_bvars to var_boundary_comm and fix error * Remove cc from bvals_cc_in_one * Format and lint * Update changelog * Add multipointer * Allow for arbitrary changes of ParArrayND rank and move all array to ParArrayND * Fix particle indexing * Simple sparse pack based interface for face, edge and node variables * Format and lint * Update copyright * update year * Allow for variables w/o topology metadata in packs * Allow ParArrayGeneric::Get to take pairs as args * Remove extraneous ParArray types * Forgot to add new header to CMake * Actually add array to tuple * cleanup, format, and lint * Allow for packing other topological types * Add small test of packing face variables * format and lint * Update changelog * try to fix HIP error and small change * Act on brryan's comments * Small * Format and lint * Update domain functions to deal with varying topological elements * Add Indexer class template * Use Indexer for buffer packing * Start rewriting index calculations * Remove GetTensor stuff * Comment on chosen numbering * Add some comments * Update sparse pack interface to take topological elements * Switch to multiple indexers and format * Remove var_boundary_com namespace * Start switching to new indexing * working for sparse advection * Replace buffer index calculation routines and generalize * In theory should communication should work now for all field types * Avoid duplication * Actually set the number of topological elements for a buffer * Update buffer size calculation * Small * Fix out of range error * Change internal storage so that we don't trigger HtoD * Fix bug in assertion * Switch to std::array since that seems to give no complaints on device * Remove explicit indices from BndInfo * Remove defaulted functions * Switch to enum for choosing boundaries * Add template option for ghost restriction range * Switch to separate indexers for prolongation and restriction * higher precision error output * Stop using RestrictGhostHalos * Remove RestrictGhostHalos * Remove BoundaryType::restricted * Switch to prolongation in one * Format and lint * Remove old boundary prolongation machinery * A temporary fix for MPI runs * format and lint * Use indexers directly * Fix device compilation errors * Format and lint * fix typo * Add generalized volume function * Add function for recovering generic central positions * Generalize restriction * Still working... * Still passing... * Generalize shared restriction * Generalize prolongation and restriction registration to work for all topological types * Index indexers correctly * add element index * Base routine for doing general internal prolongation * Move submanifold routine and add routines to check if operations are required * Add OperationRequired checks to tests * Add all possible cell element types * Remove error causing loop * Only select submanifolds, not the current manifold * Add internal prolongation calls * format and lint * format and lint... * Add internal prolongation, part 2 * fix debug issue * Add some clarifying comments * Switch to using loop indexers required for internal prolongation * Add some comments and a diagram * small * Fix indexing bug when switching to new indexing * Add debug requirement suggested by brryan * fix bug * Move metadata constructor to cpp file * Rename topological elements * format and lint * Rename and catch some bugs * format and lint * Only update cell variables * Remove cell boundary communication checks * fix sparse pack bug for face and edge fields * format and lint * format and lint * Fix indexing bugs * Fill boundary arrays correctly for face and edge variables * Remove weird merge duplication * format and lint * Generalize boundary conditions to non-cell centered fields * Revert sparse pack with fluxes behavior * Format and lint * Rename refinement ops * update changelog * Initialize size_ to zero * respond to comments * Correctly fill vector components * Remove unused dims_ * Switch to sparse packs for GenericBoundary conditions * Add constant derivative BC and move GenericBCs to a new header * Format and lint * update changelog * Add fixed value boundary condition * MeshData boundary conditions * Update to new spase pack interface --- CHANGELOG.md | 1 + src/CMakeLists.txt | 1 + src/bvals/boundary_conditions.cpp | 117 ++++--------------- src/bvals/boundary_conditions.hpp | 5 + src/bvals/boundary_conditions_generic.hpp | 130 ++++++++++++++++++++++ src/interface/sparse_pack_base.cpp | 14 ++- src/interface/sparse_pack_base.hpp | 1 - 7 files changed, 169 insertions(+), 100 deletions(-) create mode 100644 src/bvals/boundary_conditions_generic.hpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 5327a4f374b2..3aa52e36dfbd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR 884]](https://github.com/parthenon-hpc-lab/parthenon/pull/884) Add constant derivative BC and expose GenericBC - [[PR 892]](https://github.com/parthenon-hpc-lab/parthenon/pull/892) Cost-based load balancing and memory diagnostics - [[PR 889]](https://github.com/parthenon-hpc-lab/parthenon/pull/889) Add PreCommFillDerived - [[PR 872]](https://github.com/parthenon-hpc-lab/parthenon/pull/872) Boundary communication for non-cell centered fields diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0d02aa6957e1..90e4240538ea 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -105,6 +105,7 @@ add_library(parthenon bvals/comms/tag_map.cpp bvals/comms/tag_map.hpp + bvals/boundary_conditions_generic.hpp bvals/boundary_conditions.cpp bvals/boundary_conditions.hpp diff --git a/src/bvals/boundary_conditions.cpp b/src/bvals/boundary_conditions.cpp index 391649b00aa8..52625bf3c008 100644 --- a/src/bvals/boundary_conditions.cpp +++ b/src/bvals/boundary_conditions.cpp @@ -15,6 +15,7 @@ #include #include "bvals/boundary_conditions.hpp" +#include "bvals/boundary_conditions_generic.hpp" #include "bvals/bvals_interfaces.hpp" #include "defs.hpp" #include "interface/meshblock_data.hpp" @@ -52,136 +53,66 @@ TaskStatus ApplyBoundaryConditionsOnCoarseOrFine(std::shared_ptr -void GenericBC(std::shared_ptr> &rc, bool coarse, - TopologicalElement el = TopologicalElement::CC) { - // make sure DIR is X[123]DIR so we don't have to check again - static_assert(DIR == X1DIR || DIR == X2DIR || DIR == X3DIR, "DIR must be X[123]DIR"); - - // convenient shorthands - constexpr bool X1 = (DIR == X1DIR); - constexpr bool X2 = (DIR == X2DIR); - constexpr bool X3 = (DIR == X3DIR); - constexpr bool INNER = (SIDE == BCSide::Inner); +TaskStatus ApplyBoundaryConditionsMD(std::shared_ptr> &pmd) { + for (int b = 0; b < pmd->NumBlocks(); ++b) + ApplyBoundaryConditions(pmd->GetBlockData(b)); + return TaskStatus::complete; +} - std::shared_ptr pmb = rc->GetBlockPointer(); - const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds; - - const auto &range = X1 ? bounds.GetBoundsI(IndexDomain::interior, el) - : (X2 ? bounds.GetBoundsJ(IndexDomain::interior, el) - : bounds.GetBoundsK(IndexDomain::interior, el)); - const int ref = INNER ? range.s : range.e; - - std::vector flags{Metadata::FillGhost}; - if (GetTopologicalType(el) == TopologicalType::Cell) flags.push_back(Metadata::Cell); - if (GetTopologicalType(el) == TopologicalType::Face) flags.push_back(Metadata::Face); - if (GetTopologicalType(el) == TopologicalType::Edge) flags.push_back(Metadata::Edge); - if (GetTopologicalType(el) == TopologicalType::Node) flags.push_back(Metadata::Node); - - auto q = rc->PackVariables(flags, coarse); - auto nb = IndexRange{0, q.GetDim(4) - 1}; - - std::string label = (TYPE == BCType::Reflect ? "Reflect" : "Outflow"); - label += (INNER ? "Inner" : "Outer"); - label += "X" + std::to_string(DIR); - - constexpr IndexDomain domain = - INNER ? (X1 ? IndexDomain::inner_x1 - : (X2 ? IndexDomain::inner_x2 : IndexDomain::inner_x3)) - : (X1 ? IndexDomain::outer_x1 - : (X2 ? IndexDomain::outer_x2 : IndexDomain::outer_x3)); - - // used for reflections - const int offset = 2 * ref + (INNER ? -1 : 1); - - pmb->par_for_bndry( - label, nb, domain, el, coarse, - KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { - if (!q.IsAllocated(l)) return; - if (TYPE == BCType::Reflect) { - const bool reflect = (q.VectorComponent(l) == DIR); - q(el, l, k, j, i) = - (reflect ? -1.0 : 1.0) * - q(el, l, X3 ? offset - k : k, X2 ? offset - j : j, X1 ? offset - i : i); - } else { - q(el, l, k, j, i) = q(el, l, X3 ? ref : k, X2 ? ref : j, X1 ? ref : i); - } - }); +inline TaskStatus +ApplyBoundaryConditionsOnCoarseOrFineMD(std::shared_ptr> &pmd, + bool coarse) { + for (int b = 0; b < pmd->NumBlocks(); ++b) + ApplyBoundaryConditionsOnCoarseOrFine(pmd->GetBlockData(b), coarse); + return TaskStatus::complete; } void OutflowInnerX1(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void OutflowOuterX1(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void OutflowInnerX2(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void OutflowOuterX2(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void OutflowInnerX3(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void OutflowOuterX3(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void ReflectInnerX1(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void ReflectOuterX1(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void ReflectInnerX2(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void ReflectOuterX2(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void ReflectInnerX3(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } void ReflectOuterX3(std::shared_ptr> &rc, bool coarse) { - using TE = TopologicalElement; - for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) - GenericBC(rc, coarse, el); + GenericBC(rc, coarse); } } // namespace BoundaryFunction diff --git a/src/bvals/boundary_conditions.hpp b/src/bvals/boundary_conditions.hpp index 5c00ecf19f74..3776bac7f321 100644 --- a/src/bvals/boundary_conditions.hpp +++ b/src/bvals/boundary_conditions.hpp @@ -38,6 +38,11 @@ inline TaskStatus ApplyBoundaryConditions(std::shared_ptr> & return ApplyBoundaryConditionsOnCoarseOrFine(rc, false); } +TaskStatus ApplyBoundaryConditionsMD(std::shared_ptr> &pmd); + +TaskStatus ApplyBoundaryConditionsOnCoarseOrFineMD(std::shared_ptr> &pmd, + bool coarse); + namespace BoundaryFunction { void OutflowInnerX1(std::shared_ptr> &rc, bool coarse); diff --git a/src/bvals/boundary_conditions_generic.hpp b/src/bvals/boundary_conditions_generic.hpp new file mode 100644 index 000000000000..e361a932f2a3 --- /dev/null +++ b/src/bvals/boundary_conditions_generic.hpp @@ -0,0 +1,130 @@ +//======================================================================================== +// (C) (or copyright) 2020-2022. 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 BVALS_BOUNDARY_CONDITIONS_GENERIC_HPP_ +#define BVALS_BOUNDARY_CONDITIONS_GENERIC_HPP_ + +#include +#include +#include +#include +#include + +#include "basic_types.hpp" +#include "interface/meshblock_data.hpp" +#include "interface/sparse_pack.hpp" +#include "mesh/domain.hpp" +#include "mesh/mesh.hpp" +#include "mesh/meshblock.hpp" + +namespace parthenon { +namespace BoundaryFunction { + +enum class BCSide { Inner, Outer }; +enum class BCType { Outflow, Reflect, ConstantDeriv, Fixed }; + +template +void GenericBC(std::shared_ptr> &rc, bool coarse, + TopologicalElement el, Real val) { + // make sure DIR is X[123]DIR so we don't have to check again + static_assert(DIR == X1DIR || DIR == X2DIR || DIR == X3DIR, "DIR must be X[123]DIR"); + + // convenient shorthands + constexpr bool X1 = (DIR == X1DIR); + constexpr bool X2 = (DIR == X2DIR); + constexpr bool X3 = (DIR == X3DIR); + constexpr bool INNER = (SIDE == BCSide::Inner); + + std::vector flags{Metadata::FillGhost}; + if (GetTopologicalType(el) == TopologicalType::Cell) flags.push_back(Metadata::Cell); + if (GetTopologicalType(el) == TopologicalType::Face) flags.push_back(Metadata::Face); + if (GetTopologicalType(el) == TopologicalType::Edge) flags.push_back(Metadata::Edge); + if (GetTopologicalType(el) == TopologicalType::Node) flags.push_back(Metadata::Node); + + std::set opts; + if (coarse) opts = {PDOpt::Coarse}; + auto desc = MakePackDescriptor( + rc->GetBlockPointer()->pmy_mesh->resolved_packages.get(), flags, opts); + auto q = desc.GetPack(rc.get()); + const int b = 0; + const int lstart = q.GetLowerBoundHost(b); + const int lend = q.GetUpperBoundHost(b); + if (lend < lstart) return; + auto nb = IndexRange{lstart, lend}; + + std::shared_ptr pmb = rc->GetBlockPointer(); + const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds; + + const auto &range = X1 ? bounds.GetBoundsI(IndexDomain::interior, el) + : (X2 ? bounds.GetBoundsJ(IndexDomain::interior, el) + : bounds.GetBoundsK(IndexDomain::interior, el)); + const int ref = INNER ? range.s : range.e; + + std::string label = (TYPE == BCType::Reflect ? "Reflect" : "Outflow"); + label += (INNER ? "Inner" : "Outer"); + label += "X" + std::to_string(DIR); + + constexpr IndexDomain domain = + INNER ? (X1 ? IndexDomain::inner_x1 + : (X2 ? IndexDomain::inner_x2 : IndexDomain::inner_x3)) + : (X1 ? IndexDomain::outer_x1 + : (X2 ? IndexDomain::outer_x2 : IndexDomain::outer_x3)); + + // used for reflections + const int offset = 2 * ref + (INNER ? -1 : 1); + + // used for derivatives + const int offsetin = INNER; + const int offsetout = !INNER; + pmb->par_for_bndry( + label, nb, domain, el, coarse, + KOKKOS_LAMBDA(const int &l, const int &k, const int &j, const int &i) { + if (TYPE == BCType::Reflect) { + const bool reflect = (q(b, el, l).vector_component == DIR); + q(b, el, l, k, j, i) = + (reflect ? -1.0 : 1.0) * + q(b, el, l, X3 ? offset - k : k, X2 ? offset - j : j, X1 ? offset - i : i); + } else if (TYPE == BCType::ConstantDeriv) { + Real dq = q(b, el, l, X3 ? ref + offsetin : k, X2 ? ref + offsetin : j, + X1 ? ref + offsetin : i) - + q(b, el, l, X3 ? ref - offsetout : k, X2 ? ref - offsetout : j, + X1 ? ref - offsetout : i); + Real delta = 0.0; + if (X1) { + delta = i - ref; + } else if (X2) { + delta = j - ref; + } else { + delta = k - ref; + } + q(b, el, l, k, j, i) = + q(b, el, l, X3 ? ref : k, X2 ? ref : j, X1 ? ref : i) + delta * dq; + } else if (TYPE == BCType::Fixed) { + q(b, el, l, k, j, i) = val; + } else { + q(b, el, l, k, j, i) = q(b, el, l, X3 ? ref : k, X2 ? ref : j, X1 ? ref : i); + } + }); +} + +template +void GenericBC(std::shared_ptr> &rc, bool coarse, Real val = 0.0) { + using TE = TopologicalElement; + for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN}) + GenericBC(rc, coarse, el, val); +} + +} // namespace BoundaryFunction +} // namespace parthenon + +#endif // BVALS_BOUNDARY_CONDITIONS_GENERIC_HPP_ diff --git a/src/interface/sparse_pack_base.cpp b/src/interface/sparse_pack_base.cpp index cd5bc5fd36ec..aa8384e7635c 100644 --- a/src/interface/sparse_pack_base.cpp +++ b/src/interface/sparse_pack_base.cpp @@ -228,6 +228,11 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc) { pack_h(1, b, idx) = pv->data.Get(1, t, u, v); pack_h(2, b, idx) = pv->data.Get(2, t, u, v); } + if (pv->IsSet(Metadata::Vector)) { + pack_h(0, b, idx).vector_component = X1DIR; + pack_h(1, b, idx).vector_component = X2DIR; + pack_h(2, b, idx).vector_component = X3DIR; + } } else { // This is a cell, node, or a variable that doesn't have // topology information if (pack.coarse_) { @@ -235,6 +240,9 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc) { } else { pack_h(0, b, idx) = pv->data.Get(0, t, u, v); } + if (pv->IsSet(Metadata::Vector)) + pack_h(0, b, idx).vector_component = v + 1; + if (desc.with_fluxes && pv->IsSet(Metadata::WithFluxes)) { pack_h(1, b, idx) = pv->flux[X1DIR].Get(0, t, u, v); pack_h(2, b, idx) = pv->flux[X2DIR].Get(0, t, u, v); @@ -266,15 +274,9 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc) { // Record the maximum for easy access pack.bounds_h_(1, block, nvar) = idx - 1; }); - Kokkos::deep_copy(pack.pack_, pack_h); Kokkos::deep_copy(pack.bounds_, pack.bounds_h_); Kokkos::deep_copy(pack.coords_, coords_h); - pack.dims_[1] = pack.nblocks_; - pack.dims_[2] = -1; // Not allowed to ask for the ragged dimension anyway - pack.dims_[3] = pack_h(0, 0, 0).extent_int(0); - pack.dims_[4] = pack_h(0, 0, 0).extent_int(2); - pack.dims_[5] = pack_h(0, 0, 0).extent_int(3); return pack; } diff --git a/src/interface/sparse_pack_base.hpp b/src/interface/sparse_pack_base.hpp index fb5869f2c11b..c17882bd7ba9 100644 --- a/src/interface/sparse_pack_base.hpp +++ b/src/interface/sparse_pack_base.hpp @@ -117,7 +117,6 @@ class SparsePackBase { bool coarse_; bool flat_; int nblocks_; - int dims_[6]; int nvar_; int size_; };