Register boundary conditions by package

example/poisson_gmg/main.cpp

@@ -11,29 +11,10 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
-#include "bvals/boundary_conditions_generic.hpp"
 #include "parthenon_manager.hpp"
 
 #include "poisson_driver.hpp"
 
-using namespace parthenon;
-using namespace parthenon::BoundaryFunction;
-// We need to register FixedFace boundary conditions by hand since they can't
-// be chosen in the parameter input file. FixedFace boundary conditions assume
-// Dirichlet booundary conditions on the face of the domain and linearly extrapolate
-// into the ghosts to ensure the linear reconstruction on the block face obeys the
-// chosen boundary condition. Just setting the ghost zones of CC variables to a fixed
-// value results in poor MG convergence because the effective BC at the face
-// changes with MG level.
-template <CoordinateDirection DIR, BCSide SIDE>
-auto GetBoundaryCondition() {
-  return [](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
-    using namespace parthenon;
-    using namespace parthenon::BoundaryFunction;
-    GenericBC<DIR, SIDE, BCType::FixedFace, variable_names::any>(rc, coarse, 0.0);
-  };
-}
-
 int main(int argc, char *argv[]) {
   using parthenon::ParthenonManager;
   using parthenon::ParthenonStatus;
@@ -56,19 +37,6 @@ int main(int argc, char *argv[]) {
   // Now that ParthenonInit has been called and setup succeeded, the code can now
   // make use of MPI and Kokkos
 
-  // Set boundary conditions
-  pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x1] =
-      GetBoundaryCondition<X1DIR, BCSide::Inner>();
-  pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x2] =
-      GetBoundaryCondition<X2DIR, BCSide::Inner>();
-  pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x3] =
-      GetBoundaryCondition<X3DIR, BCSide::Inner>();
-  pman.app_input->boundary_conditions[parthenon::BoundaryFace::outer_x1] =
-      GetBoundaryCondition<X1DIR, BCSide::Outer>();
-  pman.app_input->boundary_conditions[parthenon::BoundaryFace::outer_x2] =
-      GetBoundaryCondition<X2DIR, BCSide::Outer>();
-  pman.app_input->boundary_conditions[parthenon::BoundaryFace::outer_x3] =
-      GetBoundaryCondition<X3DIR, BCSide::Outer>();
   pman.ParthenonInitPackagesAndMesh();
 
   // This needs to be scoped so that the driver object is destructed before Finalize

example/poisson_gmg/poisson_package.cpp

@@ -18,6 +18,7 @@
 #include <string>
 #include <vector>
 
+#include <bvals/boundary_conditions_generic.hpp>
 #include <coordinates/coordinates.hpp>
 #include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
@@ -34,9 +35,47 @@ using namespace parthenon::package::prelude;
 using parthenon::HostArray1D;
 namespace poisson_package {
 
+using namespace parthenon;
+using namespace parthenon::BoundaryFunction;
+// We need to register FixedFace boundary conditions by hand since they can't
+// be chosen in the parameter input file. FixedFace boundary conditions assume
+// Dirichlet booundary conditions on the face of the domain and linearly extrapolate
+// into the ghosts to ensure the linear reconstruction on the block face obeys the
+// chosen boundary condition. Just setting the ghost zones of CC variables to a fixed
+// value results in poor MG convergence because the effective BC at the face
+// changes with MG level.
+
+// Build type that selects only variables within the poisson namespace. Internal solver
+// variables have the namespace of input variables prepended, so they will also be
+// selected by this type.
+struct any_poisson : public parthenon::variable_names::base_t<true> {
+  template <class... Ts>
+  KOKKOS_INLINE_FUNCTION any_poisson(Ts &&...args)
+      : base_t<true>(std::forward<Ts>(args)...) {}
+  static std::string name() { return "poisson[.].*"; }
+};
+
+template <CoordinateDirection DIR, BCSide SIDE>
+auto GetBC() {
+  return [](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
+    using namespace parthenon;
+    using namespace parthenon::BoundaryFunction;
+    GenericBC<DIR, SIDE, BCType::FixedFace, any_poisson>(rc, coarse, 0.0);
+  };
+}
+
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto pkg = std::make_shared<StateDescriptor>("poisson_package");
 
+  // Set boundary conditions for Poisson variables
+  using BF = parthenon::BoundaryFace;
+  pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetBC<X1DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetBC<X2DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x3].push_back(GetBC<X3DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::outer_x1].push_back(GetBC<X1DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x2].push_back(GetBC<X2DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x3].push_back(GetBC<X3DIR, BCSide::Outer>());
+
   int max_poisson_iterations = pin->GetOrAddInteger("poisson", "max_iterations", 10000);
   pkg->AddParam<>("max_iterations", max_poisson_iterations);
 

src/bvals/boundary_conditions.cpp

@@ -44,6 +44,9 @@ TaskStatus ApplyBoundaryConditionsOnCoarseOrFine(std::shared_ptr<MeshBlockData<R
       PARTHENON_DEBUG_REQUIRE(pmesh->MeshBndryFnctn[i] != nullptr,
                               "boundary function must not be null");
       pmesh->MeshBndryFnctn[i](rc, coarse);
+      for (auto &bnd_func : pmesh->UserBoundaryFunctions[i]) {
+        bnd_func(rc, coarse);
+      }
     }
   }
 

src/interface/state_descriptor.cpp

@@ -425,6 +425,12 @@ StateDescriptor::CreateResolvedStateDescriptor(Packages_t &packages) {
     // sort
     field_tracker.CategorizeCollection(name, field_dict, &field_provider);
     swarm_tracker.CategorizeCollection(name, package->AllSwarms(), &swarm_provider);
+
+    // Add package registered boundary conditions
+    for (int i = 0; i < 6; ++i)
+      state->UserBoundaryFunctions[i].insert(state->UserBoundaryFunctions[i].end(),
+                                             package->UserBoundaryFunctions[i].begin(),
+                                             package->UserBoundaryFunctions[i].end());
   }
 
   // check that dependent variables are provided somewhere

src/interface/state_descriptor.hpp

@@ -44,6 +44,8 @@ class MeshBlockData;
 template <typename T>
 class MeshData;
 
+using BValFunc = std::function<void(std::shared_ptr<MeshBlockData<Real>> &, bool)>;
+
 /// A little container class owning refinement function properties
 /// needed for the state descriptor.
 /// Note using VarID here implies that custom prolongation/restriction
@@ -429,6 +431,8 @@ class StateDescriptor {
 
   friend std::ostream &operator<<(std::ostream &os, const StateDescriptor &sd);
 
+  std::array<std::vector<BValFunc>, 6> UserBoundaryFunctions;
+
  private:
   void InvertControllerMap();
 

src/mesh/mesh.cpp

@@ -406,6 +406,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
 
   resolved_packages = ResolvePackages(packages);
 
+  // Register user defined boundary conditions
+  UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions;
+
   // Setup unique comms for each variable and swarm
   SetupMPIComms();
 
@@ -663,6 +666,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr,
 
   resolved_packages = ResolvePackages(packages);
 
+  // Register user defined boundary conditions
+  UserBoundaryFunctions = resolved_packages->UserBoundaryFunctions;
+
   // Setup unique comms for each variable and swarm
   SetupMPIComms();
 

src/mesh/mesh.hpp

@@ -154,6 +154,7 @@ class Mesh {
   // Boundary Functions
   BValFunc MeshBndryFnctn[6];
   SBValFunc SwarmBndryFnctn[6];
+  std::array<std::vector<BValFunc>, 6> UserBoundaryFunctions;
 
   // defined in either the prob file or default_pgen.cpp in ../pgen/
   std::function<void(Mesh *, ParameterInput *, MeshData<Real> *)> ProblemGenerator =

tst/regression/test_suites/poisson_gmg/parthinput.poisson

@@ -21,14 +21,14 @@ multigrid = true
 nx1 = 64
 x1min = -1.0
 x1max = 1.0
-ix1_bc = user
-ox1_bc = user
+ix1_bc = outflow
+ox1_bc = outflow
 
 nx2 = 64
 x2min = -1.0
 x2max = 1.0
-ix2_bc = user
-ox2_bc = user
+ix2_bc = outflow
+ox2_bc = outflow
 
 nx3 = 1
 x3min = 0.0