aboutsummaryrefslogtreecommitdiff
path: root/Source/FieldSolver/FiniteDifferenceSolver
diff options
context:
space:
mode:
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver')
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp131
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp206
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp113
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H35
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/Make.package3
5 files changed, 488 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
new file mode 100644
index 000000000..6282892ab
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
@@ -0,0 +1,131 @@
+/* Copyright 2020 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
+#ifdef WARPX_DIM_RZ
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
+#endif
+#include "BoundaryConditions/PMLComponent.H"
+#include <AMReX_Gpu.H>
+
+using namespace amrex;
+
+/**
+ * \brief Update the B field, over one timestep
+ */
+void FiniteDifferenceSolver::EvolveBPML (
+ std::array< amrex::MultiFab*, 3 > Bfield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt ) {
+
+ // 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("PML are not implemented in cylindrical geometry.");
+#else
+ if (m_do_nodal) {
+
+ EvolveBPMLCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, dt );
+
+ } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+
+ EvolveBPMLCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, dt );
+
+ } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
+
+ EvolveBPMLCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, dt );
+
+ } else {
+ amrex::Abort("Unknown algorithm");
+ }
+#endif
+}
+
+
+#ifndef WARPX_DIM_RZ
+
+template<typename T_Algo>
+void FiniteDifferenceSolver::EvolveBPMLCartesian (
+ std::array< amrex::MultiFab*, 3 > Bfield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt ) {
+
+ // Loop through the grids, and over the tiles within each grid
+#ifdef _OPENMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+ for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+
+ // Extract field data for this grid/tile
+ 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& Ex = Efield[0]->array(mfi);
+ Array4<Real> const& Ey = Efield[1]->array(mfi);
+ Array4<Real> const& Ez = Efield[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& tbx = mfi.tilebox(Bfield[0]->ixType().ixType());
+ Box const& tby = mfi.tilebox(Bfield[1]->ixType().ixType());
+ Box const& tbz = mfi.tilebox(Bfield[2]->ixType().ixType());
+
+ // Loop over the cells and update the fields
+ amrex::ParallelFor(tbx, tby, tbz,
+
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ Bx(i, j, k, PMLComp::xz) += dt * (
+ T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx)
+ + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yy)
+ + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) );
+ Bx(i, j, k, PMLComp::xy) -= dt * (
+ T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zx)
+ + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zy)
+ + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) );
+ },
+
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ By(i, j, k, PMLComp::yx) += dt * (
+ T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zx)
+ + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zy)
+ + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) );
+ By(i, j, k, PMLComp::yz) -= dt * (
+ T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xx)
+ + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xy)
+ + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) );
+ },
+
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ Bz(i, j, k, PMLComp::zy) += dt * (
+ T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xx)
+ + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xy)
+ + T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) );
+ Bz(i, j, k, PMLComp::zx) -= dt * (
+ T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yx)
+ + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yy)
+ + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) );
+ }
+
+ );
+
+ }
+
+}
+
+#endif // corresponds to ifndef WARPX_DIM_RZ
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
new file mode 100644
index 000000000..34144ff0d
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
@@ -0,0 +1,206 @@
+/* Copyright 2020 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
+#ifdef WARPX_DIM_RZ
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
+#endif
+#include "BoundaryConditions/PML.H"
+#include "BoundaryConditions/PML_current.H"
+#include "BoundaryConditions/PMLComponent.H"
+#include "Utils/WarpXConst.H"
+#include <AMReX_Gpu.H>
+
+using namespace amrex;
+
+/**
+ * \brief Update the E field, over one timestep
+ */
+void FiniteDifferenceSolver::EvolveEPML (
+ std::array< amrex::MultiFab*, 3 > Efield,
+ std::array< amrex::MultiFab*, 3 > const Bfield,
+ std::array< amrex::MultiFab*, 3 > const Jfield,
+ amrex::MultiFab* const Ffield,
+ MultiSigmaBox const& sigba,
+ amrex::Real const dt, bool pml_has_particles ) {
+
+ // 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("PML are not implemented in cylindrical geometry.");
+#else
+ if (m_do_nodal) {
+
+ EvolveEPMLCartesian <CartesianNodalAlgorithm> (
+ Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
+
+ } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+
+ EvolveEPMLCartesian <CartesianYeeAlgorithm> (
+ Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
+
+ } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
+
+ EvolveEPMLCartesian <CartesianCKCAlgorithm> (
+ Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
+
+ } else {
+ amrex::Abort("Unknown algorithm");
+ }
+#endif
+}
+
+
+#ifndef WARPX_DIM_RZ
+
+template<typename T_Algo>
+void FiniteDifferenceSolver::EvolveEPMLCartesian (
+ std::array< amrex::MultiFab*, 3 > Efield,
+ std::array< amrex::MultiFab*, 3 > const Bfield,
+ std::array< amrex::MultiFab*, 3 > const Jfield,
+ amrex::MultiFab* const Ffield,
+ MultiSigmaBox const& sigba,
+ amrex::Real const dt, bool pml_has_particles ) {
+
+ Real constexpr c2 = PhysConst::c * PhysConst::c;
+
+ // 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().ixType());
+ Box const& tey = mfi.tilebox(Efield[1]->ixType().ixType());
+ Box const& tez = mfi.tilebox(Efield[2]->ixType().ixType());
+
+ // 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, PMLComp::xz) -= c2 * dt * (
+ T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, PMLComp::yx)
+ + T_Algo::DownwardDz(By, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) );
+ Ex(i, j, k, PMLComp::xy) += c2 * dt * (
+ T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k, PMLComp::zx)
+ + T_Algo::DownwardDy(Bz, coefs_y, n_coefs_y, i, j, k, PMLComp::zy) );
+ },
+
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ Ey(i, j, k, PMLComp::yx) -= c2 * dt * (
+ T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, PMLComp::zx)
+ + T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k, PMLComp::zy) );
+ Ey(i, j, k, PMLComp::yz) += c2 * dt * (
+ T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k, PMLComp::xy)
+ + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) );
+ },
+
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ Ez(i, j, k, PMLComp::zy) -= c2 * dt * (
+ T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, PMLComp::xy)
+ + T_Algo::DownwardDy(Bx, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) );
+ Ez(i, j, k, PMLComp::zx) += c2 * dt * (
+ T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k, PMLComp::yx)
+ + T_Algo::DownwardDx(By, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) );
+ }
+
+ );
+
+ // If F is not a null pointer, further update E using the grad(F) term
+ // (hyperbolic correction for errors in charge conservation)
+ if (Ffield) {
+ // Extract field data for this grid/tile
+ Array4<Real> const& F = Ffield->array(mfi);
+
+ // 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, PMLComp::xx) += c2 * dt * (
+ T_Algo::UpwardDx(F, coefs_x, n_coefs_x, i, j, k, PMLComp::x)
+ + T_Algo::UpwardDx(F, coefs_x, n_coefs_x, i, j, k, PMLComp::y)
+ + T_Algo::UpwardDx(F, coefs_x, n_coefs_x, i, j, k, PMLComp::z) );
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ Ey(i, j, k, PMLComp::yy) += c2 * dt * (
+ T_Algo::UpwardDy(F, coefs_y, n_coefs_y, i, j, k, PMLComp::x)
+ + T_Algo::UpwardDy(F, coefs_y, n_coefs_y, i, j, k, PMLComp::y)
+ + T_Algo::UpwardDy(F, coefs_y, n_coefs_y, i, j, k, PMLComp::z) );
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+ Ez(i, j, k, PMLComp::zz) += c2 * dt * (
+ T_Algo::UpwardDz(F, coefs_z, n_coefs_z, i, j, k, PMLComp::x)
+ + T_Algo::UpwardDz(F, coefs_z, n_coefs_z, i, j, k, PMLComp::y)
+ + T_Algo::UpwardDz(F, coefs_z, n_coefs_z, i, j, k, PMLComp::z) );
+ }
+ );
+ }
+
+ // Update the E field in the PML, using the current
+ // deposited by the particles in the PML
+ if (pml_has_particles) {
+
+ // Extract field data for this grid/tile
+ Array4<Real> const& Jx = Jfield[0]->array(mfi);
+ Array4<Real> const& Jy = Jfield[1]->array(mfi);
+ Array4<Real> const& Jz = Jfield[2]->array(mfi);
+ const Real* sigmaj_x = sigba[mfi].sigma[0].data();
+ const Real* sigmaj_y = sigba[mfi].sigma[1].data();
+ const Real* sigmaj_z = sigba[mfi].sigma[2].data();
+ int const x_lo = sigba[mfi].sigma[0].lo();
+#if (AMREX_SPACEDIM == 3)
+ int const y_lo = sigba[mfi].sigma[1].lo();
+ int const z_lo = sigba[mfi].sigma[2].lo();
+#else
+ int const y_lo = 0;
+ int const z_lo = sigba[mfi].sigma[1].lo();
+#endif
+ const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt;
+
+ amrex::ParallelFor( tex, tey, tez,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ push_ex_pml_current(i, j, k, Ex, Jx,
+ sigmaj_y, sigmaj_z, y_lo, z_lo, mu_c2_dt);
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ push_ey_pml_current(i, j, k, Ey, Jy,
+ sigmaj_x, sigmaj_z, x_lo, z_lo, mu_c2_dt);
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ push_ez_pml_current(i, j, k, Ez, Jz,
+ sigmaj_x, sigmaj_y, x_lo, y_lo, mu_c2_dt);
+ }
+ );
+ }
+
+ }
+
+}
+
+#endif // corresponds to ifndef WARPX_DIM_RZ
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp
new file mode 100644
index 000000000..b36c50c35
--- /dev/null
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp
@@ -0,0 +1,113 @@
+/* Copyright 2020 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "Utils/WarpXAlgorithmSelection.H"
+#include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
+#ifdef WARPX_DIM_RZ
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#else
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H"
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H"
+# include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H"
+#endif
+#include "BoundaryConditions/PMLComponent.H"
+#include <AMReX_Gpu.H>
+
+using namespace amrex;
+
+/**
+ * \brief Update the E field, over one timestep
+ */
+void FiniteDifferenceSolver::EvolveFPML (
+ amrex::MultiFab* Ffield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt ) {
+
+ // 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("PML are not implemented in cylindrical geometry.");
+#else
+ if (m_do_nodal) {
+
+ EvolveFPMLCartesian <CartesianNodalAlgorithm> ( Ffield, Efield, dt );
+
+ } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+
+ EvolveFPMLCartesian <CartesianYeeAlgorithm> ( Ffield, Efield, dt );
+
+ } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
+
+ EvolveFPMLCartesian <CartesianCKCAlgorithm> ( Ffield, Efield, dt );
+
+ } else {
+ amrex::Abort("Unknown algorithm");
+ }
+#endif
+}
+
+
+#ifndef WARPX_DIM_RZ
+
+template<typename T_Algo>
+void FiniteDifferenceSolver::EvolveFPMLCartesian (
+ amrex::MultiFab* Ffield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt ) {
+
+ // Loop through the grids, and over the tiles within each grid
+#ifdef _OPENMP
+#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+ for ( MFIter mfi(*Ffield, TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+
+ // Extract field data for this grid/tile
+ Array4<Real> const& F = Ffield->array(mfi);
+ Array4<Real> const& Ex = Efield[0]->array(mfi);
+ Array4<Real> const& Ey = Efield[1]->array(mfi);
+ Array4<Real> const& Ez = Efield[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& tf = mfi.tilebox(Ffield->ixType().ixType());
+
+ // Loop over the cells and update the fields
+ amrex::ParallelFor(tf,
+
+ [=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+ F(i, j, k, PMLComp::x) += dt * (
+ T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k, PMLComp::xx)
+ + T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k, PMLComp::xy)
+ + T_Algo::DownwardDx(Ex, coefs_x, n_coefs_x, i, j, k, PMLComp::xz) );
+
+ F(i, j, k, PMLComp::y) += dt * (
+ T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k, PMLComp::yx)
+ + T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k, PMLComp::yy)
+ + T_Algo::DownwardDy(Ey, coefs_y, n_coefs_y, i, j, k, PMLComp::yz) );
+
+ F(i, j, k, PMLComp::z) += dt * (
+ T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k, PMLComp::zx)
+ + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k, PMLComp::zy)
+ + T_Algo::DownwardDz(Ez, coefs_z, n_coefs_z, i, j, k, PMLComp::zz) );
+
+ }
+
+ );
+
+ }
+
+}
+
+#endif // corresponds to ifndef WARPX_DIM_RZ
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index b23e3bbed..e60ceed3e 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 "BoundaryConditions/PML.H"
/**
* \brief Top-level class for the electromagnetic finite-difference solver
@@ -53,6 +54,20 @@ class FiniteDifferenceSolver
void ComputeDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
amrex::MultiFab& divE );
+ void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt );
+
+ void EvolveEPML ( std::array< amrex::MultiFab*, 3 > Efield,
+ std::array< amrex::MultiFab*, 3 > const Bfield,
+ std::array< amrex::MultiFab*, 3 > const Jfield,
+ amrex::MultiFab* const Ffield,
+ MultiSigmaBox const& sigba,
+ amrex::Real const dt, bool pml_has_particles );
+
+ void EvolveFPML ( amrex::MultiFab* Ffield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt );
private:
@@ -130,6 +145,26 @@ class FiniteDifferenceSolver
const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
amrex::MultiFab& divE );
+ template< typename T_Algo >
+ void EvolveBPMLCartesian (
+ std::array< amrex::MultiFab*, 3 > Bfield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt );
+
+ template< typename T_Algo >
+ void EvolveEPMLCartesian (
+ std::array< amrex::MultiFab*, 3 > Efield,
+ std::array< amrex::MultiFab*, 3 > const Bfield,
+ std::array< amrex::MultiFab*, 3 > const Jfield,
+ amrex::MultiFab* const Ffield,
+ MultiSigmaBox const& sigba,
+ amrex::Real const dt, bool pml_has_particles );
+
+ template< typename T_Algo >
+ void EvolveFPMLCartesian ( amrex::MultiFab* Ffield,
+ std::array< amrex::MultiFab*, 3 > const Efield,
+ amrex::Real const dt );
+
#endif
};
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/Make.package b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
index 9d95e04cc..b267463bd 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/Make.package
+++ b/Source/FieldSolver/FiniteDifferenceSolver/Make.package
@@ -3,5 +3,8 @@ CEXE_sources += EvolveB.cpp
CEXE_sources += EvolveE.cpp
CEXE_sources += EvolveF.cpp
CEXE_sources += ComputeDivE.cpp
+CEXE_sources += EvolveBPML.cpp
+CEXE_sources += EvolveEPML.cpp
+CEXE_sources += EvolveFPML.cpp
VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver