Adding staircased EB in the PMLs (#2693)

* Adding EB to PMLs

* Fix to guard cells

* Initializing EB data only if AMREX_USE_EB

* Not compiling at all EB initialization when EB is off

* Bug fix (wrong order of PML init arguments)

3 files changed, 57 insertions, 8 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
index 5b182812a..f2eb555e1 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
@@ -42,26 +42,27 @@ using namespace amrex;
 void FiniteDifferenceSolver::EvolveBPML (
     std::array< amrex::MultiFab*, 3 > Bfield,
     std::array< amrex::MultiFab*, 3 > const Efield,
+    std::array< amrex::MultiFab*, 3 > const face_areas,
     amrex::Real const dt,
     const bool dive_cleaning) {
 
    // Select algorithm (The choice of algorithm is a runtime option,
    // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
-    amrex::ignore_unused(Bfield, Efield, dt, dive_cleaning);
+    amrex::ignore_unused(Bfield, Efield, dt, dive_cleaning, face_areas);
     amrex::Abort("PML are not implemented in cylindrical geometry.");
 #else
     if (m_do_nodal) {
 
-        EvolveBPMLCartesian <CartesianNodalAlgorithm> (Bfield, Efield, dt, dive_cleaning);
+        EvolveBPMLCartesian <CartesianNodalAlgorithm> (Bfield, Efield, face_areas, dt, dive_cleaning);
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) {
 
-        EvolveBPMLCartesian <CartesianYeeAlgorithm> (Bfield, Efield, dt, dive_cleaning);
+        EvolveBPMLCartesian <CartesianYeeAlgorithm> (Bfield, Efield, face_areas, dt, dive_cleaning);
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveBPMLCartesian <CartesianCKCAlgorithm> (Bfield, Efield, dt, dive_cleaning);
+        EvolveBPMLCartesian <CartesianCKCAlgorithm> (Bfield, Efield, face_areas, dt, dive_cleaning);
 
     } else {
         amrex::Abort("EvolveBPML: Unknown algorithm");
@@ -76,6 +77,7 @@ template<typename T_Algo>
 void FiniteDifferenceSolver::EvolveBPMLCartesian (
     std::array< amrex::MultiFab*, 3 > Bfield,
     std::array< amrex::MultiFab*, 3 > const Efield,
+    std::array< amrex::MultiFab*, 3 > const face_areas,
     amrex::Real const dt,
     const bool dive_cleaning) {
 
@@ -93,6 +95,12 @@ void FiniteDifferenceSolver::EvolveBPMLCartesian (
         Array4<Real> const& Ey = Efield[1]->array(mfi);
         Array4<Real> const& Ez = Efield[2]->array(mfi);
 
+#ifdef AMREX_USE_EB
+        Array4<Real> const& Sx = face_areas[0]->array(mfi);
+        Array4<Real> const& Sy = face_areas[1]->array(mfi);
+        Array4<Real> const& Sz = face_areas[2]->array(mfi);
+#endif
+
         // Extract stencil coefficients
         Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
         int const n_coefs_x = m_stencil_coefs_x.size();
@@ -110,6 +118,9 @@ void FiniteDifferenceSolver::EvolveBPMLCartesian (
         amrex::ParallelFor(tbx, tby, tbz,
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                if(Sx(i, j, k) <= 0) return;
+#endif
 
                 amrex::Real UpwardDz_Ey_yy = 0._rt;
                 amrex::Real UpwardDy_Ez_zz = 0._rt;
@@ -131,6 +142,9 @@ void FiniteDifferenceSolver::EvolveBPMLCartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                if(Sy(i, j, k) <= 0) return;
+#endif
 
                 amrex::Real UpwardDx_Ez_zz = 0._rt;
                 amrex::Real UpwardDz_Ex_xx = 0._rt;
@@ -152,6 +166,9 @@ void FiniteDifferenceSolver::EvolveBPMLCartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                if(Sz(i, j, k) <= 0) return;
+#endif
 
                 amrex::Real UpwardDy_Ex_xx = 0._rt;
                 amrex::Real UpwardDx_Ey_yy = 0._rt;
@@ -176,6 +193,10 @@ void FiniteDifferenceSolver::EvolveBPMLCartesian (
 
     }
 
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(face_areas);
+#endif
+
 }
 
 #endif // corresponds to ifndef WARPX_DIM_RZ
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
index b8c43190f..8258ceb4b 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
@@ -46,6 +46,7 @@ void FiniteDifferenceSolver::EvolveEPML (
     std::array< amrex::MultiFab*, 3 > Efield,
     std::array< amrex::MultiFab*, 3 > const Bfield,
     std::array< amrex::MultiFab*, 3 > const Jfield,
+    std::array< amrex::MultiFab*, 3 > const edge_lengths,
     amrex::MultiFab* const Ffield,
     MultiSigmaBox const& sigba,
     amrex::Real const dt, bool pml_has_particles ) {
@@ -53,23 +54,23 @@ void FiniteDifferenceSolver::EvolveEPML (
    // Select algorithm (The choice of algorithm is a runtime option,
    // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
-    amrex::ignore_unused(Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles);
+    amrex::ignore_unused(Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles, edge_lengths);
     amrex::Abort("PML are not implemented in cylindrical geometry.");
 #else
     if (m_do_nodal) {
 
         EvolveEPMLCartesian <CartesianNodalAlgorithm> (
-            Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
+            Efield, Bfield, Jfield, edge_lengths, Ffield, sigba, dt, pml_has_particles );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee || m_fdtd_algo == MaxwellSolverAlgo::ECT) {
 
         EvolveEPMLCartesian <CartesianYeeAlgorithm> (
-            Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
+            Efield, Bfield, Jfield,  edge_lengths, Ffield, sigba, dt, pml_has_particles );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
         EvolveEPMLCartesian <CartesianCKCAlgorithm> (
-            Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
+            Efield, Bfield, Jfield,  edge_lengths, Ffield, sigba, dt, pml_has_particles );
 
     } else {
         amrex::Abort("EvolveEPML: Unknown algorithm");
@@ -85,6 +86,7 @@ void FiniteDifferenceSolver::EvolveEPMLCartesian (
     std::array< amrex::MultiFab*, 3 > Efield,
     std::array< amrex::MultiFab*, 3 > const Bfield,
     std::array< amrex::MultiFab*, 3 > const Jfield,
+    std::array< amrex::MultiFab*, 3 > const edge_lengths,
     amrex::MultiFab* const Ffield,
     MultiSigmaBox const& sigba,
     amrex::Real const dt, bool pml_has_particles ) {
@@ -105,6 +107,12 @@ void FiniteDifferenceSolver::EvolveEPMLCartesian (
         Array4<Real> const& By = Bfield[1]->array(mfi);
         Array4<Real> const& Bz = Bfield[2]->array(mfi);
 
+#ifdef AMREX_USE_EB
+        Array4<Real> const& lx = edge_lengths[0]->array(mfi);
+        Array4<Real> const& ly = edge_lengths[1]->array(mfi);
+        Array4<Real> const& lz = edge_lengths[2]->array(mfi);
+#endif
+
         // Extract stencil coefficients
         Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
         int const n_coefs_x = m_stencil_coefs_x.size();
@@ -122,6 +130,10 @@ void FiniteDifferenceSolver::EvolveEPMLCartesian (
         amrex::ParallelFor(tex, tey, tez,
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                if(lx(i, j, k) <= 0) return;
+#endif
+
                 Ex(i, j, k, PMLComp::xz) -= c2 * dt * (
                     T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, PMLComp::yx)
                   + T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) );
@@ -131,6 +143,10 @@ void FiniteDifferenceSolver::EvolveEPMLCartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                if(ly(i, j, k)<=0) return;
+#endif
+
                 Ey(i, j, k, PMLComp::yx) -= c2 * dt * (
                     T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, PMLComp::zx)
                   + T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, PMLComp::zy) );
@@ -140,6 +156,10 @@ void FiniteDifferenceSolver::EvolveEPMLCartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                if(lz(i, j, k) <= 0) return;
+#endif
+
                 Ez(i, j, k, PMLComp::zy) -= c2 * dt * (
                     T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, PMLComp::xy)
                   + T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) );
@@ -219,6 +239,10 @@ void FiniteDifferenceSolver::EvolveEPMLCartesian (
 
     }
 
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(edge_lengths);
+#endif
+
 }
 
 #endif // corresponds to ifndef WARPX_DIM_RZ
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index ab9b31e4d..f1f9b4837 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -107,12 +107,14 @@ class FiniteDifferenceSolver
 
         void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield,
                       std::array< amrex::MultiFab*, 3 > const Efield,
+                      std::array< amrex::MultiFab*, 3 > const face_areas,
                       amrex::Real const dt,
                       const bool dive_cleaning);
 
        void EvolveEPML ( std::array< amrex::MultiFab*, 3 > Efield,
                       std::array< amrex::MultiFab*, 3 > const Bfield,
                       std::array< amrex::MultiFab*, 3 > const Jfield,
+                      std::array< amrex::MultiFab*, 3 > const edge_lengths,
                       amrex::MultiFab* const Ffield,
                       MultiSigmaBox const& sigba,
                       amrex::Real const dt, bool pml_has_particles );
@@ -242,6 +244,7 @@ class FiniteDifferenceSolver
         void EvolveBPMLCartesian (
             std::array< amrex::MultiFab*, 3 > Bfield,
             std::array< amrex::MultiFab*, 3 > const Efield,
+            std::array< amrex::MultiFab*, 3 > const face_areas,
             amrex::Real const dt,
             const bool dive_cleaning);
 
@@ -250,6 +253,7 @@ class FiniteDifferenceSolver
             std::array< amrex::MultiFab*, 3 > Efield,
             std::array< amrex::MultiFab*, 3 > const Bfield,
             std::array< amrex::MultiFab*, 3 > const Jfield,
+            std::array< amrex::MultiFab*, 3 > const edge_lengths,
             amrex::MultiFab* const Ffield,
             MultiSigmaBox const& sigba,
             amrex::Real const dt, bool pml_has_particles );