Macroscopic Maxwell solver: do not update fields in EB (stair-case approximation) (#2899)

* Macroscopic Maxwell solver: do not update fields in EB
* Correct unused variable
* Add test with macroscopic solver
* Add automated test

6 files changed, 72 insertions, 10 deletions

diff --git a/Examples/Modules/embedded_boundary_cube/analysis_fields.py b/Examples/Modules/embedded_boundary_cube/analysis_fields.py
index b81f72d9d..13bd32d73 100755
--- a/Examples/Modules/embedded_boundary_cube/analysis_fields.py
+++ b/Examples/Modules/embedded_boundary_cube/analysis_fields.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python3
 
 import os
+import re
 import sys
 
 import numpy as np
@@ -41,6 +42,16 @@ filename = sys.argv[1]
 ds = yt.load(filename)
 data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
 
+# Parse test name and check whether this use the macroscopic solver
+# (i.e. solving the equation in a dielectric)
+macroscopic = True if re.search( 'macroscopic', filename ) else False
+
+# Calculate frequency of the mode oscillation
+omega = np.sqrt( h_2 ) * c
+if macroscopic:
+    # Relative permittivity used in this test: epsilon_r = 1.5
+    omega *= 1./np.sqrt(1.5)
+
 t = ds.current_time.to_value()
 
 # Compute the analytic solution
@@ -60,8 +71,7 @@ for i in range(ncells[0]):
                                                                     (-Lx/2 <= x < Lx/2) *
                                                                     (-Ly/2 <= y < Ly/2) *
                                                                     (-Lz/2 <= z < Lz/2) *
-                                                                    np.cos(np.sqrt(2) *
-                                                                    np.pi / Lx * c * t))
+                                                                    np.cos(omega * t))
 
             x = i*dx + lo[0]
             y = j*dy + lo[1]
@@ -72,7 +82,7 @@ for i in range(ncells[0]):
                                    (-Lx/2 <= x < Lx/2) *
                                    (-Ly/2 <= y < Ly/2) *
                                    (-Lz/2 <= z < Lz/2) *
-                                   np.cos(np.sqrt(2) * np.pi / Lx * c * t))
+                                   np.cos(omega * t))
 
 rel_tol_err = 1e-1
 
diff --git a/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json b/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json
new file mode 100644
index 000000000..453b83f5e
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/embedded_boundary_cube_macroscopic.json
@@ -0,0 +1,10 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.005101824310293575,
+    "Bz": 0.005101824310293573,
+    "Ex": 4414725.184731115,
+    "Ey": 0.0,
+    "Ez": 0.0
+  }
+}
\ No newline at end of file
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 90a578824..13c77f6f6 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -2709,6 +2709,23 @@ doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields.py
 
+[embedded_boundary_cube_macroscopic]
+buildDir = .
+inputFile = Examples/Modules/embedded_boundary_cube/inputs_3d
+runtime_params = algo.em_solver_medium=macroscopic macroscopic.epsilon=1.5*8.8541878128e-12  macroscopic.sigma=0 macroscopic.mu=1.25663706212e-06
+dim = 3
+addToCompileString = USE_EB=TRUE
+cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_EB=ON
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Modules/embedded_boundary_cube/analysis_fields.py
+
 [embedded_boundary_cube_2d]
 buildDir = .
 inputFile = Examples/Modules/embedded_boundary_cube/inputs_2d
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index 5b3ef3930..fb774678f 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -108,6 +108,7 @@ class FiniteDifferenceSolver
         void MacroscopicEvolveE ( 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,
                       amrex::Real const dt,
                       std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
 
@@ -241,6 +242,7 @@ class FiniteDifferenceSolver
             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,
             amrex::Real const dt,
             std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
 
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
index ef30b9b51..bc9576cfd 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
@@ -37,6 +37,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
     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,
     amrex::Real const dt,
     std::unique_ptr<MacroscopicProperties> const& macroscopic_properties)
 {
@@ -44,7 +45,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
    // 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, dt, macroscopic_properties);
+    amrex::ignore_unused(Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties);
     amrex::Abort("currently macro E-push does not work for RZ");
 #else
     if (m_do_nodal) {
@@ -55,13 +56,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
 
             MacroscopicEvolveECartesian <CartesianYeeAlgorithm, LaxWendroffAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
+                       ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties);
 
         }
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
 
             MacroscopicEvolveECartesian <CartesianYeeAlgorithm, BackwardEulerAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
+                       ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties);
 
         }
 
@@ -72,12 +73,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
 
             MacroscopicEvolveECartesian <CartesianCKCAlgorithm, LaxWendroffAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
+                       ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties);
 
         } else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
 
             MacroscopicEvolveECartesian <CartesianCKCAlgorithm, BackwardEulerAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
+                       ( Efield, Bfield, Jfield, edge_lengths, dt, macroscopic_properties);
 
         }
 
@@ -96,9 +97,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
     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,
     amrex::Real const dt,
     std::unique_ptr<MacroscopicProperties> const& macroscopic_properties)
 {
+#ifndef AMREX_USE_EB
+    amrex::ignore_unused(edge_lengths);
+#endif
 
     amrex::MultiFab& sigma_mf = macroscopic_properties->getsigma_mf();
     amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf();
@@ -130,6 +135,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
         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
+
         // material prop //
         amrex::Array4<amrex::Real> const& sigma_arr = sigma_mf.array(mfi);
         amrex::Array4<amrex::Real> const& eps_arr = epsilon_mf.array(mfi);
@@ -159,6 +170,10 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
         // Loop over the cells and update the fields
         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
                 // Interpolate conductivity, sigma, to Ex position on the grid
                 amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag,
                                            Ex_stag, macro_cr, i, j, k, scomp);
@@ -174,6 +189,10 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
             },
 
             [=] 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
                 // Interpolate conductivity, sigma, to Ey position on the grid
                 amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag,
                                            Ey_stag, macro_cr, i, j, k, scomp);
@@ -190,6 +209,10 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
             },
 
             [=] 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
                 // Interpolate conductivity, sigma, to Ez position on the grid
                 amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag,
                                            Ez_stag, macro_cr, i, j, k, scomp);
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 2487ffa0c..82ca52c4c 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -741,8 +741,8 @@ void
 WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) {
     if (patch_type == PatchType::fine) {
         m_fdtd_solver_fp[lev]->MacroscopicEvolveE( Efield_fp[lev], Bfield_fp[lev],
-                                             current_fp[lev], a_dt,
-                                             m_macroscopic_properties);
+                                             current_fp[lev], m_edge_lengths[lev],
+                                             a_dt, m_macroscopic_properties);
     }
     else {
         amrex::Abort("Macroscopic EvolveE is not implemented for lev > 0, yet.");