Cell Center Macroscopic Properties (#2530)

* Cell Center Macroscopic Properties

* Commit Suggestions from PR Review

* Fix Error for 2D

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

1 files changed, 53 insertions, 41 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
index 2647a1c08..ef30b9b51 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
@@ -38,13 +38,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
     amrex::Real const dt,
-    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties, int const lev)
+    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties)
 {
 
    // 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, lev);
+    amrex::ignore_unused(Efield, Bfield, Jfield, dt, macroscopic_properties);
     amrex::Abort("currently macro E-push does not work for RZ");
 #else
     if (m_do_nodal) {
@@ -55,13 +55,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
 
             MacroscopicEvolveECartesian <CartesianYeeAlgorithm, LaxWendroffAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
 
         }
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
 
             MacroscopicEvolveECartesian <CartesianYeeAlgorithm, BackwardEulerAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
 
         }
 
@@ -72,12 +72,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
 
             MacroscopicEvolveECartesian <CartesianCKCAlgorithm, LaxWendroffAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
 
         } else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
 
             MacroscopicEvolveECartesian <CartesianCKCAlgorithm, BackwardEulerAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties);
 
         }
 
@@ -97,22 +97,21 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
     amrex::Real const dt,
-    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties, int const lev)
+    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties)
 {
 
+    amrex::MultiFab& sigma_mf = macroscopic_properties->getsigma_mf();
+    amrex::MultiFab& epsilon_mf = macroscopic_properties->getepsilon_mf();
+    amrex::MultiFab& mu_mf = macroscopic_properties->getmu_mf();
+
+    // Index type required for calling CoarsenIO::Interp to interpolate macroscopic
+    // properties from their respective staggering to the Ex, Ey, Ez locations
+    amrex::GpuArray<int, 3> const& sigma_stag = macroscopic_properties->sigma_IndexType;
+    amrex::GpuArray<int, 3> const& epsilon_stag = macroscopic_properties->epsilon_IndexType;
+    amrex::GpuArray<int, 3> const& macro_cr     = macroscopic_properties->macro_cr_ratio;
     amrex::GpuArray<int, 3> const& Ex_stag = macroscopic_properties->Ex_IndexType;
     amrex::GpuArray<int, 3> const& Ey_stag = macroscopic_properties->Ey_IndexType;
     amrex::GpuArray<int, 3> const& Ez_stag = macroscopic_properties->Ez_IndexType;
-    amrex::GpuArray<int, 3> const& Bx_stag = macroscopic_properties->Bx_IndexType;
-    amrex::GpuArray<int, 3> const& By_stag = macroscopic_properties->By_IndexType;
-    amrex::GpuArray<int, 3> const& Bz_stag = macroscopic_properties->Bz_IndexType;
-
-    const auto getSigma = GetSigmaMacroparameter();
-    const auto getEpsilon = GetEpsilonMacroparameter();
-    const auto getMu = GetMuMacroparameter();
-    auto &warpx = WarpX::GetInstance();
-    const auto problo = warpx.Geom(lev).ProbLoArray();
-    const auto dx = warpx.Geom(lev).CellSizeArray();
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef AMREX_USE_OMP
@@ -131,6 +130,11 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
         Array4<Real> const& jy = Jfield[1]->array(mfi);
         Array4<Real> const& jz = Jfield[2]->array(mfi);
 
+        // material prop //
+        amrex::Array4<amrex::Real> const& sigma_arr = sigma_mf.array(mfi);
+        amrex::Array4<amrex::Real> const& eps_arr = epsilon_mf.array(mfi);
+        amrex::Array4<amrex::Real> const& mu_arr = mu_mf.array(mfi);
+
         // 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();
@@ -139,24 +143,30 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
         Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr();
         int const n_coefs_z = m_stencil_coefs_z.size();
 
-        FieldAccessorMacroscopic<GetMuMacroparameter> const Hx(Bx, getMu, Bx_stag, problo, dx);
-        FieldAccessorMacroscopic<GetMuMacroparameter> const Hy(By, getMu, By_stag, problo, dx);
-        FieldAccessorMacroscopic<GetMuMacroparameter> const Hz(Bz, getMu, Bz_stag, problo, dx);
+        // This functor computes Hx = Bx/mu
+        // Note that mu is cell-centered here and will be interpolated/averaged
+        // to the location where the B-field and H-field are defined
+        FieldAccessorMacroscopic const Hx(Bx, mu_arr);
+        FieldAccessorMacroscopic const Hy(By, mu_arr);
+        FieldAccessorMacroscopic const Hz(Bz, mu_arr);
 
         // Extract tileboxes for which to loop
         Box const& tex  = mfi.tilebox(Efield[0]->ixType().toIntVect());
         Box const& tey  = mfi.tilebox(Efield[1]->ixType().toIntVect());
         Box const& tez  = mfi.tilebox(Efield[2]->ixType().toIntVect());
+        // starting component to interpolate macro properties to Ex, Ey, Ez locations
+        const int scomp = 0;
         // Loop over the cells and update the fields
         amrex::ParallelFor(tex, tey, tez,
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                amrex::Real x, y, z;
-                WarpXUtilAlgo::getCellCoordinates (i, j, k, Ex_stag, problo, dx,
-                                                   x, y, z );
-                amrex::Real const sigma = getSigma(x, y, z);
-                amrex::Real const epsilon = getEpsilon(x, y, z);
-                amrex::Real alpha = T_MacroAlgo::alpha( sigma, epsilon, dt);
-                amrex::Real beta = T_MacroAlgo::beta( sigma, epsilon, dt);
+                // 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);
+                // Interpolated permittivity, epsilon, to Ex position on the grid
+                amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag,
+                                           Ex_stag, macro_cr, i, j, k, scomp);
+                amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt);
+                amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt);
                 Ex(i, j, k) = alpha * Ex(i, j, k)
                             + beta * ( - T_Algo::DownwardDz(Hy, coefs_z, n_coefs_z, i, j, k,0)
                                        + T_Algo::DownwardDy(Hz, coefs_y, n_coefs_y, i, j, k,0)
@@ -164,13 +174,14 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                amrex::Real x, y, z;
-                WarpXUtilAlgo::getCellCoordinates (i, j, k, Ey_stag, problo, dx,
-                                                   x, y, z );
-                amrex::Real const sigma = getSigma(x, y, z);
-                amrex::Real const epsilon = getEpsilon(x, y, z);
-                amrex::Real alpha = T_MacroAlgo::alpha( sigma, epsilon, dt);
-                amrex::Real beta = T_MacroAlgo::beta( sigma, epsilon, dt);
+                // 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);
+                // Interpolated permittivity, epsilon, to Ey position on the grid
+                amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag,
+                                           Ey_stag, macro_cr, i, j, k, scomp);
+                amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt);
+                amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt);
 
                 Ey(i, j, k) = alpha * Ey(i, j, k)
                             + beta * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0)
@@ -179,13 +190,14 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                amrex::Real x, y, z;
-                WarpXUtilAlgo::getCellCoordinates (i, j, k, Ez_stag, problo, dx,
-                                                   x, y, z );
-                amrex::Real const sigma = getSigma(x, y, z);
-                amrex::Real const epsilon = getEpsilon(x, y, z);
-                amrex::Real alpha = T_MacroAlgo::alpha( sigma, epsilon, dt);
-                amrex::Real beta = T_MacroAlgo::beta( sigma, epsilon, dt);
+                // 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);
+                // Interpolated permittivity, epsilon, to Ez position on the grid
+                amrex::Real const epsilon_interp = CoarsenIO::Interp( eps_arr, epsilon_stag,
+                                           Ez_stag, macro_cr, i, j, k, scomp);
+                amrex::Real alpha = T_MacroAlgo::alpha( sigma_interp, epsilon_interp, dt);
+                amrex::Real beta = T_MacroAlgo::beta( sigma_interp, epsilon_interp, dt);
 
                 Ez(i, j, k) = alpha * Ez(i, j, k)
                             + beta * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0)