diff options
Diffstat (limited to 'Source/FieldSolver')
8 files changed, 461 insertions, 85 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H index ef6d53416..b688110b7 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H @@ -143,9 +143,10 @@ struct CartesianCKCAlgorithm { /** * Perform derivative along x on a nodal grid, from a cell-centered field `F` */ + template< typename T_Field> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDx ( - amrex::Array4<amrex::Real> const& F, + T_Field const& F, amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { @@ -189,15 +190,16 @@ struct CartesianCKCAlgorithm { /** * Perform derivative along y on a nodal grid, from a cell-centered field `F` */ + template< typename T_Field> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDy ( - amrex::Array4<amrex::Real> const& F, + T_Field const& F, amrex::Real const * const coefs_y, int const n_coefs_y, int const i, int const j, int const k, int const ncomp=0 ) { using namespace amrex; #if defined WARPX_DIM_3D - Real const inv_dy = coefs_y[0]; + amrex::Real const inv_dy = coefs_y[0]; return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); amrex::ignore_unused(n_coefs_y); #elif (defined WARPX_DIM_XZ) @@ -248,11 +250,12 @@ struct CartesianCKCAlgorithm { /** * Perform derivative along z on a nodal grid, from a cell-centered field `F` */ + template< typename T_Field> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDz ( - amrex::Array4<amrex::Real> const& F, + T_Field const& F, amrex::Real const * const coefs_z, int const /*n_coefs_z*/, - int const i, int const j, int const k, int const ncomp=0 ) { + int const i, int const j, int const k, int const ncomp=0) { amrex::Real const inv_dz = coefs_z[0]; #if defined WARPX_DIM_3D diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H index 305fb3507..2ec6a3d22 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -9,6 +9,7 @@ #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_ #include "Utils/WarpXConst.H" +#include "FieldAccessorFunctors.H" #include <AMReX.H> #include <AMReX_REAL.H> @@ -68,9 +69,10 @@ struct CartesianYeeAlgorithm { /** * Perform derivative along x on a nodal grid, from a cell-centered field `F`*/ + template< typename T_Field> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDx ( - amrex::Array4<amrex::Real> const& F, + T_Field const& F, amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { @@ -103,9 +105,10 @@ struct CartesianYeeAlgorithm { /** * Perform derivative along y on a nodal grid, from a cell-centered field `F`*/ + template< typename T_Field> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDy ( - amrex::Array4<amrex::Real> const& F, + T_Field const& F, amrex::Real const * const coefs_y, int const n_coefs_y, int const i, int const j, int const k, int const ncomp=0 ) { @@ -143,9 +146,10 @@ struct CartesianYeeAlgorithm { /** * Perform derivative along z on a nodal grid, from a cell-centered field `F`*/ + template< typename T_Field> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE static amrex::Real DownwardDz ( - amrex::Array4<amrex::Real> const& F, + T_Field const& F, amrex::Real const * const coefs_z, int const /*n_coefs_z*/, int const i, int const j, int const k, int const ncomp=0 ) { diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H new file mode 100644 index 000000000..3ba6de335 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/FieldAccessorFunctors.H @@ -0,0 +1,40 @@ +#ifndef WARPX_FIELD_ACCESSOR_FUNCTORS_H +#define WARPX_FIELD_ACCESSOR_FUNCTORS_H +#include "WarpX.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). + */ +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) {} + + /** + * \brief return field value at (i,j,k,ncomp) scaled by (1/m_parameters(i,j,k,ncomp)) + * \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. + * + * \return m_field/m_paramter 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) ) ; + } +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; +}; + + +#endif diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index 047a8cb98..0c37e3fa8 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -164,7 +164,7 @@ class FiniteDifferenceSolver const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); - template< typename T_Algo > + template< typename T_Algo, typename T_MacroAlgo > void MacroscopicEvolveECartesian ( std::array< std::unique_ptr< amrex::MultiFab>, 3>& Efield, std::array< std::unique_ptr< amrex::MultiFab>, 3> const &Bfield, 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 ¯oscopic_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 diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H index 673eae76e..17ab616ae 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -2,11 +2,14 @@ #define WARPX_MACROSCOPICPROPERTIES_H_ +#include "Parser/WarpXParserWrapper.H" +#include "Utils/WarpXConst.H" #include <AMReX_MultiFab.H> + /** - * \brief This class contains the macroscopic parameters of the medium needed to + * \brief This class contains the macroscopic properties of the medium needed to * evaluate macroscopic Maxwell equation. */ class @@ -14,22 +17,130 @@ MacroscopicProperties { public: MacroscopicProperties (); // constructor - /** \brief Read user-defined macroscopic properties. Called in constructor. */ + /** Read user-defined macroscopic properties. Called in constructor. */ void ReadParameters (); - /** return Real, sigma (conductivity) of the medium. */ - amrex::Real sigma () const noexcept {return m_sigma;} - /** return Real, epsilon (permittivity) of the medium. */ - amrex::Real epsilon () const noexcept {return m_epsilon;} - /** return Real, mu (permeability) of the medium. */ - amrex::Real mu () const noexcept {return m_mu;} + /** 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, + ParserWrapper<3> *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; private: /** Conductivity, sigma, of the medium */ - amrex::Real m_sigma; + amrex::Real m_sigma = 0.0; /** Permittivity, epsilon, of the medium */ - amrex::Real m_epsilon; + amrex::Real m_epsilon = PhysConst::ep0; /** Permeability, mu, of the medium */ - amrex::Real m_mu; + 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 */ + std::string m_epsilon_s = "constant"; + /** Stores initialization type for permeability : constant or parser */ + std::string m_mu_s = "constant"; + + /** string for storing parser function */ + std::string m_str_sigma_function; + std::string m_str_epsilon_function; + std::string m_str_mu_function; + /** Parser Wrappers */ + std::unique_ptr<ParserWrapper<3> > m_sigma_parser; + std::unique_ptr<ParserWrapper<3> > m_epsilon_parser; + std::unique_ptr<ParserWrapper<3> > m_mu_parser; +}; + +/** + * \brief + * This struct contains only static functions to compute the co-efficients for the + * Lax-Wendroff scheme of macroscopic Maxwells equations using + * macroscopic properties, namely, conductivity (sigma), permittivity (epsilon). + * Permeability of the material, mu, is used as (beta/mu) for the E-update + * defined in MacroscopicEvolveECartesian(). + */ +struct LaxWendroffAlgo { + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real alpha (amrex::Real const sigma, + amrex::Real const epsilon, + amrex::Real dt) { + using namespace amrex; + amrex::Real fac1 = 0.5_rt * sigma * dt / epsilon; + amrex::Real alpha = (1._rt - fac1)/(1._rt + fac1); + return alpha; + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real beta (amrex::Real const sigma, + amrex::Real const epsilon, + amrex::Real dt) { + using namespace amrex; + amrex::Real fac1 = 0.5_rt * sigma * dt / epsilon; + amrex::Real beta = dt / ( epsilon * (1._rt + fac1) ); + return beta; + }; + +}; + +/** + * \brief + * This struct contains only static functions to compute the co-efficients for the + * BackwardEuler scheme of macroscopic Maxwells equations using + * macroscopic properties, namely, conductivity (sigma) and permittivity (epsilon). + * Permeability of the material, mu, is used as (beta/mu) for the E-update + * defined in MacroscopicEvolveECartesian(). + */ +struct BackwardEulerAlgo { + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real alpha (amrex::Real const sigma, + amrex::Real const epsilon, + amrex::Real dt) { + using namespace amrex; + amrex::Real fac1 = sigma * dt / epsilon; + amrex::Real alpha = (1._rt)/(1._rt + fac1); + return alpha; + }; + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real beta (amrex::Real const sigma, + amrex::Real const epsilon, + amrex::Real dt) { + using namespace amrex; + amrex::Real fac1 = sigma * dt / epsilon; + amrex::Real beta = dt / ( epsilon * (1._rt + fac1) ); + return beta; + }; + }; #endif // WARPX_MACROSCOPIC_PROPERTIES_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp index 7bb1911fd..d2f72ad6c 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -1,5 +1,6 @@ #include "MacroscopicProperties.H" #include <AMReX_ParmParse.H> +#include "WarpX.H" using namespace amrex; @@ -12,10 +13,185 @@ void MacroscopicProperties::ReadParameters () { ParmParse pp("macroscopic"); - // Since macroscopic maxwell solve is turned on, user must define sigma, mu, and epsilon // - pp.get("sigma", m_sigma); - pp.get("mu", m_mu); - pp.get("epsilon", m_epsilon); + // Since macroscopic maxwell solve is turned on, + // user-defined sigma, mu, and epsilon are queried. + // The vacuum values are used as default for the macroscopic parameters + // with a warning message to the user to indicate that no value was specified. + + // Query input for material conductivity, sigma. + bool sigma_specified = false; + if (pp.query("sigma", m_sigma)) { + m_sigma_s = "constant"; + sigma_specified = true; + } + if (pp.query("sigma_function(x,y,z)", m_str_sigma_function) ) { + m_sigma_s = "parse_sigma_function"; + sigma_specified = true; + } + if (!sigma_specified) { + amrex::Print() << "WARNING: Material conductivity is not specified. Using default vacuum value of " << m_sigma << " in the simulation\n"; + } + // initialization of sigma (conductivity) with parser + if (m_sigma_s == "parse_sigma_function") { + Store_parserString(pp, "sigma_function(x,y,z)", m_str_sigma_function); + m_sigma_parser.reset(new ParserWrapper<3>( + makeParser(m_str_sigma_function,{"x","y","z"}) ) ); + } + + bool epsilon_specified = false; + if (pp.query("epsilon", m_epsilon)) { + m_epsilon_s = "constant"; + epsilon_specified = true; + } + if (pp.query("epsilon_function(x,y,z)", m_str_epsilon_function) ) { + m_epsilon_s = "parse_epsilon_function"; + epsilon_specified = true; + } + if (!epsilon_specified) { + amrex::Print() << "WARNING: Material permittivity is not specified. Using default vacuum value of " << m_epsilon << " in the simulation\n"; + } + + // initialization of epsilon (permittivity) with parser + if (m_epsilon_s == "parse_epsilon_function") { + Store_parserString(pp, "epsilon_function(x,y,z)", m_str_epsilon_function); + m_epsilon_parser.reset(new ParserWrapper<3>( + makeParser(m_str_epsilon_function,{"x","y","z"}) ) ); + } + + // Query input for material permittivity, epsilon. + bool mu_specified = false; + if (pp.query("mu", m_mu)) { + m_mu_s = "constant"; + mu_specified = true; + } + if (pp.query("mu_function(x,y,z)", m_str_mu_function) ) { + m_mu_s = "parse_mu_function"; + mu_specified = true; + } + if (!mu_specified) { + amrex::Print() << "WARNING: Material permittivity is not specified. Using default vacuum value of " << m_mu << " in the simulation\n"; + } + + // initialization of mu (permeability) with parser + if (m_mu_s == "parse_mu_function") { + Store_parserString(pp, "mu_function(x,y,z)", m_str_mu_function); + m_mu_parser.reset(new ParserWrapper<3>( + makeParser(m_str_mu_function,{"x","y","z"}) ) ); + } + +} + +void +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); + int ng = 3; + // Define material property multifabs using ba and dmap from WarpX instance + // sigma is cell-centered MultiFab + m_sigma_mf = std::make_unique<MultiFab>(amrex::convert(ba,IntVect::TheUnitVector()), dmap, 1, ng); + // epsilon is cell-centered MultiFab + m_eps_mf = std::make_unique<MultiFab>(amrex::convert(ba,IntVect::TheUnitVector()), dmap, 1, ng); + // mu is cell-centered MultiFab + m_mu_mf = std::make_unique<MultiFab>(amrex::convert(ba,IntVect::TheUnitVector()), 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.get(), 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.get(), 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.get(), 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(); + + 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; + } +#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; +#endif + + +} + +void +MacroscopicProperties::InitializeMacroMultiFabUsingParser ( + MultiFab *macro_mf, ParserWrapper<3> *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(); + IntVect grown_iv = iv ; + for ( MFIter mfi(*macro_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + // Initialize ghost cells in addition to valid cells + + const Box& tb = mfi.growntilebox(grown_iv); + 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; + + 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; + + // initialize the macroparameter + macro_fab(i,j,k) = (*macro_parser)(x,y,z); + }); + + } + } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6c14757ea..7043ad809 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -311,8 +311,20 @@ WarpX::MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { else { amrex::Abort("Macroscopic EvolveE is not implemented for lev > 0, yet."); } - if (do_pml) { - amrex::Abort("Macroscopic EvolveE is not implemented for pml boundary condition, yet"); + if (do_pml && pml[lev]->ok()) { + if (patch_type == PatchType::fine) { + m_fdtd_solver_fp[lev]->EvolveEPML( + pml[lev]->GetE_fp(), pml[lev]->GetB_fp(), + pml[lev]->Getj_fp(), pml[lev]->GetF_fp(), + pml[lev]->GetMultiSigmaBox_fp(), + a_dt, pml_has_particles ); + } else { + m_fdtd_solver_cp[lev]->EvolveEPML( + pml[lev]->GetE_cp(), pml[lev]->GetB_cp(), + pml[lev]->Getj_cp(), pml[lev]->GetF_cp(), + pml[lev]->GetMultiSigmaBox_cp(), + a_dt, pml_has_particles ); + } } } |