Interpolate mu to B-location when computing H (#2252)

* interpolate mu to B-location when computing H

* add comments

* modify comment to reflect interpolation

* eol

* add struct to compute macroparameters at x,y,z, and util function to compute cellcoordinates

* fix eol

* remove unused variables, multifabs, and initialization of macro multifabs

* host device function

* fix typo

* doxygen comments

5 files changed, 229 insertions, 180 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H
index 0b5245228..c6cb7f9db 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H
@@ -1,41 +1,66 @@
+/* Copyright 2020 Revathi Jambunathan
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
 #ifndef WARPX_FIELD_ACCESSOR_FUNCTORS_H
 #define WARPX_FIELD_ACCESSOR_FUNCTORS_H
 
 #include "WarpX.H"
+#include "Utils/CoarsenIO.H"
+#include "Utils/WarpXUtil.H"
+#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
+#include <AMReX_Array.H>
 
 /**
  * \brief Functor that returns the division of the source m_field Array4 value
-          by macroparameter, m_parameter value at the respective (i,j,k,ncomp).
+          by macroparameter obtained using functor, m_getParameter,
+          at the respective (i,j,k,ncomp).
  */
+template< typename T_GetMacroparameter>
 struct FieldAccessorMacroscopic
 {
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     FieldAccessorMacroscopic ( amrex::Array4<amrex::Real const> const a_field,
-                               amrex::Array4<amrex::Real const> const a_parameter )
-        : m_field(a_field), m_parameter(a_parameter) {}
+                               T_GetMacroparameter const& a_getParameter,
+                               amrex::GpuArray<int,3> const& a_field_stag,
+                               amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const a_domain_lo,
+                               amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const a_dx )
+        : m_field(a_field), m_getParameter(a_getParameter), m_field_stag(a_field_stag),
+          m_domain_lo(a_domain_lo), m_dx(a_dx) {}
 
     /**
-     * \brief return field value at (i,j,k,ncomp) scaled by (1/m_parameters(i,j,k,ncomp))
+     * \brief return field value at (i,j,k,ncomp) scaled by (1/m_getParameter(x,y,z))
+     *
      * \param[in] i      index along x of the Array4, m_field and m_parameter.
      * \param[in] j      index along y of the Array4, m_field and m_parameter.
      * \param[in] k      index along z of the Array4, m_field and m_parameter.
      * \param[in] ncomp  index along fourth component of the Array4, containing field-data
-                         to be returned after diving with zero-th component
-                         of m_paramter.
+                         to be returned after dividing by the macroparameter.
      *
-     * \return           m_field/m_paramter at (i,j,k,ncomp)
+     * \return           m_field/m_getParameter(x,y,z) at (i,j,k,ncomp)
      */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     amrex::Real operator() (int const i, int const j,
                             int const k, int const ncomp) const noexcept
     {
-        return ( m_field(i, j, k, ncomp)  /  m_parameter(i, j, k, 0) ) ;
+        amrex::Real x, y, z;
+        WarpXUtilAlgo::getCellCoordinates(i, j, k, m_field_stag, m_domain_lo, m_dx, x, y, z);
+        return ( m_field(i, j, k, ncomp)  /  m_getParameter(x, y, z) ) ;
     }
 private:
     /** Array4 of the source field to be scaled and returned by the operator() */
     amrex::Array4<amrex::Real const> const m_field;
-    /** Array4 of the macroscopic parameter used to divide m_field in the operator() */
-    amrex::Array4<amrex::Real const> const m_parameter;
+    /** Functor object to return the macroparameter at a given position on the grid.*/
+    T_GetMacroparameter const m_getParameter;
+    /** Staggering of the field multifab, m_field */
+    amrex::GpuArray<int,3> const m_field_stag;
+    /** Lower physical coordinates of the simulation domain. */
+    amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const m_domain_lo;
+    /** Cell-size array */
+    amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const m_dx;
 };
 
 
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index 2a96ae44d..a32c50ab7 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -100,10 +100,11 @@ 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,
-                            amrex::Real const dt,
-                            std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
+                      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);
 
         void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield,
                       std::array< amrex::MultiFab*, 3 > const Efield,
@@ -236,7 +237,8 @@ class FiniteDifferenceSolver
             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);
+            std::unique_ptr<MacroscopicProperties> const& macroscopic_properties,
+            int const lev);
 
         template< typename T_Algo >
         void EvolveBPMLCartesian (
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
index f2a0cc9c7..2647a1c08 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
@@ -10,6 +10,7 @@
 #include "MacroscopicProperties/MacroscopicProperties.H"
 #include "Utils/CoarsenIO.H"
 #include "Utils/WarpXAlgorithmSelection.H"
+#include "Utils/WarpXUtil.H"
 #include "WarpX.H"
 
 #include <AMReX.H>
@@ -36,12 +37,14 @@ 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,
-    amrex::Real const dt, std::unique_ptr<MacroscopicProperties> const& macroscopic_properties ) {
+    amrex::Real const dt,
+    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties, int const lev)
+{
 
    // 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, dt, macroscopic_properties, lev);
     amrex::Abort("currently macro E-push does not work for RZ");
 #else
     if (m_do_nodal) {
@@ -52,13 +55,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
 
             MacroscopicEvolveECartesian <CartesianYeeAlgorithm, LaxWendroffAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
 
         }
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
 
             MacroscopicEvolveECartesian <CartesianYeeAlgorithm, BackwardEulerAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
 
         }
 
@@ -69,12 +72,12 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::LaxWendroff) {
 
             MacroscopicEvolveECartesian <CartesianCKCAlgorithm, LaxWendroffAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
 
         } else if (WarpX::macroscopic_solver_algo == MacroscopicSolverAlgo::BackwardEuler) {
 
             MacroscopicEvolveECartesian <CartesianCKCAlgorithm, BackwardEulerAlgo>
-                       ( Efield, Bfield, Jfield, dt, macroscopic_properties );
+                       ( Efield, Bfield, Jfield, dt, macroscopic_properties, lev );
 
         }
 
@@ -93,21 +96,23 @@ 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 ) {
+    amrex::Real const dt,
+    std::unique_ptr<MacroscopicProperties> const& macroscopic_properties, int const lev)
+{
 
-    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& 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;
+    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
@@ -126,11 +131,6 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
         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();
         int const n_coefs_x = m_stencil_coefs_x.size();
@@ -139,27 +139,24 @@ 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);
+        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);
 
         // 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){
-                //// 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);
+                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);
                 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)
@@ -167,12 +164,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int 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);
+                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);
 
                 Ey(i, j, k) = alpha * Ey(i, j, k)
                             + beta * ( - T_Algo::DownwardDx(Hz, coefs_x, n_coefs_x, i, j, k,0)
@@ -181,12 +179,13 @@ void FiniteDifferenceSolver::MacroscopicEvolveECartesian (
             },
 
             [=] AMREX_GPU_DEVICE (int i, int j, int 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);
+                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);
 
                 Ez(i, j, k) = alpha * Ez(i, j, k)
                             + beta * ( - T_Algo::DownwardDy(Hx, coefs_y, n_coefs_y, i, j, k,0)
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
index e638be1c5..040fe2a67 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
@@ -1,3 +1,10 @@
+/* Copyright 2020 Revathi Jambunathan
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
 #ifndef WARPX_MACROSCOPICPROPERTIES_H_
 #define WARPX_MACROSCOPICPROPERTIES_H_
 
@@ -15,6 +22,81 @@
 #include <memory>
 #include <string>
 
+
+enum MacroparameterInitType { ConstantValue, ParserFunction};
+
+/**
+ * \brief Functor to return macropameter, either constant value, m_value, or
+ *        spatially varying scalar value computed using the parser function, m_parser.
+ */
+
+struct GetMacroparameter
+{
+    /* Type of initialization for macroparameter, constant or parser function */
+    MacroparameterInitType m_type;
+    /* Constant value of the macroparameter. */
+    amrex::Real m_value;
+    /* Parser funtion of the spatially varying macroparameter*/
+    amrex::ParserExecutor<3> m_parser;
+    /**
+     * \brief Functor call. This method returns the value of the macroparameter,
+     *        or property of the medium needed for the macroscopic Maxwell solver,
+     *        at a given location (x,y,z) in the domain.
+     *
+     * @param[in] x x-coordinate of a given location
+     * @param[in] y y-coordinate of a given location
+     * @param[in] z z-coordinate of a given location
+     * @return value of the macroparameter at (x,y,z).
+     *         m_value if init-type is constant
+     *         m_parser(x,y,z) if init-type is parser function
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real operator () (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept
+    {
+        using namespace amrex::literals;
+        if (m_type == ConstantValue)
+        {
+            return m_value;
+        }
+        else if (m_type == ParserFunction)
+        {
+            return m_parser(x,y,z);
+        }
+        else
+        {
+            amrex::Abort("macroparameter init type not valid.");
+            return 0.;
+        }
+        return 0.;
+    }
+};
+
+/**
+ * \brief Functor for conductivity, sigma, of the medium.
+ */
+struct GetSigmaMacroparameter : GetMacroparameter
+{
+    /** Constructor to store the type of intialization, m_type, and the value or parser function. */
+    GetSigmaMacroparameter () noexcept;
+};
+
+/**
+ * \brief Functor for permeability, mu, of the medium.
+ */
+struct GetMuMacroparameter : GetMacroparameter
+{
+    /** Constructor to store the type of intialization, m_type, and the value or parser function. */
+    GetMuMacroparameter () noexcept;
+};
+
+/**
+ * \brief Functor for permittivity, epsilon, of the medium.
+ */
+struct GetEpsilonMacroparameter : GetMacroparameter
+{
+    /** Constructor to store the type of intialization, m_type, and the value or parser function. */
+    GetEpsilonMacroparameter () noexcept;
+};
 /**
  * \brief This class contains the macroscopic properties of the medium needed to
  * evaluate macroscopic Maxwell equation.
@@ -29,47 +111,25 @@ public:
      /** Initialize multifabs storing macroscopic multifabs */
      void InitData ();
 
-     /** return MultiFab, sigma (conductivity) of the medium. */
-     amrex::MultiFab& getsigma_mf  () {return (*m_sigma_mf);}
-     /** return MultiFab, epsilon (permittivity) of the medium. */
-     amrex::MultiFab& getepsilon_mf  () {return (*m_eps_mf);}
-     /** return MultiFab, mu (permeability) of the medium. */
-     amrex::MultiFab& getmu_mf  () {return (*m_mu_mf);}
-
-     /** Initializes the Multifabs storing macroscopic properties
-      *  with user-defined functions(x,y,z).
-      */
-     void InitializeMacroMultiFabUsingParser (amrex::MultiFab *macro_mf,
-                                              amrex::ParserExecutor<3> const& macro_parser,
-                                              int lev);
-     /** Gpu Vector with index type of the conductivity multifab */
-     amrex::GpuArray<int, 3> sigma_IndexType;
-     /** Gpu Vector with index type of the permittivity multifab */
-     amrex::GpuArray<int, 3> epsilon_IndexType;
-     /** Gpu Vector with index type of the permeability multifab */
-     amrex::GpuArray<int, 3> mu_IndexType;
      /** Gpu Vector with index type of the Ex multifab */
      amrex::GpuArray<int, 3> Ex_IndexType;
      /** Gpu Vector with index type of the Ey multifab */
      amrex::GpuArray<int, 3> Ey_IndexType;
      /** Gpu Vector with index type of the Ez multifab */
      amrex::GpuArray<int, 3> Ez_IndexType;
-     /** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */
-     amrex::GpuArray<int, 3> macro_cr_ratio;
+     /** Gpu Vector with index type of the Bx multifab */
+     amrex::GpuArray<int, 3> Bx_IndexType;
+     /** Gpu Vector with index type of the By multifab */
+     amrex::GpuArray<int, 3> By_IndexType;
+     /** Gpu Vector with index type of the Bz multifab */
+     amrex::GpuArray<int, 3> Bz_IndexType;
 
-private:
      /** Conductivity, sigma, of the medium */
      amrex::Real m_sigma = 0.0;
      /** Permittivity, epsilon, of the medium */
      amrex::Real m_epsilon = PhysConst::ep0;
      /** Permeability, mu, of the medium */
      amrex::Real m_mu = PhysConst::mu0;
-     /** Multifab for m_sigma */
-     std::unique_ptr<amrex::MultiFab> m_sigma_mf;
-     /** Multifab for m_epsilon */
-     std::unique_ptr<amrex::MultiFab> m_eps_mf;
-     /** Multifab for m_mu */
-     std::unique_ptr<amrex::MultiFab> m_mu_mf;
      /** Stores initialization type for conductivity : constant or parser */
      std::string m_sigma_s = "constant";
      /** Stores initialization type for permittivity : constant or parser */
@@ -85,6 +145,7 @@ private:
      std::unique_ptr<amrex::Parser> m_sigma_parser;
      std::unique_ptr<amrex::Parser> m_epsilon_parser;
      std::unique_ptr<amrex::Parser> m_mu_parser;
+
 };
 
 /**
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp
index cfa3479da..7166eb99d 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp
@@ -22,6 +22,49 @@
 
 using namespace amrex;
 
+GetSigmaMacroparameter::GetSigmaMacroparameter () noexcept
+{
+    auto& warpx = WarpX::GetInstance();
+    auto& macroscopic_properties = warpx.GetMacroscopicProperties();
+    if (macroscopic_properties.m_sigma_s == "constant") {
+        m_type = ConstantValue;
+        m_value = macroscopic_properties.m_sigma;
+    }
+    else if (macroscopic_properties.m_sigma_s == "parse_sigma_function") {
+        m_type = ParserFunction;
+        m_parser = macroscopic_properties.m_sigma_parser->compile<3>();
+    }
+}
+
+GetMuMacroparameter::GetMuMacroparameter () noexcept
+{
+    auto& warpx = WarpX::GetInstance();
+    auto& macroscopic_properties = warpx.GetMacroscopicProperties();
+    if (macroscopic_properties.m_mu_s == "constant") {
+        m_type = ConstantValue;
+        m_value = macroscopic_properties.m_mu;
+    }
+    else if (macroscopic_properties.m_mu_s == "parse_mu_function") {
+        m_type = ParserFunction;
+        m_parser = macroscopic_properties.m_mu_parser->compile<3>();
+    }
+}
+
+GetEpsilonMacroparameter::GetEpsilonMacroparameter () noexcept
+{
+    auto& warpx = WarpX::GetInstance();
+    auto& macroscopic_properties = warpx.GetMacroscopicProperties();
+    if (macroscopic_properties.m_epsilon_s == "constant") {
+        m_type = ConstantValue;
+        m_value = macroscopic_properties.m_epsilon;
+    }
+    else if (macroscopic_properties.m_epsilon_s == "parse_epsilon_function") {
+        m_type = ParserFunction;
+        m_parser = macroscopic_properties.m_epsilon_parser->compile<3>();
+    }
+}
+
+
 MacroscopicProperties::MacroscopicProperties ()
 {
     ReadParameters();
@@ -105,112 +148,31 @@ MacroscopicProperties::InitData ()
     amrex::Print() << "we are in init data of macro \n";
     auto & warpx = WarpX::GetInstance();
 
-    // Get BoxArray and DistributionMap of warpx instant.
-    int lev = 0;
-    BoxArray ba = warpx.boxArray(lev);
-    DistributionMapping dmap = warpx.DistributionMap(lev);
-    const amrex::IntVect ng = warpx.getngE();
-    // Define material property multifabs using ba and dmap from WarpX instance
-    // sigma is cell-centered MultiFab
-    m_sigma_mf = std::make_unique<MultiFab>(ba, dmap, 1, ng);
-    // epsilon is cell-centered MultiFab
-    m_eps_mf = std::make_unique<MultiFab>(ba, dmap, 1, ng);
-    // mu is cell-centered MultiFab
-    m_mu_mf = std::make_unique<MultiFab>(ba, dmap, 1, ng);
-    // Initialize sigma
-    if (m_sigma_s == "constant") {
-
-        m_sigma_mf->setVal(m_sigma);
-
-    } else if (m_sigma_s == "parse_sigma_function") {
-
-        InitializeMacroMultiFabUsingParser(m_sigma_mf.get(), m_sigma_parser->compile<3>(), lev);
-    }
-    // Initialize epsilon
-    if (m_epsilon_s == "constant") {
-
-        m_eps_mf->setVal(m_epsilon);
-
-    } else if (m_epsilon_s == "parse_epsilon_function") {
-
-        InitializeMacroMultiFabUsingParser(m_eps_mf.get(), m_epsilon_parser->compile<3>(), lev);
-
-    }
-    // Initialize mu
-    if (m_mu_s == "constant") {
-
-        m_mu_mf->setVal(m_mu);
-
-    } else if (m_mu_s == "parse_mu_function") {
-
-        InitializeMacroMultiFabUsingParser(m_mu_mf.get(), m_mu_parser->compile<3>(), lev);
-
-    }
-
-
-    IntVect sigma_stag = m_sigma_mf->ixType().toIntVect();
-    IntVect epsilon_stag = m_eps_mf->ixType().toIntVect();
-    IntVect mu_stag = m_mu_mf->ixType().toIntVect();
     IntVect Ex_stag = warpx.getEfield_fp(0,0).ixType().toIntVect();
     IntVect Ey_stag = warpx.getEfield_fp(0,1).ixType().toIntVect();
     IntVect Ez_stag = warpx.getEfield_fp(0,2).ixType().toIntVect();
+    IntVect Bx_stag = warpx.getBfield_fp(0,0).ixType().toIntVect();
+    IntVect By_stag = warpx.getBfield_fp(0,1).ixType().toIntVect();
+    IntVect Bz_stag = warpx.getBfield_fp(0,2).ixType().toIntVect();
+
 
     for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-        sigma_IndexType[idim]   = sigma_stag[idim];
-        epsilon_IndexType[idim] = epsilon_stag[idim];
-        mu_IndexType[idim]      = mu_stag[idim];
         Ex_IndexType[idim]      = Ex_stag[idim];
         Ey_IndexType[idim]      = Ey_stag[idim];
         Ez_IndexType[idim]      = Ez_stag[idim];
-        macro_cr_ratio[idim]    = 1;
+        Bx_IndexType[idim]      = Bx_stag[idim];
+        By_IndexType[idim]      = By_stag[idim];
+        Bz_IndexType[idim]      = Bz_stag[idim];
     }
 #if (AMREX_SPACEDIM==2)
-        sigma_IndexType[2]   = 0;
-        epsilon_IndexType[2] = 0;
-        mu_IndexType[2]      = 0;
         Ex_IndexType[2]      = 0;
         Ey_IndexType[2]      = 0;
         Ez_IndexType[2]      = 0;
-        macro_cr_ratio[2]    = 1;
+        Bx_IndexType[2]      = 0;
+        By_IndexType[2]      = 0;
+        Bz_IndexType[2]      = 0;
 #endif
 
 
 }
 
-void
-MacroscopicProperties::InitializeMacroMultiFabUsingParser (
-                       MultiFab *macro_mf, ParserExecutor<3> const& macro_parser,
-                       int lev)
-{
-    auto& warpx = WarpX::GetInstance();
-    const auto dx_lev = warpx.Geom(lev).CellSizeArray();
-    const RealBox& real_box = warpx.Geom(lev).ProbDomain();
-    IntVect iv = macro_mf->ixType().toIntVect();
-    for ( MFIter mfi(*macro_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-        // Initialize ghost cells in addition to valid cells
-
-        const Box& tb = mfi.tilebox(iv, macro_mf->nGrowVect());
-        auto const& macro_fab =  macro_mf->array(mfi);
-        amrex::ParallelFor (tb,
-            [=] AMREX_GPU_DEVICE (int i, int j, int k) {
-                // Shift x, y, z position based on index type
-                Real fac_x = (1._rt - iv[0]) * dx_lev[0] * 0.5_rt;
-                Real x = i * dx_lev[0] + real_box.lo(0) + fac_x;
-#if (AMREX_SPACEDIM==2)
-                amrex::Real y = 0._rt;
-                Real fac_z = (1._rt - iv[1]) * dx_lev[1] * 0.5_rt;
-                Real z = j * dx_lev[1] + real_box.lo(1) + fac_z;
-#else
-                Real fac_y = (1._rt - iv[1]) * dx_lev[1] * 0.5_rt;
-                Real y = j * dx_lev[1] + real_box.lo(1) + fac_y;
-                Real fac_z = (1._rt - iv[2]) * dx_lev[2] * 0.5_rt;
-                Real z = k * dx_lev[2] + real_box.lo(2) + fac_z;
-#endif
-                // initialize the macroparameter
-                macro_fab(i,j,k) = macro_parser(x,y,z);
-        });
-
-    }
-
-
-}