Variation of macroscopic properties for E-update (#1016)

* Adding macro-E Push and new file

* Add macroEvolveE, call it, and include algo selection in utils

* fix eol

* Fixing bug in macroE for sigma method 1

* changing MacroEvolveE to MacroscopicEvolveE

* add class for macroscopicproperties and an object in WarpX class

* fix eol

* adding templated ckc call with comment that EvolveE is same for yee and ckc

* add header file pointing to ckc algorithm

* adding obejct m_macroscopic_properties to access sigma,eps,mu

* some cleaning

* Adding comments

* adding documentation

* spelling wandroff to wendroff

* fixing eol

* eol

* const in the right place. Thanks bot!

* profiler for macroscopic evolveE

* MultiFab macroproperties with constant init, templated macroEvolveE,

* call macroparameter init

* eol fix

* add parser for macroscopic properties

* fix eol

* adding input file

* __device__ lambda cannot be private in class

* fix grown tilebox declaration for init data

* [skip ci] some more merge conflicts

* some comments

* removing redundant calls to Macroscopic EvolveE

* [skip ci] fix growntilebox for initializing macro mf

* clean and fix BackwardEuler call for ckc

* commenting out old alpha and beta implementations

* temporarily commiting local changes with calls in Interp.

* fixing a typo

* clean and add documentation

* remove the test input file

* fix typo

* eol fix

* Update Docs/source/running_cpp/parameters.rst

Co-authored-by: Andy Nonaka <AJNonaka@lbl.gov>

* PR suggestions

* Update Source/FieldSolver/WarpXPushFieldsEM.cpp

Co-authored-by: Andy Nonaka <AJNonaka@lbl.gov>

* removing unnecessary includes

* adding 2D initialization for stag arrays

* removing init_style input parameter for material properties

* eol fix

* Adding dE/dt eq curl of (B/mu)

* PhysConst ep0 and mu0

* Add functor for field access for macro and vacuum (B/mu)

* fix eol

* Update Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H

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

* Apply suggestions from code review

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

* Update Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H

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

* Apply suggestions from code review

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

* fixing compilation errors and removing Gpu ManagedVector

Co-authored-by: Andy Nonaka <AJNonaka@lbl.gov>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

1 files changed, 89 insertions, 59 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
index b7f2e26da..7212b6276 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
@@ -5,8 +5,10 @@
 #else
 #   include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
 #   include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+#   include "FiniteDifferenceAlgorithms/FieldAccessorFunctors.H"
 #endif
 #include "Utils/WarpXConst.H"
+#include "Utils/CoarsenIO.H"
 #include <AMReX_Gpu.H>
 #include <WarpX.H>
 
@@ -28,15 +30,35 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        MacroscopicEvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, dt,
-                                                              macroscopic_properties );
+        if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
+
+            MacroscopicEvolveECartesian <CartesianYeeAlgorithm, LaxWendroffAlgo>
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+
+        }
+        if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
+
+            MacroscopicEvolveECartesian <CartesianYeeAlgorithm, BackwardEulerAlgo>
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+
+        }
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
         // Note : EvolveE is the same for CKC and Yee.
         // In the templated Yee and CKC calls, the core operations for EvolveE is the same.
-        MacroscopicEvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, dt,
-                                                              macroscopic_properties );
+        if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
+
+            MacroscopicEvolveECartesian <CartesianCKCAlgorithm, LaxWendroffAlgo>
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+
+        } else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
+
+            MacroscopicEvolveECartesian <CartesianCKCAlgorithm, BackwardEulerAlgo>
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+
+        }
+
     } else {
         amrex::Abort("Unknown algorithm");
     }
@@ -47,34 +69,27 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
 
 #ifndef WARPX_DIM_RZ
 
-template<typename T_Algo>
+template<typename T_Algo, typename T_MacroAlgo>
 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,
     amrex::Real const dt, std::unique_ptr<MacroscopicProperties> const& macroscopic_properties ) {
 
-    const int &macroscopic_solver_algo = WarpX::macroscopic_solver_algo;
-    Real sigma = macroscopic_properties->sigma();
-    Real const mu = macroscopic_properties->mu();
-    Real const epsilon = macroscopic_properties->epsilon();
-
-    Real alpha = 0._rt;
-    Real beta  = 0._rt;
-    Real fac1 = 0._rt;
-    Real inv_fac = 0._rt;
-
-    if (macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
-        fac1 = sigma * dt / epsilon;
-        inv_fac = 1._rt / ( 1._rt + fac1);
-        alpha = inv_fac;
-        beta  = dt * inv_fac / epsilon;
-    } else if (macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
-        fac1 = 0.5_rt * sigma * dt / epsilon;
-        inv_fac = 1._rt / ( 1._rt + fac1);
-        alpha = (1.0_rt - fac1) * inv_fac;
-        beta  = dt * inv_fac / epsilon;
-    }
+    auto& sigma_mf = macroscopic_properties->getsigma_mf();
+    auto& epsilon_mf = macroscopic_properties->getepsilon_mf();
+    auto& 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& mu_stag = macroscopic_properties->mu_IndexType;
+    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& macro_cr     = macroscopic_properties->macro_cr_ratio;
+
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef _OPENMP
@@ -89,6 +104,14 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
         Array4<Real> const& Bx = Bfield[0]->array(mfi);
         Array4<Real> const& By = Bfield[1]->array(mfi);
         Array4<Real> const& Bz = Bfield[2]->array(mfi);
+        Array4<Real> const& jx = Jfield[0]->array(mfi);
+        Array4<Real> const& jy = Jfield[1]->array(mfi);
+        Array4<Real> const& jz = Jfield[2]->array(mfi);
+
+        // material prop //
+        Array4<Real> const& sigma_arr = sigma_mf.array(mfi);
+        Array4<Real> const& eps_arr = epsilon_mf.array(mfi);
+        Array4<Real> const& mu_arr = mu_mf.array(mfi);
 
         // Extract stencil coefficients
         Real const * const AMREX_RESTRICT coefs_x = m_stencil_coefs_x.dataPtr();
@@ -98,55 +121,62 @@ 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 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){
-                Ex(i, j, k) = alpha * Ex(i, j, k) + (beta/mu)
-                     * ( - T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k)
-                         + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k));
+                //// 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)
+                                     ) - beta * jx(i, j, k);
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                Ey(i, j, k) = alpha * Ey(i, j, k) + (beta/mu)
-                     * ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k)
-                         + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k));
+                amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag,
+                                           Ey_stag, macro_cr, i, j, k, scomp);
+                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)
+                                       + T_Algo::DownwardDz(Hx, coefs_z, n_coefs_z, i, j, k,0)
+                                     ) - beta * jy(i, j, k);
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int k){
-                Ez(i, j, k) = alpha * Ez(i, j, k) + (beta/mu)
-                     * ( - T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k)
-                         + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k));
+                amrex::Real const sigma_interp = CoarsenIO::Interp( sigma_arr, sigma_stag,
+                                           Ez_stag, macro_cr, i, j, k, scomp);
+                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)
+                                       + T_Algo::DownwardDx(Hy, coefs_x, n_coefs_x, i, j, k,0)
+                                     ) - beta * jz(i, j, k);
             }
-
         );
-
-        // update E using J, if source currents are specified.
-        if (Jfield[0]) {
-            Array4<Real> const& jx = Jfield[0]->array(mfi);
-            Array4<Real> const& jy = Jfield[1]->array(mfi);
-            Array4<Real> const& jz = Jfield[2]->array(mfi);
-
-            amrex::ParallelFor(tex, tey, tez,
-                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                    Ex(i, j, k) += -beta * jx(i, j, k);
-                },
-                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                    Ey(i, j, k) += -beta * jy(i, j, k);
-                },
-                [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                    Ez(i, j, k) += -beta * jz(i, j, k);
-                }
-            );
-        }
     }
-
 }
 
 #endif // corresponds to ifndef WARPX_DIM_RZ