diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
6 files changed, 241 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index e60ceed3e..31f8008bf 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -9,6 +9,7 @@ #define WARPX_FINITE_DIFFERENCE_SOLVER_H_ #include <AMReX_MultiFab.H> +#include "MacroscopicProperties/MacroscopicProperties.H" #include "BoundaryConditions/PML.H" /** @@ -54,6 +55,24 @@ class FiniteDifferenceSolver void ComputeDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ); + /** + * \brief Macroscopic E-update for non-vacuum medium using the user-selected + * finite-difference algorithm and macroscopic sigma-method defined in + * WarpXAlgorithmSelection.H + * + * \param[out] Efield vector of electric field MultiFabs updated at a given level + * \param[in] Bfield vector of magnetic field MultiFabs at a given level + * \param[in] Jfield vector of current density MultiFabs at a given level + * \param[in] dt timestep of the simulation + * \param[in] macroscopic_properties contains user-defined properties of the medium. + */ + + 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); + void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield, std::array< amrex::MultiFab*, 3 > const Efield, amrex::Real const dt ); @@ -146,6 +165,14 @@ class FiniteDifferenceSolver amrex::MultiFab& divE ); template< typename T_Algo > + void 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); + + template< typename T_Algo > void EvolveBPMLCartesian ( std::array< amrex::MultiFab*, 3 > Bfield, std::array< amrex::MultiFab*, 3 > const Efield, diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp new file mode 100644 index 000000000..b7f2e26da --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp @@ -0,0 +1,152 @@ +#include "Utils/WarpXAlgorithmSelection.H" +#include "FiniteDifferenceSolver.H" +#ifdef WARPX_DIM_RZ + // currently works only for 3D +#else +# include "FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H" +# include "FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H" +#endif +#include "Utils/WarpXConst.H" +#include <AMReX_Gpu.H> +#include <WarpX.H> + +using namespace amrex; + +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 ) { + + // Select algorithm (The choice of algorithm is a runtime option, + // but we compile code for each algorithm, using templates) +#ifdef WARPX_DIM_RZ + amrex::Abort("currently macro E-push does not work for RZ"); +#else + if (m_do_nodal) { + amrex::Abort(" macro E-push does not work for nodal "); + + } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) { + + MacroscopicEvolveECartesian <CartesianYeeAlgorithm> ( 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 ); + } else { + amrex::Abort("Unknown algorithm"); + } +#endif + +} + + +#ifndef WARPX_DIM_RZ + +template<typename T_Algo> +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; + } + + // Loop through the grids, and over the tiles within each grid +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + + // Extract field data for this grid/tile + Array4<Real> const& Ex = Efield[0]->array(mfi); + Array4<Real> const& Ey = Efield[1]->array(mfi); + Array4<Real> const& Ez = Efield[2]->array(mfi); + Array4<Real> const& Bx = Bfield[0]->array(mfi); + Array4<Real> const& By = Bfield[1]->array(mfi); + Array4<Real> const& Bz = Bfield[2]->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(); + Real const * const AMREX_RESTRICT coefs_y = m_stencil_coefs_y.dataPtr(); + int const n_coefs_y = m_stencil_coefs_y.size(); + Real const * const AMREX_RESTRICT coefs_z = m_stencil_coefs_z.dataPtr(); + int const n_coefs_z = m_stencil_coefs_z.size(); + + // 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()); + + + // 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)); + }, + + [=] 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_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)); + } + + ); + + // 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 new file mode 100644 index 000000000..673eae76e --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H @@ -0,0 +1,35 @@ +#ifndef WARPX_MACROSCOPICPROPERTIES_H_ +#define WARPX_MACROSCOPICPROPERTIES_H_ + + +#include <AMReX_MultiFab.H> + + +/** + * \brief This class contains the macroscopic parameters of the medium needed to + * evaluate macroscopic Maxwell equation. + */ +class +MacroscopicProperties +{ +public: + MacroscopicProperties (); // constructor + /** \brief 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;} + +private: + /** Conductivity, sigma, of the medium */ + amrex::Real m_sigma; + /** Permittivity, epsilon, of the medium */ + amrex::Real m_epsilon; + /** Permeability, mu, of the medium */ + amrex::Real m_mu; +}; + +#endif // WARPX_MACROSCOPIC_PROPERTIES_H_ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp new file mode 100644 index 000000000..7bb1911fd --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.cpp @@ -0,0 +1,21 @@ +#include "MacroscopicProperties.H" +#include <AMReX_ParmParse.H> + +using namespace amrex; + +MacroscopicProperties::MacroscopicProperties () +{ + ReadParameters(); +} + +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); + +} + diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/Make.package new file mode 100644 index 000000000..9548d24b8 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/Make.package @@ -0,0 +1,3 @@ +CEXE_sources += MacroscopicProperties.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package index b267463bd..31bccd1e4 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package +++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package @@ -3,8 +3,11 @@ CEXE_sources += EvolveB.cpp CEXE_sources += EvolveE.cpp CEXE_sources += EvolveF.cpp CEXE_sources += ComputeDivE.cpp +CEXE_sources += MacroscopicEvolveE.cpp CEXE_sources += EvolveBPML.cpp CEXE_sources += EvolveEPML.cpp CEXE_sources += EvolveFPML.cpp +include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/Make.package + VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver |