Add constant derivative boundary condition and expose GenericBC

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

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
 

src/bvals/boundary_conditions.cpp

@@ -15,6 +15,7 @@
 #include <vector>
 
 #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<MeshBlockData<R
 
 namespace BoundaryFunction {
 
-enum class BCSide { Inner, Outer };
-enum class BCType { Outflow, Reflect };
-
-template <CoordinateDirection DIR, BCSide SIDE, BCType TYPE>
-void GenericBC(std::shared_ptr<MeshBlockData<Real>> &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<MeshData<Real>> &pmd) {
+  for (int b = 0; b < pmd->NumBlocks(); ++b)
+    ApplyBoundaryConditions(pmd->GetBlockData(b));
+  return TaskStatus::complete;
+}
 
-  std::shared_ptr<MeshBlock> 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<MetadataFlag> 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<MeshData<Real>> &pmd,
+                                        bool coarse) {
+  for (int b = 0; b < pmd->NumBlocks(); ++b)
+    ApplyBoundaryConditionsOnCoarseOrFine(pmd->GetBlockData(b), coarse);
+  return TaskStatus::complete;
 }
 
 void OutflowInnerX1(std::shared_ptr<MeshBlockData<Real>> &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<X1DIR, BCSide::Inner, BCType::Outflow>(rc, coarse, el);
+  GenericBC<X1DIR, BCSide::Inner, BCType::Outflow, variable_names::any>(rc, coarse);
 }
 
 void OutflowOuterX1(std::shared_ptr<MeshBlockData<Real>> &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<X1DIR, BCSide::Outer, BCType::Outflow>(rc, coarse, el);
+  GenericBC<X1DIR, BCSide::Outer, BCType::Outflow, variable_names::any>(rc, coarse);
 }
 
 void OutflowInnerX2(std::shared_ptr<MeshBlockData<Real>> &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<X2DIR, BCSide::Inner, BCType::Outflow>(rc, coarse, el);
+  GenericBC<X2DIR, BCSide::Inner, BCType::Outflow, variable_names::any>(rc, coarse);
 }
 
 void OutflowOuterX2(std::shared_ptr<MeshBlockData<Real>> &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<X2DIR, BCSide::Outer, BCType::Outflow>(rc, coarse, el);
+  GenericBC<X2DIR, BCSide::Outer, BCType::Outflow, variable_names::any>(rc, coarse);
 }
 
 void OutflowInnerX3(std::shared_ptr<MeshBlockData<Real>> &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<X3DIR, BCSide::Inner, BCType::Outflow>(rc, coarse, el);
+  GenericBC<X3DIR, BCSide::Inner, BCType::Outflow, variable_names::any>(rc, coarse);
 }
 
 void OutflowOuterX3(std::shared_ptr<MeshBlockData<Real>> &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<X3DIR, BCSide::Outer, BCType::Outflow>(rc, coarse, el);
+  GenericBC<X3DIR, BCSide::Outer, BCType::Outflow, variable_names::any>(rc, coarse);
 }
 
 void ReflectInnerX1(std::shared_ptr<MeshBlockData<Real>> &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<X1DIR, BCSide::Inner, BCType::Reflect>(rc, coarse, el);
+  GenericBC<X1DIR, BCSide::Inner, BCType::Reflect, variable_names::any>(rc, coarse);
 }
 
 void ReflectOuterX1(std::shared_ptr<MeshBlockData<Real>> &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<X1DIR, BCSide::Outer, BCType::Reflect>(rc, coarse, el);
+  GenericBC<X1DIR, BCSide::Outer, BCType::Reflect, variable_names::any>(rc, coarse);
 }
 
 void ReflectInnerX2(std::shared_ptr<MeshBlockData<Real>> &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<X2DIR, BCSide::Inner, BCType::Reflect>(rc, coarse, el);
+  GenericBC<X2DIR, BCSide::Inner, BCType::Reflect, variable_names::any>(rc, coarse);
 }
 
 void ReflectOuterX2(std::shared_ptr<MeshBlockData<Real>> &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<X2DIR, BCSide::Outer, BCType::Reflect>(rc, coarse, el);
+  GenericBC<X2DIR, BCSide::Outer, BCType::Reflect, variable_names::any>(rc, coarse);
 }
 
 void ReflectInnerX3(std::shared_ptr<MeshBlockData<Real>> &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<X3DIR, BCSide::Inner, BCType::Reflect>(rc, coarse, el);
+  GenericBC<X3DIR, BCSide::Inner, BCType::Reflect, variable_names::any>(rc, coarse);
 }
 
 void ReflectOuterX3(std::shared_ptr<MeshBlockData<Real>> &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<X3DIR, BCSide::Outer, BCType::Reflect>(rc, coarse, el);
+  GenericBC<X3DIR, BCSide::Outer, BCType::Reflect, variable_names::any>(rc, coarse);
 }
 
 } // namespace BoundaryFunction

src/bvals/boundary_conditions.hpp

@@ -38,6 +38,11 @@ inline TaskStatus ApplyBoundaryConditions(std::shared_ptr<MeshBlockData<Real>> &
   return ApplyBoundaryConditionsOnCoarseOrFine(rc, false);
 }
 
+TaskStatus ApplyBoundaryConditionsMD(std::shared_ptr<MeshData<Real>> &pmd);
+
+TaskStatus ApplyBoundaryConditionsOnCoarseOrFineMD(std::shared_ptr<MeshData<Real>> &pmd,
+                                                   bool coarse);
+
 namespace BoundaryFunction {
 
 void OutflowInnerX1(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse);

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 <functional>
+#include <memory>
+#include <set>
+#include <string>
+#include <vector>
+
+#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 <CoordinateDirection DIR, BCSide SIDE, BCType TYPE, class... var_ts>
+void GenericBC(std::shared_ptr<MeshBlockData<Real>> &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<MetadataFlag> 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<PDOpt> opts;
+  if (coarse) opts = {PDOpt::Coarse};
+  auto desc = MakePackDescriptor<var_ts...>(
+      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<MeshBlock> 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 <CoordinateDirection DIR, BCSide SIDE, BCType TYPE, class... var_ts>
+void GenericBC(std::shared_ptr<MeshBlockData<Real>> &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<DIR, SIDE, TYPE, var_ts...>(rc, coarse, el, val);
+}
+
+} // namespace BoundaryFunction
+} // namespace parthenon
+
+#endif // BVALS_BOUNDARY_CONDITIONS_GENERIC_HPP_

src/interface/sparse_pack_base.cpp

@@ -228,13 +228,21 @@ 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_) {
                       pack_h(0, b, idx) = pv->coarse_s.Get(0, t, u, v);
                     } 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;
 }

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_;
 };