From 549d8a0a89f6d9c40444acdc7d11f16d5a87858f Mon Sep 17 00:00:00 2001
From: Lorenzo Giacomel <47607756+lgiacome@users.noreply.github.com>
Date: Thu, 23 Dec 2021 16:04:44 +0100
Subject: 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)
---
 .../FiniteDifferenceSolver/EvolveBPML.cpp          | 29 +++++++++++++++++++---
 1 file changed, 25 insertions(+), 4 deletions(-)

(limited to 'Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp')

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
-- 
cgit v1.2.3