EB: Introduce Runtime Parameter

Docs/source/usage/parameters.rst

@@ -490,13 +490,18 @@ Additional PML parameters
 .. _running-cpp-parameters-eb:
 
 Embedded Boundary Conditions
-----------------------------
+---------------------------
 
 * ``warpx.eb_implicit_function`` (`string`)
     A function of `x`, `y`, `z` that defines the surface of the embedded
     boundary. That surface lies where the function value is 0 ;
     the physics simulation area is where the function value is negative ;
-    the interior of the embeddded boundary is where the function value is positive.
+    the interior of the embedded boundary is where the function value is positive.
+
+* ``eb2.geom_type = stl`` (`string`, default: empty) and ``eb2.stl_file`` (`string` with a filepath to a STL file)
+    Alternatively to defining an embedded boundary via an implicit function, one can load
+    a computer-aided design (CAD) file in the `STL file format <https://en.wikipedia.org/wiki/STL_(file_format)>`__.
+    `See the AMReX documentation for more details <https://amrex-codes.github.io/amrex/docs_html/EB.html>`__.
 
 * ``warpx.eb_potential(x,y,z,t)`` (`string`)
     Gives the value of the electric potential at the surface of the embedded boundary,

Source/BoundaryConditions/PML.H

@@ -146,6 +146,7 @@ public:
          bool do_pml_dive_cleaning, bool do_pml_divb_cleaning,
          const amrex::IntVect& fill_guards_fields,
          const amrex::IntVect& fill_guards_current,
+         bool eb_enabled,
          int max_guard_EB, amrex::Real v_sigma_sb,
          amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(),
          amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector());

Source/BoundaryConditions/PML.cpp

@@ -554,6 +554,7 @@ PML::PML (const int lev, const BoxArray& grid_ba,
           const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning,
           const amrex::IntVect& fill_guards_fields,
           const amrex::IntVect& fill_guards_current,
+          bool eb_enabled,
           int max_guard_EB, const amrex::Real v_sigma_sb,
           const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
     : m_dive_cleaning(do_pml_dive_cleaning),
@@ -563,6 +564,10 @@ PML::PML (const int lev, const BoxArray& grid_ba,
       m_geom(geom),
       m_cgeom(cgeom)
 {
+#ifndef AMREX_USE_EB
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!eb_enabled, "PML: eb_enabled is true but was not compiled in.");
+#endif
+
     // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain or fine patch(es)
     // (instead of extending `ncell` outside of the physical domain or fine patch(es))
     // In order to implement this, we define a new reduced Box Array ensuring that it does not
@@ -669,9 +674,11 @@ PML::PML (const int lev, const BoxArray& grid_ba,
     }
 
 #ifdef AMREX_USE_EB
-    pml_field_factory = amrex::makeEBFabFactory(*geom, ba, dm,
-                                              {max_guard_EB, max_guard_EB, max_guard_EB},
-                                              amrex::EBSupport::full);
+    if (eb_enabled) {
+        pml_field_factory = amrex::makeEBFabFactory(*geom, ba, dm,
+                                                  {max_guard_EB, max_guard_EB, max_guard_EB},
+                                                  amrex::EBSupport::full);
+    }
 #else
     amrex::ignore_unused(max_guard_EB);
     pml_field_factory = std::make_unique<FArrayBoxFactory>();
@@ -703,20 +710,22 @@ PML::PML (const int lev, const BoxArray& grid_ba,
     WarpX::AllocInitMultiFab(pml_j_fp[2], ba_jz, dm, 1, ngb, lev, "pml_j_fp[z]", 0.0_rt);
 
 #ifdef AMREX_USE_EB
-    const amrex::IntVect max_guard_EB_vect = amrex::IntVect(max_guard_EB);
-    WarpX::AllocInitMultiFab(pml_edge_lengths[0], ba_Ex, dm, WarpX::ncomps, max_guard_EB_vect, lev, "pml_edge_lengths[x]", 0.0_rt);
-    WarpX::AllocInitMultiFab(pml_edge_lengths[1], ba_Ey, dm, WarpX::ncomps, max_guard_EB_vect, lev, "pml_edge_lengths[y]", 0.0_rt);
-    WarpX::AllocInitMultiFab(pml_edge_lengths[2], ba_Ez, dm, WarpX::ncomps, max_guard_EB_vect, lev, "pml_edge_lengths[z]", 0.0_rt);
+    if (eb_enabled) {
+        const amrex::IntVect max_guard_EB_vect = amrex::IntVect(max_guard_EB);
+        WarpX::AllocInitMultiFab(pml_edge_lengths[0], ba_Ex, dm, WarpX::ncomps, max_guard_EB_vect, lev, "pml_edge_lengths[x]", 0.0_rt);
+        WarpX::AllocInitMultiFab(pml_edge_lengths[1], ba_Ey, dm, WarpX::ncomps, max_guard_EB_vect, lev, "pml_edge_lengths[y]", 0.0_rt);
+        WarpX::AllocInitMultiFab(pml_edge_lengths[2], ba_Ez, dm, WarpX::ncomps, max_guard_EB_vect, lev, "pml_edge_lengths[z]", 0.0_rt);
 
-    if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee ||
-        WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::CKC ||
-        WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) {
+        if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::Yee ||
+            WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::CKC ||
+            WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT) {
 
-        auto const eb_fact = fieldEBFactory();
+            auto const eb_fact = fieldEBFactory();
 
-        WarpX::ComputeEdgeLengths(pml_edge_lengths, eb_fact);
-        WarpX::ScaleEdges(pml_edge_lengths, WarpX::CellSize(lev));
+            WarpX::ComputeEdgeLengths(pml_edge_lengths, eb_fact);
+            WarpX::ScaleEdges(pml_edge_lengths, WarpX::CellSize(lev));
 
+        }
     }
 #endif
 

Source/BoundaryConditions/WarpXEvolvePML.cpp

@@ -269,13 +269,16 @@ WarpX::DampJPML (int lev, PatchType patch_type)
             const Real* sigma_star_cumsum_fac_j_z = sigba[mfi].sigma_star_cumsum_fac[1].data();
 #endif
 
-#ifdef AMREX_USE_EB
-            const auto& pml_edge_lenghts = pml[lev]->Get_edge_lengths();
-
-            auto const& pml_lxfab = pml_edge_lenghts[0]->array(mfi);
-            auto const& pml_lyfab = pml_edge_lenghts[1]->array(mfi);
-            auto const& pml_lzfab = pml_edge_lenghts[2]->array(mfi);
-#endif
+            amrex::Array4<amrex::Real> pml_lxfab, pml_lyfab, pml_lzfab;
+            if (m_eb_enabled) {
+                const auto &pml_edge_lenghts = pml[lev]->Get_edge_lengths();
+
+                pml_lxfab = pml_edge_lenghts[0]->array(mfi);
+                pml_lyfab = pml_edge_lenghts[1]->array(mfi);
+                pml_lzfab = pml_edge_lenghts[2]->array(mfi);
+            } else {
+                amrex::ignore_unused(pml_lxfab, pml_lyfab, pml_lzfab);
+            }
 
             const Box& tjx  = mfi.tilebox( pml_j[0]->ixType().toIntVect() );
             const Box& tjy  = mfi.tilebox( pml_j[1]->ixType().toIntVect() );
@@ -299,29 +302,24 @@ WarpX::DampJPML (int lev, PatchType patch_type)
             int const zs_lo = sigba[mfi].sigma_star_cumsum_fac[1].lo();
 #endif
 
+            bool const eb_enabled = m_eb_enabled;
             amrex::ParallelFor( tjx, tjy, tjz,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-#ifdef AMREX_USE_EB
-                    if(pml_lxfab(i, j, k) <= 0) return;
-#endif
+                    if (eb_enabled && pml_lxfab(i, j, k) <= 0) return;
 
                     damp_jx_pml(i, j, k, pml_jxfab, sigma_star_cumsum_fac_j_x,
                                 sigma_cumsum_fac_j_y, sigma_cumsum_fac_j_z,
                                 xs_lo,y_lo, z_lo);
                 },
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-#ifdef AMREX_USE_EB
-                    if(pml_lyfab(i, j, k) <= 0) return;
-#endif
+                    if (eb_enabled && pml_lyfab(i, j, k) <= 0) return;
 
                     damp_jy_pml(i, j, k, pml_jyfab, sigma_cumsum_fac_j_x,
                                 sigma_star_cumsum_fac_j_y, sigma_cumsum_fac_j_z,
                                 x_lo,ys_lo, z_lo);
                 },
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-#ifdef AMREX_USE_EB
-                    if(pml_lzfab(i, j, k)<=0) return;
-#endif
+                    if (eb_enabled && pml_lzfab(i, j, k)<=0) return;
 
                     damp_jz_pml(i, j, k, pml_jzfab, sigma_cumsum_fac_j_x,
                                 sigma_cumsum_fac_j_y, sigma_star_cumsum_fac_j_z,

Source/Diagnostics/BoundaryScrapingDiagnostics.cpp

@@ -6,6 +6,7 @@
  */
 
 #include "BoundaryScrapingDiagnostics.H"
+#include "EmbeddedBoundary/Enabled.H"
 #include "ComputeDiagFunctors/ComputeDiagFunctor.H"
 #include "Diagnostics/Diagnostics.H"
 #include "Diagnostics/FlushFormats/FlushFormat.H"
@@ -39,11 +40,11 @@ BoundaryScrapingDiagnostics::ReadParameters ()
 
     // num_buffers corresponds to the number of boundaries
     // (upper/lower domain boundary in each dimension)
-    // + the EB boundary if available
     m_num_buffers = AMREX_SPACEDIM*2;
-#ifdef AMREX_USE_EB
-    m_num_buffers += 1;
-#endif
+
+    // + the EB boundary if available
+    bool const eb_enabled = EB::enabled();
+    if (eb_enabled) { m_num_buffers += 1; }
 
     // Do a few checks
 #ifndef WARPX_USE_OPENPMD

Source/Diagnostics/ReducedDiags/ChargeOnEB.cpp

@@ -8,6 +8,7 @@
 #include "ChargeOnEB.H"
 
 #include "Diagnostics/ReducedDiags/ReducedDiags.H"
+#include "EmbeddedBoundary/Enabled.H"
 #include "FieldSolver/Fields.H"
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXConst.H"
@@ -24,11 +25,13 @@
 
 #include <algorithm>
 #include <fstream>
+#include <stdexcept>
 #include <vector>
 
 using namespace amrex;
 using namespace warpx::fields;
 
+
 // constructor
 ChargeOnEB::ChargeOnEB (const std::string& rd_name)
 : ReducedDiags{rd_name}
@@ -44,6 +47,10 @@ ChargeOnEB::ChargeOnEB (const std::string& rd_name)
         "ChargeOnEB reduced diagnostics only works when compiling with EB support");
 #endif
 
+    if (!EB::enabled()) {
+        throw std::runtime_error("ChargeOnEB reduced diagnostics only works when EBs are enabled at runtime");
+    }
+
     // resize data array
     m_data.resize(1, 0.0_rt);
 
@@ -87,6 +94,9 @@ void ChargeOnEB::ComputeDiags (const int step)
     // Judge whether the diags should be done
     if (!m_intervals.contains(step+1)) { return; }
 
+    if (!EB::enabled()) {
+        throw std::runtime_error("ComputeDiags only works when EBs are enabled at runtime");
+    }
 #if ((defined WARPX_DIM_3D) && (defined AMREX_USE_EB))
     // get a reference to WarpX instance
     auto & warpx = WarpX::GetInstance();

Source/Diagnostics/WarpXIO.cpp

@@ -393,7 +393,7 @@ WarpX::InitFromCheckpoint ()
         }
     }
 
-    InitializeEBGridData(maxLevel());
+    if (m_eb_enabled) { InitializeEBGridData(maxLevel()); }
 
     // Initialize particles
     mypc->AllocData();

Source/EmbeddedBoundary/CMakeLists.txt

@@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS)
     warpx_set_suffix_dims(SD ${D})
     target_sources(lib_${SD}
       PRIVATE
+        Enabled.cpp
         WarpXInitEB.cpp
         WarpXFaceExtensions.cpp
         WarpXFaceInfoBox.H

Source/EmbeddedBoundary/DistanceToEB.H

@@ -14,8 +14,6 @@
 #include <AMReX_RealVect.H>
 #include <AMReX_Array.H>
 
-#ifdef AMREX_USE_EB
-
 namespace DistanceToEB
 {
 
@@ -120,13 +118,11 @@ amrex::RealVect interp_normal (int i, int j, int k, const amrex::Real W[AMREX_SP
 
 #else
     amrex::ignore_unused(i, j, k, ic, jc, kc, W, Wc, phi, dxi);
-    amrex::RealVect normal{0.0, 0.0};
+    amrex::RealVect normal({0.0});
     WARPX_ABORT_WITH_MESSAGE("Error: interp_distance not yet implemented in 1D");
 
 #endif
     return normal;
 }
 }
-
-#endif // AMREX_USE_EB
 #endif // WARPX_DISTANCETOEB_H_

Source/EmbeddedBoundary/Enabled.H

@@ -0,0 +1,18 @@
+/* Copyright 2024 Axel Huebl
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_EB_ENABLED_H_
+#define WARPX_EB_ENABLED_H_
+
+#include <AMReX_ParmParse.H>
+
+namespace EB
+{
+    /** Are embedded boundaries enabled? */
+    bool enabled ();
+
+} // namespace EB
+#endif // WARPX_EB_ENABLED_H_

Source/EmbeddedBoundary/Enabled.cpp

@@ -0,0 +1,45 @@
+/* Copyright 2024 Axel Huebl
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "Enabled.H"
+
+#ifdef AMREX_USE_EB
+#include <AMReX_ParmParse.H>
+#endif
+#if defined(AMREX_USE_EB) && defined(WARPX_DIM_RZ)
+#include <stdexcept>
+#endif
+
+
+namespace EB
+{
+    bool enabled ()
+    {
+#ifndef AMREX_USE_EB
+        return false;
+#else
+        amrex::ParmParse const pp_warpx("warpx");
+        amrex::ParmParse const pp_eb2("eb2");
+
+        // test various runtime options to enable EBs
+        std::string eb_implicit_function;
+        bool eb_enabled = pp_warpx.query("eb_implicit_function", eb_implicit_function);
+
+        // https://amrex-codes.github.io/amrex/docs_html/EB.html
+        std::string eb_stl;
+        eb_enabled |= pp_eb2.query("geom_type", eb_stl);
+
+#if defined(WARPX_DIM_RZ)
+        if (eb_enabled) {
+            throw std::runtime_error("RZ Geometry does not yet support EBs, but EBs are enabled in runtime inputs.");
+        }
+#endif
+
+        return eb_enabled;
+#endif
+    }
+
+} // namespace EB