Staircased embedded boundaries in the YEE solver (#1881)

* Added staircased embedded boundaris to the YEE solver

* adding spherical resonating cavity test

* adding functions for fields initialization

* style adjustments

* fixing tabs

* fixed name of analysis script

* fixed name of analysis script

* fixed a few wrong preprocessor directives

* workaround for missing boost

* Revert "workaround for missing boost"

This reverts commit 601f9eb2ec6f8c2100304379b2bea1c6cf9d1851.

* another workaround for missing boost

* getting rid of boost by depending on c++17

* Removed a few unused variables

* adding USE_EB to addToCompileString for EB testing

* removed tabs

* fixing the inputs name for EB sphere test

* shortened the test

* zero padding the names of the images

* adjusted two for loops

* removed some unused variables

* improving the fields initialization

* removed the sphere test and implemented the cube test

* fixed edges lengths computation and added comments

* Fixed the case of all_regular geometries

* fixing a bug that was breaking some tests

* adding test folder

* fixed the default values for the EB cube test

* simplified the analysis script

* fixed cubic resonator default results

* inputting the plot file name from command line

* fixing the diag name

* Fixed a bug in edges initialization

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Adding comments to the staircased yee solver (thanks Remi)

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* fixed the cube resonator test

* removed an unused import

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

1 files changed, 23 insertions, 3 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
index cd1e3bdf9..3cb433e97 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
@@ -28,6 +28,7 @@ void FiniteDifferenceSolver::EvolveE (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
     int lev, amrex::Real const dt ) {
 
@@ -41,15 +42,15 @@ void FiniteDifferenceSolver::EvolveE (
 #else
     if (m_do_nodal) {
 
-        EvolveECartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
+        EvolveECartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        EvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
+        EvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
+        EvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, edge_lengths, Ffield, lev, dt );
 
 #endif
     } else {
@@ -66,6 +67,7 @@ void FiniteDifferenceSolver::EvolveECartesian (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& edge_lengths,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
     int lev, amrex::Real const dt ) {
 
@@ -94,6 +96,12 @@ void FiniteDifferenceSolver::EvolveECartesian (
         Array4<Real> const& jy = Jfield[1]->array(mfi);
         Array4<Real> const& jz = Jfield[2]->array(mfi);
 
+#ifdef AMREX_USE_EB
+        amrex::Array4<amrex::Real> const& lx = edge_lengths[0]->array(mfi);
+        amrex::Array4<amrex::Real> const& ly = edge_lengths[1]->array(mfi);
+        amrex::Array4<amrex::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();
@@ -111,6 +119,10 @@ void FiniteDifferenceSolver::EvolveECartesian (
         amrex::ParallelFor(tex, tey, tez,
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip field push if this cell is fully covered by embedded boundaries
+                if (lx(i, j, k) <= 0) return;
+#endif
                 Ex(i, j, k) += c2 * dt * (
                     - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k)
                     + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k)
@@ -118,6 +130,10 @@ void FiniteDifferenceSolver::EvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip field push if this cell is fully covered by embedded boundaries
+                if (ly(i,j,k) <= 0) return;
+#endif
                 Ey(i, j, k) += c2 * dt * (
                     - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k)
                     + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k)
@@ -125,6 +141,10 @@ void FiniteDifferenceSolver::EvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
+#ifdef AMREX_USE_EB
+                // Skip field push if this cell is fully covered by embedded boundaries
+                if (lz(i,j,k) <= 0) return;
+#endif
                 Ez(i, j, k) += c2 * dt * (
                     - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k)
                     + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k)