diff options
author | 2021-11-16 10:15:53 -0800 | |
---|---|---|
committer | 2021-11-16 10:15:53 -0800 | |
commit | d4eaa378c40725a5219d58c8df1b4ede37dc9418 (patch) | |
tree | 6430b7546bac6556c99bc4ebc455952694716820 /Source/FieldSolver/FiniteDifferenceSolver | |
parent | 13f00cab1465d6e8ed22b2417421c0f76cbb7d64 (diff) | |
download | WarpX-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')
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); + }); + } } |