aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver/FiniteDifferenceSolver
diff options
context:
space:
mode:
authorGravatar Revathi Jambunathan <41089244+RevathiJambunathan@users.noreply.github.com> 2021-11-16 10:15:53 -0800
committerGravatar GitHub <noreply@github.com> 2021-11-16 10:15:53 -0800
commitd4eaa378c40725a5219d58c8df1b4ede37dc9418 (patch)
tree6430b7546bac6556c99bc4ebc455952694716820 /Source/FieldSolver/FiniteDifferenceSolver
parent13f00cab1465d6e8ed22b2417421c0f76cbb7d64 (diff)
downloadWarpX-d4eaa378c40725a5219d58c8df1b4ede37dc9418.tar.gz
WarpX-d4eaa378c40725a5219d58c8df1b4ede37dc9418.tar.zst
WarpX-d4eaa378c40725a5219d58c8df1b4ede37dc9418.zip
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>
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H28
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H6
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp94
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H111
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp152
5 files changed, 187 insertions, 204 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H
index c6cb7f9db..8b39296f9 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H
@@ -16,20 +16,14 @@
/**
* \brief Functor that returns the division of the source m_field Array4 value
- by macroparameter obtained using functor, m_getParameter,
- at the respective (i,j,k,ncomp).
+ by macroparameter obtained using m_parameter, at the respective (i,j,k).
*/
-template< typename T_GetMacroparameter>
struct FieldAccessorMacroscopic
{
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
FieldAccessorMacroscopic ( amrex::Array4<amrex::Real const> const a_field,
- 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) {}
+ amrex::Array4<amrex::Real> const& a_parameter)
+ : m_field(a_field), m_parameter(a_parameter) {}
/**
* \brief return field value at (i,j,k,ncomp) scaled by (1/m_getParameter(x,y,z))
@@ -38,7 +32,7 @@ struct FieldAccessorMacroscopic
* \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 dividing by the macroparameter.
+ * to be returned after dividing by the macroparameter.
*
* \return m_field/m_getParameter(x,y,z) at (i,j,k,ncomp)
*/
@@ -46,21 +40,13 @@ struct FieldAccessorMacroscopic
amrex::Real operator() (int const i, int const j,
int const k, int const ncomp) const noexcept
{
- 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) ) ;
+ return ( m_field(i, j, k, ncomp) / m_parameter(i, j, k) ) ;
}
private:
/** Array4 of the source field to be scaled and returned by the operator() */
amrex::Array4<amrex::Real const> const m_field;
- /** 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;
+ /** Array4 of the macroscopic parameter used to divide m_field in the operator() */
+ amrex::Array4<amrex::Real const> const m_parameter;
};
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index a32c50ab7..ab9b31e4d 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -103,8 +103,7 @@ 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,
- int const lev);
+ std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > const Efield,
@@ -237,8 +236,7 @@ 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,
- int const lev);
+ std::unique_ptr<MacroscopicProperties> const& macroscopic_properties);
template< typename T_Algo >
void EvolveBPMLCartesian (
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)
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
index 040fe2a67..b8e1db2a6 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H
@@ -23,80 +23,6 @@
#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.
@@ -111,18 +37,36 @@ 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,
+ const 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 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;
+ /** Gpu Vector with index type of coarsening ratio with default value (1,1,1) */
+ amrex::GpuArray<int, 3> macro_cr_ratio;
+
+private:
/** Conductivity, sigma, of the medium */
amrex::Real m_sigma = 0.0;
@@ -130,6 +74,13 @@ public:
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 */
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp
index e8941b6e3..3cd607392 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp
@@ -15,6 +15,7 @@
#include <AMReX_ParmParse.H>
#include <AMReX_Print.H>
#include <AMReX_RealBox.H>
+#include <AMReX_Parser.H>
#include <AMReX_BaseFwd.H>
@@ -23,49 +24,6 @@
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();
@@ -100,7 +58,7 @@ MacroscopicProperties::ReadParameters ()
// initialization of sigma (conductivity) with parser
if (m_sigma_s == "parse_sigma_function") {
Store_parserString(pp_macroscopic, "sigma_function(x,y,z)", m_str_sigma_function);
- m_sigma_parser = std::make_unique<Parser>(
+ m_sigma_parser = std::make_unique<amrex::Parser>(
makeParser(m_str_sigma_function,{"x","y","z"}));
}
@@ -124,7 +82,7 @@ MacroscopicProperties::ReadParameters ()
// initialization of epsilon (permittivity) with parser
if (m_epsilon_s == "parse_epsilon_function") {
Store_parserString(pp_macroscopic, "epsilon_function(x,y,z)", m_str_epsilon_function);
- m_epsilon_parser = std::make_unique<Parser>(
+ m_epsilon_parser = std::make_unique<amrex::Parser>(
makeParser(m_str_epsilon_function,{"x","y","z"}));
}
@@ -149,7 +107,7 @@ MacroscopicProperties::ReadParameters ()
// initialization of mu (permeability) with parser
if (m_mu_s == "parse_mu_function") {
Store_parserString(pp_macroscopic, "mu_function(x,y,z)", m_str_mu_function);
- m_mu_parser = std::make_unique<Parser>(
+ m_mu_parser = std::make_unique<amrex::Parser>(
makeParser(m_str_mu_function,{"x","y","z"}));
}
@@ -161,30 +119,108 @@ MacroscopicProperties::InitData ()
amrex::Print() << "we are in init data of macro \n";
auto & warpx = WarpX::GetInstance();
- 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();
+ // Get BoxArray and DistributionMap of warpx instance.
+ int lev = 0;
+ amrex::BoxArray ba = warpx.boxArray(lev);
+ amrex::DistributionMapping dmap = warpx.DistributionMap(lev);
+ const amrex::IntVect ng_EB_alloc = warpx.getngE();
+ // Define material property multifabs using ba and dmap from WarpX instance
+ // sigma is cell-centered MultiFab
+ m_sigma_mf = std::make_unique<amrex::MultiFab>(ba, dmap, 1, ng_EB_alloc);
+ // epsilon is cell-centered MultiFab
+ m_eps_mf = std::make_unique<amrex::MultiFab>(ba, dmap, 1, ng_EB_alloc);
+ // mu is cell-centered MultiFab
+ m_mu_mf = std::make_unique<amrex::MultiFab>(ba, dmap, 1, ng_EB_alloc);
+ // 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);
+
+ }
+
+ amrex::IntVect sigma_stag = m_sigma_mf->ixType().toIntVect();
+ amrex::IntVect epsilon_stag = m_eps_mf->ixType().toIntVect();
+ amrex::IntVect mu_stag = m_mu_mf->ixType().toIntVect();
+ amrex::IntVect Ex_stag = warpx.getEfield_fp(0,0).ixType().toIntVect();
+ amrex::IntVect Ey_stag = warpx.getEfield_fp(0,1).ixType().toIntVect();
+ amrex::IntVect Ez_stag = warpx.getEfield_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];
- Bx_IndexType[idim] = Bx_stag[idim];
- By_IndexType[idim] = By_stag[idim];
- Bz_IndexType[idim] = Bz_stag[idim];
+ macro_cr_ratio[idim] = 1;
}
#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;
- Bx_IndexType[2] = 0;
- By_IndexType[2] = 0;
- Bz_IndexType[2] = 0;
+ macro_cr_ratio[2] = 1;
#endif
+}
+void
+MacroscopicProperties::InitializeMacroMultiFabUsingParser (
+ amrex::MultiFab *macro_mf,
+ amrex::ParserExecutor<3> const& macro_parser,
+ const int lev)
+{
+ WarpX& warpx = WarpX::GetInstance();
+ const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> dx_lev = warpx.Geom(lev).CellSizeArray();
+ const amrex::RealBox& real_box = warpx.Geom(lev).ProbDomain();
+ amrex::IntVect iv = macro_mf->ixType().toIntVect();
+ for ( amrex::MFIter mfi(*macro_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+ // Initialize ghost cells in addition to valid cells
+
+ const amrex::Box& tb = mfi.tilebox( iv, macro_mf->nGrowVect());
+ amrex::Array4<amrex::Real> 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
+ amrex::Real fac_x = (1._rt - iv[0]) * dx_lev[0] * 0.5_rt;
+ amrex::Real x = i * dx_lev[0] + real_box.lo(0) + fac_x;
+#if (AMREX_SPACEDIM==2)
+ amrex::Real y = 0._rt;
+ amrex::Real fac_z = (1._rt - iv[1]) * dx_lev[1] * 0.5_rt;
+ amrex::Real z = j * dx_lev[1] + real_box.lo(1) + fac_z;
+#else
+ amrex::Real fac_y = (1._rt - iv[1]) * dx_lev[1] * 0.5_rt;
+ amrex::Real y = j * dx_lev[1] + real_box.lo(1) + fac_y;
+ amrex::Real fac_z = (1._rt - iv[2]) * dx_lev[2] * 0.5_rt;
+ amrex::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);
+ });
+ }
}