Poisson computePhi: Simplify Boundary Handler

Move the boundary handler to become an optional argument, which otherwise defaults to Dirichlet conditions, e.g., in non-EB cases. This simplifies the ImpactX implementation and fixes a linker issue with CUDA for ImpactX.
ax3l · Sep 30, 2024 · 1794437 · 1794437
1 parent b1aa846
commit 1794437
Show file tree

Hide file tree

Showing 2 changed files with 29 additions and 12 deletions.
diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
@@ -201,12 +201,12 @@ ElectrostaticSolver::computePhi (
         warpx.DistributionMap(),
         warpx.boxArray(),
         WarpX::grid_type,
-        *m_poisson_boundary_handler,
         is_solver_igf_on_lev0,
         EB::enabled(),
         WarpX::do_single_precision_comms,
         warpx.refRatio(),
         post_phi_calculation,
+        *m_poisson_boundary_handler,
         warpx.gett_new(0),
         eb_farray_box_factory
     );

diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H
@@ -164,8 +164,8 @@ inline void interpolatePhiBetweenLevels (
  *   \vec{\nabla}^2 r \phi - (\vec{\beta}\cdot\vec{\nabla})^2 r \phi = -\frac{r \rho}{\epsilon_0}
  * \f]
  *
- * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler
  * \tparam T_PostPhiCalculationFunctor a calculation per level directly after phi was calculated
+ * \tparam T_BoundaryHandler handler for boundary conditions, for example @see ElectrostaticSolver::PoissonBoundaryHandler (EB ONLY)
  * \tparam T_FArrayBoxFactory usually nothing or an amrex::EBFArrayBoxFactory (EB ONLY)
  * \param[in] rho The charge density a given species
  * \param[out] phi The potential to be computed by this function
@@ -188,8 +188,8 @@ inline void interpolatePhiBetweenLevels (
  * \param[in] eb_farray_box_factory a factory for field data, @see amrex::EBFArrayBoxFactory; required for embedded boundaries (default: none)
  */
 template<
-    typename T_BoundaryHandler,
     typename T_PostPhiCalculationFunctor = std::nullopt_t,
+    typename T_BoundaryHandler = std::nullopt_t,
     typename T_FArrayBoxFactory = void
 >
 void
@@ -205,12 +205,12 @@ computePhi (
     amrex::Vector<amrex::DistributionMapping> const& dmap,
     amrex::Vector<amrex::BoxArray> const& grids,
     utils::enums::GridType grid_type,
-    T_BoundaryHandler const boundary_handler,
     bool is_solver_igf_on_lev0,
     bool eb_enabled = false,
     bool do_single_precision_comms = false,
     std::optional<amrex::Vector<amrex::IntVect> > rel_ref_ratio = std::nullopt,
     [[maybe_unused]] T_PostPhiCalculationFunctor post_phi_calculation = std::nullopt,
+    [[maybe_unused]] T_BoundaryHandler const boundary_handler = std::nullopt, // only used for EB
     [[maybe_unused]] std::optional<amrex::Real const> current_time = std::nullopt, // only used for EB
     [[maybe_unused]] std::optional<amrex::Vector<T_FArrayBoxFactory const *> > eb_farray_box_factory = std::nullopt // only used for EB
 )
@@ -349,12 +349,18 @@ computePhi (
 #endif
 #if defined(AMREX_USE_EB)
             if (eb_enabled) {
-                // if the EB potential only depends on time, the potential can be passed
-                // as a float instead of a callable
-                if (boundary_handler.phi_EB_only_t) {
-                    linop_nodelap->setEBDirichlet(boundary_handler.potential_eb_t(current_time.value()));
-                } else {
-                    linop_nodelap->setEBDirichlet(boundary_handler.getPhiEB(current_time.value()));
+                if constexpr (!std::is_same_v<T_BoundaryHandler, std::nullopt_t>) {
+                    // if the EB potential only depends on time, the potential can be passed
+                    // as a float instead of a callable
+                    if (boundary_handler.phi_EB_only_t) {
+                        linop_nodelap->setEBDirichlet(boundary_handler.potential_eb_t(current_time.value()));
+                    } else {
+                        linop_nodelap->setEBDirichlet(boundary_handler.getPhiEB(current_time.value()));
+                    }
+                } else
+                {
+                    ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0,
+                        "EB Poisson solver enabled but no 'boundary_handler' passed!");
                 }
             }
 #endif
@@ -372,9 +378,20 @@ computePhi (
             linop = std::move(linop_tenslap);
         }
 
-        // Solve the Poisson equation
-        linop->setDomainBC(boundary_handler.lobc, boundary_handler.hibc);
+        // Level 0 domain boundary
+        if constexpr (std::is_same_v<T_BoundaryHandler, std::nullopt_t>) {
+            amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> const lobc = {AMREX_D_DECL(
+                amrex::LinOpBCType::Dirichlet,
+                amrex::LinOpBCType::Dirichlet,
+                amrex::LinOpBCType::Dirichlet
+            )};
+            amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM> const hibc = lobc;
+            linop->setDomainBC(lobc, hibc);
+        } else {
+            linop->setDomainBC(boundary_handler.lobc, boundary_handler.hibc);
+        }
 
+        // Solve the Poisson equation
         amrex::MLMG mlmg(*linop); // actual solver defined here
         mlmg.setVerbose(verbosity);
         mlmg.setMaxIter(max_iters);