From 1794437dc891c441685b47f51a63303bb46c639c Mon Sep 17 00:00:00 2001
From: Axel Huebl <axel.huebl@plasma.ninja>
Date: Mon, 30 Sep 2024 12:03:41 -0700
Subject: [PATCH] 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.
---
 .../ElectrostaticSolver.cpp                   |  2 +-
 Source/ablastr/fields/PoissonSolver.H         | 39 +++++++++++++------
 2 files changed, 29 insertions(+), 12 deletions(-)

diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
index 1ced0a07152..0b1dca675be 100644
--- 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
index 8b4f9cea9a1..d7eeecead1b 100644
--- 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);