aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Docs/source/running_cpp/parameters.rst12
-rwxr-xr-xExamples/Tests/particles_in_PML/analysis.py37
-rw-r--r--Examples/Tests/particles_in_PML/inputs56
-rw-r--r--Examples/Tests/particles_in_PML/inputs2d56
-rw-r--r--Regression/WarpX-tests.ini29
-rw-r--r--Source/BoundaryConditions/PML.H44
-rw-r--r--Source/BoundaryConditions/PML.cpp346
-rw-r--r--Source/BoundaryConditions/PML_current.H129
-rw-r--r--Source/BoundaryConditions/PML_routines.F901
-rw-r--r--Source/BoundaryConditions/WarpXEvolvePML.cpp115
-rw-r--r--Source/Evolve/WarpXEvolveEM.cpp31
-rw-r--r--Source/FieldSolver/WarpXPushFieldsEM.cpp53
-rw-r--r--Source/FortranInterface/WarpX_f.H2
-rw-r--r--Source/Initialization/WarpXInitData.cpp2
-rw-r--r--Source/Parallelization/WarpXComm.cpp37
-rw-r--r--Source/WarpX.H9
-rw-r--r--Source/WarpX.cpp6
17 files changed, 839 insertions, 126 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst
index da148fde6..9beb234ae 100644
--- a/Docs/source/running_cpp/parameters.rst
+++ b/Docs/source/running_cpp/parameters.rst
@@ -663,6 +663,18 @@ Boundary conditions
The characteristic depth, in number of cells, over which
the absorption coefficients of the PML increases.
+* ``warpx.do_pml_in_domain`` (`int`; default: 0)
+ Whether to create the PML inside the simulation area or outside. If inside,
+ it allows the user to propagate particles in PML and to use extended PML
+
+* ``warpx.do_pml_has_particles`` (`int`; default: 0)
+ Whether to propagate particles in PML or not. Can only be done if PML are in simulation domain,
+ i.e. if `warpx.do_pml_in_domain = 1`.
+
+* ``warpx.do_pml_j_damping`` (`int`; default: 0)
+ Whether to damp current in PML. Can only be used if particles are propagated in PML,
+ i.e. if `warpx.do_pml_has_particles = 1`.
+
* ``warpx.do_pml_Lo`` (`2 ints in 2D`, `3 ints in 3D`; default: `1 1 1`)
The directions along which one wants a pml boundary condition for lower boundaries on mother grid.
diff --git a/Examples/Tests/particles_in_PML/analysis.py b/Examples/Tests/particles_in_PML/analysis.py
new file mode 100755
index 000000000..9d0c95d35
--- /dev/null
+++ b/Examples/Tests/particles_in_PML/analysis.py
@@ -0,0 +1,37 @@
+#! /usr/bin/env python
+"""
+This script tests the absorption of particles in the PML.
+
+The input file inputs_2d/inputs is used: it features a positive and a
+negative particle, going in opposite direction and eventually
+leaving the box. This script tests that the field in the box
+is close to 0 once the particles have left. With regular
+PML, this test fails, since the particles leave a spurious
+charge, with associated fields, behind them.
+"""
+import sys
+import yt
+import numpy as np
+yt.funcs.mylog.setLevel(0)
+
+# Open plotfile specified in command line
+filename = sys.argv[1]
+ds = yt.load( filename )
+
+# Check that the field is low enough
+ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Ex_array = ad0['Ex'].to_ndarray()
+Ey_array = ad0['Ey'].to_ndarray()
+Ez_array = ad0['Ez'].to_ndarray()
+max_Efield = max(Ex_array.max(), Ey_array.max(), Ez_array.max())
+print( max_Efield )
+
+# The field associated with the particle does not have
+# the same amplitude in 2d and 3d
+if ds.dimensionality == 2:
+ assert max_Efield < 0.0003
+elif ds.dimensionality == 3:
+ assert max_Efield < 10
+else:
+ raise ValueError("Unknown dimensionality")
+
diff --git a/Examples/Tests/particles_in_PML/inputs b/Examples/Tests/particles_in_PML/inputs
new file mode 100644
index 000000000..32d72d452
--- /dev/null
+++ b/Examples/Tests/particles_in_PML/inputs
@@ -0,0 +1,56 @@
+max_step = 120
+amr.plot_int = 30
+amr.n_cell = 128 64 64
+
+amr.blocking_factor = 8
+amr.max_grid_size = 1024
+amr.max_level = 0
+
+warpx.dump_plotfiles = 1
+
+# Geometry
+geometry.coord_sys = 0 # 0: Cartesian
+geometry.is_periodic = 0 0 0 # Is periodic?
+geometry.prob_lo = -32.e-6 -32.e-6 -32.e-6 # physical domain
+geometry.prob_hi = 32.e-6 32.e-6 32.e-6
+
+# PML
+warpx.do_pml = 1
+warpx.pml_ncell = 12
+warpx.pml_delta = 6
+warpx.pml_has_particles = 1
+warpx.do_pml_in_domain = 1
+warpx.do_pml_j_damping = 1
+
+
+# Algorithms
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = vectorized
+algo.particle_pusher = vay
+algo.maxwell_fdtd_solver = ckc
+warpx.cfl = 1.0
+warpx.use_filter = 1
+
+# Particle species
+particles.nspecies = 2
+particles.species_names = electron proton
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "singleparticle"
+electron.single_particle_pos = 0. 0. 0.
+electron.single_particle_vel = 2. 0. 0.
+electron.single_particle_weight = 1.
+
+proton.charge = q_e
+proton.mass = m_p # Very heavy ; should not move
+proton.injection_style = "singleparticle"
+proton.single_particle_pos = 0. 0. 0.
+proton.single_particle_vel = -2. 0. 0.
+proton.single_particle_weight = 1.
+
+# Particle shape factor in each direction
+interpolation.nox = 3
+interpolation.noy = 3
+interpolation.noz = 3
diff --git a/Examples/Tests/particles_in_PML/inputs2d b/Examples/Tests/particles_in_PML/inputs2d
new file mode 100644
index 000000000..076d7a55d
--- /dev/null
+++ b/Examples/Tests/particles_in_PML/inputs2d
@@ -0,0 +1,56 @@
+max_step = 180
+amr.plot_int = 30
+amr.n_cell = 224 224
+
+amr.blocking_factor = 8
+amr.max_grid_size = 1024
+amr.max_level = 0
+
+warpx.dump_plotfiles = 1
+
+# Geometry
+geometry.coord_sys = 0 # 0: Cartesian
+geometry.is_periodic = 0 0 # Is periodic?
+geometry.prob_lo = -32.e-6 -32.e-6 # physical domain
+geometry.prob_hi = 32.e-6 32.e-6
+
+# PML
+warpx.do_pml = 1
+warpx.pml_ncell = 12
+warpx.pml_delta = 6
+warpx.pml_has_particles = 1
+warpx.do_pml_in_domain = 1
+warpx.do_pml_j_damping = 1
+
+
+# Algorithms
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = vectorized
+algo.particle_pusher = vay
+algo.maxwell_fdtd_solver = ckc
+warpx.cfl = 1.0
+warpx.use_filter = 1
+
+# Particle species
+particles.nspecies = 2
+particles.species_names = electron proton
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "singleparticle"
+electron.single_particle_pos = 0. 0. 0.
+electron.single_particle_vel = 2. 0. 0.
+electron.single_particle_weight = 1.
+
+proton.charge = q_e
+proton.mass = m_p # Very heavy ; should not move
+proton.injection_style = "singleparticle"
+proton.single_particle_pos = 0. 0. 0.
+proton.single_particle_vel = -2. 0. 0.
+proton.single_particle_weight = 1.
+
+# Particle shape factor in each direction
+interpolation.nox = 3
+interpolation.noy = 3
+interpolation.noz = 3
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 0327f378f..473c5a074 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -534,3 +534,32 @@ compileTest = 0
doVis = 0
compareParticles = 1
particleTypes = electrons
+
+[particles_in_pml_2d]
+buildDir = .
+inputFile = Examples/Tests/particles_in_PML/inputs2d
+dim = 2
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/particles_in_PML/analysis.py
+
+[particles_in_pml]
+buildDir = .
+inputFile = Examples/Tests/particles_in_PML/inputs
+dim = 3
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 2
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/particles_in_PML/analysis.py
+
diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H
index b34cbe88b..9e04322f5 100644
--- a/Source/BoundaryConditions/PML.H
+++ b/Source/BoundaryConditions/PML.H
@@ -47,10 +47,15 @@ struct SigmaBox
using SigmaVect = std::array<Sigma,AMREX_SPACEDIM>;
- SigmaVect sigma; // sigma/epsilon
- SigmaVect sigma_star; // sigma_star/mu
+ SigmaVect sigma;
+ SigmaVect sigma_cumsum;
+ SigmaVect sigma_star;
+ SigmaVect sigma_star_cumsum;
SigmaVect sigma_fac;
+ SigmaVect sigma_cumsum_fac;
SigmaVect sigma_star_fac;
+ SigmaVect sigma_star_cumsum_fac;
+
};
namespace amrex {
@@ -102,14 +107,18 @@ public:
amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
#endif
int do_dive_cleaning, int do_moving_window,
- const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi);
+ int pml_has_particles, int do_pml_in_domain,
+ const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(),
+ const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector());
void ComputePMLFactors (amrex::Real dt);
std::array<amrex::MultiFab*,3> GetE_fp ();
std::array<amrex::MultiFab*,3> GetB_fp ();
+ std::array<amrex::MultiFab*,3> Getj_fp ();
std::array<amrex::MultiFab*,3> GetE_cp ();
std::array<amrex::MultiFab*,3> GetB_cp ();
+ std::array<amrex::MultiFab*,3> Getj_cp ();
amrex::MultiFab* GetF_fp ();
amrex::MultiFab* GetF_cp ();
@@ -125,16 +134,21 @@ public:
#endif
void ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp,
- const std::array<amrex::MultiFab*,3>& B_cp);
+ const std::array<amrex::MultiFab*,3>& B_cp, int do_pml_in_domain);
void ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
- const std::array<amrex::MultiFab*,3>& E_cp);
+ const std::array<amrex::MultiFab*,3>& E_cp, int do_pml_in_domain);
+ void CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp,
+ const std::array<amrex::MultiFab*,3>& j_cp);
+
void ExchangeB (PatchType patch_type,
- const std::array<amrex::MultiFab*,3>& Bp);
+ const std::array<amrex::MultiFab*,3>& Bp, int do_pml_in_domain);
void ExchangeE (PatchType patch_type,
- const std::array<amrex::MultiFab*,3>& Ep);
+ const std::array<amrex::MultiFab*,3>& Ep, int do_pml_in_domain);
+ void CopyJtoPMLs (PatchType patch_type,
+ const std::array<amrex::MultiFab*,3>& jp);
- void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp);
- void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp);
+ void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain);
+ void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain);
void FillBoundary ();
void FillBoundaryE ();
@@ -157,9 +171,11 @@ private:
std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_fp;
std::array<std::unique_ptr<amrex::MultiFab>,3> pml_B_fp;
+ std::array<std::unique_ptr<amrex::MultiFab>,3> pml_j_fp;
std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_cp;
std::array<std::unique_ptr<amrex::MultiFab>,3> pml_B_cp;
+ std::array<std::unique_ptr<amrex::MultiFab>,3> pml_j_cp;
std::unique_ptr<amrex::MultiFab> pml_F_fp;
std::unique_ptr<amrex::MultiFab> pml_F_cp;
@@ -173,11 +189,13 @@ private:
#endif
static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom,
- const amrex::BoxArray& grid_ba, int ncell,
- const amrex::IntVect do_pml_Lo,
- const amrex::IntVect do_pml_Hi);
+ const amrex::BoxArray& grid_ba,
+ int ncell, int do_pml_in_domain,
+ const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(),
+ const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector());
- static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom);
+ static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain);
+ static void CopyToPML (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom);
};
#ifdef WARPX_USE_PSATD
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 21d348482..8f8a2608e 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -1,4 +1,3 @@
-
#include <PML.H>
#include <WarpX.H>
#include <WarpXConst.H>
@@ -16,54 +15,70 @@ using namespace amrex;
namespace
{
- static void FillLo (int idim, Sigma& sigma, Sigma& sigma_star,
+ static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+ Sigma& sigma_star, Sigma& sigma_star_cumsum,
const Box& overlap, const Box& grid, Real fac)
{
int glo = grid.smallEnd(idim);
int olo = overlap.smallEnd(idim);
int ohi = overlap.bigEnd(idim);
int slo = sigma.m_lo;
+ int shi = sigma.m_hi;
int sslo = sigma_star.m_lo;
+
for (int i = olo; i <= ohi+1; ++i)
{
Real offset = static_cast<Real>(glo-i);
sigma[i-slo] = fac*(offset*offset);
+ // sigma_cumsum is the analytical integral of sigma function at same points than sigma
+ sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
}
+
for (int i = olo; i <= ohi; ++i)
{
Real offset = static_cast<Real>(glo-i) - 0.5;
sigma_star[i-sslo] = fac*(offset*offset);
+ // sigma_star_cumsum is the analytical integral of sigma function at same points than sigma_star
+ sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
}
}
- static void FillHi (int idim, Sigma& sigma, Sigma& sigma_star,
+ static void FillHi (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+ Sigma& sigma_star, Sigma& sigma_star_cumsum,
const Box& overlap, const Box& grid, Real fac)
{
int ghi = grid.bigEnd(idim);
int olo = overlap.smallEnd(idim);
int ohi = overlap.bigEnd(idim);
int slo = sigma.m_lo;
+ int shi = sigma.m_hi;
int sslo = sigma_star.m_lo;
for (int i = olo; i <= ohi+1; ++i)
{
Real offset = static_cast<Real>(i-ghi-1);
sigma[i-slo] = fac*(offset*offset);
+ sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
}
for (int i = olo; i <= ohi; ++i)
{
Real offset = static_cast<Real>(i-ghi) - 0.5;
sigma_star[i-sslo] = fac*(offset*offset);
+ sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c;
}
}
- static void FillZero (int idim, Sigma& sigma, Sigma& sigma_star, const Box& overlap)
+ static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cumsum,
+ Sigma& sigma_star, Sigma& sigma_star_cumsum,
+ const Box& overlap)
{
int olo = overlap.smallEnd(idim);
int ohi = overlap.bigEnd(idim);
int slo = sigma.m_lo;
int sslo = sigma_star.m_lo;
std::fill(sigma.begin()+(olo-slo), sigma.begin()+(ohi+2-slo), 0.0);
+ std::fill(sigma_cumsum.begin()+(olo-slo), sigma_cumsum.begin()+(ohi+2-slo), 0.0);
std::fill(sigma_star.begin()+(olo-sslo), sigma_star.begin()+(ohi+1-sslo), 0.0);
+ std::fill(sigma_star_cumsum.begin()+(olo-sslo), sigma_star_cumsum.begin()+(ohi+1-sslo), 0.0);
}
}
@@ -77,19 +92,31 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
{
- sigma [idim].resize(sz[idim]+1);
- sigma_star [idim].resize(sz[idim] );
- sigma_fac [idim].resize(sz[idim]+1);
- sigma_star_fac[idim].resize(sz[idim] );
-
- sigma [idim].m_lo = lo[idim];
- sigma [idim].m_hi = hi[idim]+1;
- sigma_star [idim].m_lo = lo[idim];
- sigma_star [idim].m_hi = hi[idim];
- sigma_fac [idim].m_lo = lo[idim];
- sigma_fac [idim].m_hi = hi[idim]+1;
- sigma_star_fac[idim].m_lo = lo[idim];
- sigma_star_fac[idim].m_hi = hi[idim];
+ sigma [idim].resize(sz[idim]+1);
+ sigma_cumsum [idim].resize(sz[idim]+1);
+ sigma_star [idim].resize(sz[idim]);
+ sigma_star_cumsum [idim].resize(sz[idim]);
+ sigma_fac [idim].resize(sz[idim]+1);
+ sigma_cumsum_fac [idim].resize(sz[idim]+1);
+ sigma_star_fac [idim].resize(sz[idim]);
+ sigma_star_cumsum_fac[idim].resize(sz[idim]);
+
+ sigma [idim].m_lo = lo[idim];
+ sigma [idim].m_hi = hi[idim]+1;
+ sigma_cumsum [idim].m_lo = lo[idim];
+ sigma_cumsum [idim].m_hi = hi[idim]+1;
+ sigma_star [idim].m_lo = lo[idim];
+ sigma_star [idim].m_hi = hi[idim];
+ sigma_star_cumsum [idim].m_lo = lo[idim];
+ sigma_star_cumsum [idim].m_hi = hi[idim];
+ sigma_fac [idim].m_lo = lo[idim];
+ sigma_fac [idim].m_hi = hi[idim]+1;
+ sigma_cumsum_fac [idim].m_lo = lo[idim];
+ sigma_cumsum_fac [idim].m_hi = hi[idim]+1;
+ sigma_star_fac [idim].m_lo = lo[idim];
+ sigma_star_fac [idim].m_hi = hi[idim];
+ sigma_star_cumsum_fac[idim].m_lo = lo[idim];
+ sigma_star_cumsum_fac[idim].m_hi = hi[idim];
}
Array<Real,AMREX_SPACEDIM> fac;
@@ -157,7 +184,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
#endif
Box looverlap = lobox & box;
if (looverlap.ok()) {
- FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]);
+ FillLo(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim],
+ looverlap, grid_box, fac[idim]);
}
Box hibox = amrex::adjCellHi(grid_box, idim, ncell);
@@ -167,7 +196,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
#endif
Box hioverlap = hibox & box;
if (hioverlap.ok()) {
- FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]);
+ FillHi(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim],
+ hioverlap, grid_box, fac[idim]);
}
if (!looverlap.ok() && !hioverlap.ok()) {
@@ -181,8 +212,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
const Box& grid_box = grids[gid];
const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box;
if (overlap.ok()) {
- FillZero(idim, sigma[idim], sigma_star[idim], overlap);
- } else {
+ FillZero(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim], overlap);
+ }
+ else {
amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n");
}
}
@@ -194,13 +227,17 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
Box lobox = amrex::adjCellLo(grid_box, idim, ncell);
Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box;
if (looverlap.ok()) {
- FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]);
+ FillLo(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim],
+ looverlap, grid_box, fac[idim]);
}
Box hibox = amrex::adjCellHi(grid_box, idim, ncell);
Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box;
if (hioverlap.ok()) {
- FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]);
+ FillHi(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim],
+ hioverlap, grid_box, fac[idim]);
}
if (!looverlap.ok() && !hioverlap.ok()) {
@@ -218,7 +255,8 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box;
#endif
if (overlap.ok()) {
- FillZero(idim, sigma[idim], sigma_star[idim], overlap);
+ FillZero(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim], overlap);
} else {
amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n");
}
@@ -231,13 +269,17 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell);
Box looverlap = lobox & box;
if (looverlap.ok()) {
- FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]);
+ FillLo(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim],
+ looverlap, grid_box, fac[idim]);
}
const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell);
Box hioverlap = hibox & box;
if (hioverlap.ok()) {
- FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]);
+ FillHi(idim, sigma[idim], sigma_cumsum[idim],
+ sigma_star[idim], sigma_star_cumsum[idim],
+ hioverlap, grid_box, fac[idim]);
}
if (!looverlap.ok() && !hioverlap.ok()) {
@@ -251,6 +293,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n
}
}
+
void
SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt)
{
@@ -259,6 +302,7 @@ SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt)
for (int i = 0, N = sigma_star[idim].size(); i < N; ++i)
{
sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt);
+ sigma_star_cumsum_fac[idim][i] = std::exp(-sigma_star_cumsum[idim][i]*dx[idim]);
}
}
}
@@ -271,6 +315,7 @@ SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt)
for (int i = 0, N = sigma[idim].size(); i < N; ++i)
{
sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt);
+ sigma_cumsum_fac[idim][i] = std::exp(-sigma_cumsum[idim][i]*dx[idim]);
}
}
}
@@ -320,11 +365,35 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
#endif
int do_dive_cleaning, int do_moving_window,
+ int pml_has_particles, int do_pml_in_domain,
const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
: m_geom(geom),
m_cgeom(cgeom)
{
- const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi);
+
+ // When `do_pml_in_domain` is true, the PML overlap with the last `ncell` of the physical domain
+ // (instead of extending `ncell` outside of the physical domain)
+ // In order to implement this, a reduced domain is created here (decreased by ncells in all direction)
+ // and passed to `MakeBoxArray`, which surrounds it by PML boxes
+ // (thus creating the PML boxes at the right position, where they overlap with the original domain)
+ Box domain0 = geom->Domain();
+ for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+ if ( ! geom->isPeriodic(idim)) {
+ if (do_pml_Lo[idim]){
+ domain0.growLo(idim, -ncell);
+ }
+ if (do_pml_Hi[idim]){
+ domain0.growHi(idim, -ncell);
+ }
+
+ }
+ }
+ const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0));
+
+ const BoxArray& ba = (do_pml_in_domain)?
+ MakeBoxArray(*geom, grid_ba_reduced, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi) :
+ MakeBoxArray(*geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
+
if (ba.size() == 0) {
m_ok = false;
return;
@@ -366,6 +435,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
pml_B_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::By_nodal_flag), dm, 2, ngb));
pml_B_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Bz_nodal_flag), dm, 2, ngb));
+
pml_E_fp[0]->setVal(0.0);
pml_E_fp[1]->setVal(0.0);
pml_E_fp[2]->setVal(0.0);
@@ -373,15 +443,30 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
pml_B_fp[1]->setVal(0.0);
pml_B_fp[2]->setVal(0.0);
+ pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb));
+ pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb));
+ pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb));
+ pml_j_fp[0]->setVal(0.0);
+ pml_j_fp[1]->setVal(0.0);
+ pml_j_fp[2]->setVal(0.0);
+
if (do_dive_cleaning)
{
pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf));
pml_F_fp->setVal(0.0);
}
- sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta));
+ if (do_pml_in_domain){
+ sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta));
+ }
+ else {
+ sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta));
+ }
+
#ifdef WARPX_USE_PSATD
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml_in_domain==false,
+ "PSATD solver cannot be used with `do_pml_in_domain`.");
const bool in_pml = true; // Tells spectral solver to use split-PML equations
const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))};
// Get the cell-centered box, with guard cells
@@ -400,7 +485,11 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
BoxArray grid_cba = grid_ba;
grid_cba.coarsen(ref_ratio);
- const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi);
+ const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain0));
+
+ const BoxArray& cba = (do_pml_in_domain) ?
+ MakeBoxArray(*cgeom, grid_cba_reduced, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi) :
+ MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
DistributionMapping cdm{cba};
@@ -422,9 +511,20 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
{
pml_F_cp.reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()), cdm, 3, ngf));
pml_F_cp->setVal(0.0);
- }
- sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta));
+ }
+ pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb));
+ pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb));
+ pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb));
+ pml_j_cp[0]->setVal(0.0);
+ pml_j_cp[1]->setVal(0.0);
+ pml_j_cp[2]->setVal(0.0);
+
+ if (do_pml_in_domain){
+ sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta));
+ } else {
+ sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta));
+ }
#ifdef WARPX_USE_PSATD
const bool in_pml = true; // Tells spectral solver to use split-PML equations
@@ -439,7 +539,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
}
BoxArray
-PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell,
+PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
+ int ncell, int do_pml_in_domain,
const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
{
Box domain = geom.Domain();
@@ -453,14 +554,18 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
}
}
}
-
BoxList bl;
for (int i = 0, N = grid_ba.size(); i < N; ++i)
{
const Box& grid_bx = grid_ba[i];
const IntVect& grid_bx_sz = grid_bx.size();
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell,
- "Consider using larger amr.blocking_factor");
+
+ if (do_pml_in_domain == 0) {
+ // Make sure that, in the case of several distinct refinement patches,
+ // the PML cells surrounding these patches cannot overlap
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell,
+ "Consider using larger amr.blocking_factor");
+ }
Box bx = grid_bx;
bx.grow(ncell);
@@ -530,6 +635,12 @@ PML::GetB_fp ()
}
std::array<MultiFab*,3>
+PML::Getj_fp ()
+{
+ return {pml_j_fp[0].get(), pml_j_fp[1].get(), pml_j_fp[2].get()};
+}
+
+std::array<MultiFab*,3>
PML::GetE_cp ()
{
return {pml_E_cp[0].get(), pml_E_cp[1].get(), pml_E_cp[2].get()};
@@ -541,6 +652,12 @@ PML::GetB_cp ()
return {pml_B_cp[0].get(), pml_B_cp[1].get(), pml_B_cp[2].get()};
}
+std::array<MultiFab*,3>
+PML::Getj_cp ()
+{
+ return {pml_j_cp[0].get(), pml_j_cp[1].get(), pml_j_cp[2].get()};
+}
+
MultiFab*
PML::GetF_fp ()
{
@@ -555,116 +672,181 @@ PML::GetF_cp ()
void
PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp,
- const std::array<amrex::MultiFab*,3>& B_cp)
+ const std::array<amrex::MultiFab*,3>& B_cp,
+ int do_pml_in_domain)
{
- ExchangeB(PatchType::fine, B_fp);
- ExchangeB(PatchType::coarse, B_cp);
+ ExchangeB(PatchType::fine, B_fp, do_pml_in_domain);
+ ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain);
}
void
PML::ExchangeB (PatchType patch_type,
- const std::array<amrex::MultiFab*,3>& Bp)
+ const std::array<amrex::MultiFab*,3>& Bp,
+ int do_pml_in_domain)
{
if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0])
{
- Exchange(*pml_B_fp[0], *Bp[0], *m_geom);
- Exchange(*pml_B_fp[1], *Bp[1], *m_geom);
- Exchange(*pml_B_fp[2], *Bp[2], *m_geom);
+ Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain);
+ Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain);
+ Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain);
}
else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0])
{
- Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom);
- Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom);
- Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom);
+ Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain);
+ Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain);
+ Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain);
}
}
void
PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
- const std::array<amrex::MultiFab*,3>& E_cp)
+ const std::array<amrex::MultiFab*,3>& E_cp,
+ int do_pml_in_domain)
{
- ExchangeE(PatchType::fine, E_fp);
- ExchangeE(PatchType::coarse, E_cp);
+ ExchangeE(PatchType::fine, E_fp, do_pml_in_domain);
+ ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain);
}
void
PML::ExchangeE (PatchType patch_type,
- const std::array<amrex::MultiFab*,3>& Ep)
+ const std::array<amrex::MultiFab*,3>& Ep,
+ int do_pml_in_domain)
{
if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0])
{
- Exchange(*pml_E_fp[0], *Ep[0], *m_geom);
- Exchange(*pml_E_fp[1], *Ep[1], *m_geom);
- Exchange(*pml_E_fp[2], *Ep[2], *m_geom);
+ Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain);
+ Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain);
+ Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain);
}
else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0])
{
- Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom);
- Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom);
- Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom);
+ Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain);
+ Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain);
+ Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain);
}
}
void
-PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp)
+PML::CopyJtoPMLs (PatchType patch_type,
+ const std::array<amrex::MultiFab*,3>& jp)
+{
+ if (patch_type == PatchType::fine && pml_j_fp[0] && jp[0])
+ {
+ CopyToPML(*pml_j_fp[0], *jp[0], *m_geom);
+ CopyToPML(*pml_j_fp[1], *jp[1], *m_geom);
+ CopyToPML(*pml_j_fp[2], *jp[2], *m_geom);
+ }
+ else if (patch_type == PatchType::coarse && pml_j_cp[0] && jp[0])
+ {
+ CopyToPML(*pml_j_cp[0], *jp[0], *m_cgeom);
+ CopyToPML(*pml_j_cp[1], *jp[1], *m_cgeom);
+ CopyToPML(*pml_j_cp[2], *jp[2], *m_cgeom);
+ }
+}
+
+void
+PML::CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp,
+ const std::array<amrex::MultiFab*,3>& j_cp)
+{
+ CopyJtoPMLs(PatchType::fine, j_fp);
+ CopyJtoPMLs(PatchType::coarse, j_cp);
+}
+
+
+void
+PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain)
{
- ExchangeF(PatchType::fine, F_fp);
- ExchangeF(PatchType::coarse, F_cp);
+ ExchangeF(PatchType::fine, F_fp, do_pml_in_domain);
+ ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain);
}
void
-PML::ExchangeF (PatchType patch_type, MultiFab* Fp)
+PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain)
{
if (patch_type == PatchType::fine && pml_F_fp && Fp) {
- Exchange(*pml_F_fp, *Fp, *m_geom);
+ Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain);
} else if (patch_type == PatchType::coarse && pml_F_cp && Fp) {
- Exchange(*pml_F_cp, *Fp, *m_cgeom);
+ Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain);
}
}
+
void
-PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom)
+PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom,
+ int do_pml_in_domain)
{
+ BL_PROFILE("PML::Exchange");
+
const IntVect& ngr = reg.nGrowVect();
const IntVect& ngp = pml.nGrowVect();
const int ncp = pml.nComp();
const auto& period = geom.periodicity();
+ // Create temporary MultiFab to copy to and from the PML
MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr);
- if (ngp.max() > 0) // Copy from pml to the ghost cells of regular data
- {
- MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0);
- MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0);
- if (ncp == 3) {
- MultiFab::Add(totpmlmf,pml,2,0,1,0);
- }
-
- MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr);
- tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period);
+ // Create the sum of the split fields, in the PML
+ MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); // Allocate
+ MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); // Sum
+ if (ncp == 3) {
+ MultiFab::Add(totpmlmf,pml,2,0,1,0); // Sum the third split component
+ }
+ // Copy from the sum of PML split field to valid cells of regular grid
+ if (do_pml_in_domain){
+ // Valid cells of the PML and of the regular grid overlap
+ // Copy from valid cells of the PML to valid cells of the regular grid
+ reg.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), IntVect(0), period);
+ } else {
+ // Valid cells of the PML only overlap with guard cells of regular grid
+ // (and outermost valid cell of the regular grid, for nodal direction)
+ // Copy from valid cells of PML to ghost cells of regular grid
+ // but avoid updating the outermost valid cell
+ if (ngr.max() > 0) {
+ MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr);
+ tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period);
#ifdef _OPENMP
#pragma omp parallel
#endif
- for (MFIter mfi(reg); mfi.isValid(); ++mfi)
- {
- const FArrayBox& src = tmpregmf[mfi];
- FArrayBox& dst = reg[mfi];
- const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox());
- for (const Box& bx : bl)
+ for (MFIter mfi(reg); mfi.isValid(); ++mfi)
{
- dst.copy(src, bx, 0, bx, 0, 1);
+ const FArrayBox& src = tmpregmf[mfi];
+ FArrayBox& dst = reg[mfi];
+ const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox());
+ // boxDiff avoids the outermost valid cell
+ for (const Box& bx : bl) {
+ dst.copy(src, bx, 0, bx, 0, 1);
+ }
}
}
}
- // Copy from regular data to PML's first component
+ // Copy from valid cells of the regular grid to guard cells of the PML
+ // (and outermost valid cell in the nodal direction)
+ // More specifically, copy from regular data to PML's first component
// Zero out the second (and third) component
- MultiFab::Copy(tmpregmf,reg,0,0,1,0);
- tmpregmf.setVal(0.0, 1, ncp-1, 0);
+ MultiFab::Copy(tmpregmf,reg,0,0,1,0); // Fill first component of tmpregmf
+ tmpregmf.setVal(0.0, 1, ncp-1, 0); // Zero out the second (and third) component
+ if (do_pml_in_domain){
+ // Where valid cells of tmpregmf overlap with PML valid cells,
+ // copy the PML (this is order to avoid overwriting PML valid cells,
+ // in the next `ParallelCopy`)
+ tmpregmf.ParallelCopy(pml,0, 0, ncp, IntVect(0), IntVect(0), period);
+ }
pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period);
}
+
+void
+PML::CopyToPML (MultiFab& pml, MultiFab& reg, const Geometry& geom)
+{
+ const IntVect& ngr = reg.nGrowVect();
+ const IntVect& ngp = pml.nGrowVect();
+ const auto& period = geom.periodicity();
+
+ pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period);
+}
+
void
PML::FillBoundary ()
{
diff --git a/Source/BoundaryConditions/PML_current.H b/Source/BoundaryConditions/PML_current.H
new file mode 100644
index 000000000..910186a96
--- /dev/null
+++ b/Source/BoundaryConditions/PML_current.H
@@ -0,0 +1,129 @@
+#ifndef PML_CURRENT_H_
+#define PML_CURRENT_H_
+
+#include <AMReX_FArrayBox.H>
+
+using namespace amrex;
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void push_ex_pml_current (int j, int k, int l, Array4<Real> const& Ex,
+ Array4<Real const> const& jx,
+ Real const * const sigjy,
+ Real const * const sigjz,
+ int ylo, int zlo,
+ Real mu_c2_dt)
+{
+#if (AMREX_SPACEDIM == 3)
+ Real alpha_xy, alpha_xz;
+ if (sigjy[k-ylo]+sigjz[l-zlo] == 0){
+ alpha_xy = 0.5;
+ alpha_xz = 0.5;
+ }
+ else {
+ alpha_xy = sigjy[k-ylo]/(sigjy[k-ylo]+sigjz[l-zlo]);
+ alpha_xz = sigjz[l-zlo]/(sigjy[k-ylo]+sigjz[l-zlo]);
+ }
+ Ex(j,k,l,0) = Ex(j,k,l,0) - mu_c2_dt * alpha_xy * jx(j,k,l);
+ Ex(j,k,l,1) = Ex(j,k,l,1) - mu_c2_dt * alpha_xz * jx(j,k,l);
+#else
+ Ex(j,k,l,1) = Ex(j,k,l,1) - mu_c2_dt * jx(j,k,l);
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void push_ey_pml_current (int j, int k, int l, Array4<Real> const& Ey,
+ Array4<Real const> const& jy,
+ Real const * const sigjx,
+ Real const * const sigjz,
+ int xlo, int zlo,
+ Real mu_c2_dt)
+{
+#if (AMREX_SPACEDIM == 3)
+ Real alpha_yx, alpha_yz;
+ if (sigjx[j-xlo]+sigjz[l-zlo] == 0){
+ alpha_yx = 0.5;
+ alpha_yz = 0.5;
+ }
+ else {
+ alpha_yx = sigjx[j-xlo]/(sigjx[j-xlo]+sigjz[l-zlo]);
+ alpha_yz = sigjz[l-zlo]/(sigjx[j-xlo]+sigjz[l-zlo]);
+ }
+ Ey(j,k,l,0) = Ey(j,k,l,0) - mu_c2_dt * alpha_yx * jy(j,k,l);
+ Ey(j,k,l,1) = Ey(j,k,l,1) - mu_c2_dt * alpha_yz * jy(j,k,l);
+#else
+ Ey(j,k,l,0) = Ey(j,k,l,0) - 0.5 * mu_c2_dt * jy(j,k,l);
+ Ey(j,k,l,1) = Ey(j,k,l,1) - 0.5 * mu_c2_dt * jy(j,k,l);
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void push_ez_pml_current (int j, int k, int l, Array4<Real> const& Ez,
+ Array4<Real const> const& jz,
+ Real const * const sigjx,
+ Real const * const sigjy,
+ int xlo, int ylo,
+ Real mu_c2_dt)
+{
+#if (AMREX_SPACEDIM == 3)
+ Real alpha_zx, alpha_zy;
+ if (sigjx[j-xlo]+sigjy[k-ylo]==0){
+ alpha_zx = 0.5;
+ alpha_zy = 0.5;
+ }
+ else {
+ alpha_zx = sigjx[j-xlo]/(sigjx[j-xlo]+sigjy[k-ylo]);
+ alpha_zy = sigjy[k-ylo]/(sigjx[j-xlo]+sigjy[k-ylo]);
+ }
+ Ez(j,k,l,0) = Ez(j,k,l,0) - mu_c2_dt * alpha_zx * jz(j,k,l);
+ Ez(j,k,l,1) = Ez(j,k,l,1) - mu_c2_dt * alpha_zy * jz(j,k,l);
+#else
+ Ez(j,k,l,0) = Ez(j,k,l,0) - mu_c2_dt * jz(j,k,l);
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void damp_jx_pml (int j, int k, int l,
+ Array4<Real> const& jx,
+ Real const* const sigsjx,
+ Real const* const sigjy,
+ Real const* const sigjz,
+ int xlo, int ylo, int zlo)
+{
+#if (AMREX_SPACEDIM == 3)
+ jx(j,k,l) = jx(j,k,l) * sigsjx[j-xlo] * sigjy[k-ylo] * sigjz[l-zlo];
+#else
+ jx(j,k,l) = jx(j,k,l) * sigsjx[j-xlo] * sigjz[k-zlo];
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void damp_jy_pml (int j, int k, int l,
+ Array4<Real> const& jy,
+ Real const * const sigjx,
+ Real const * const sigsjy,
+ Real const * const sigjz,
+ int xlo, int ylo, int zlo)
+{
+#if (AMREX_SPACEDIM == 3)
+ jy(j,k,l) = jy(j,k,l) * sigjx[j-xlo] * sigsjy[k-ylo] * sigjz[l-zlo];
+#else
+ jy(j,k,l) = jy(j,k,l) * sigjx[j-xlo] * sigjz[k-zlo];
+#endif
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void damp_jz_pml (int j, int k, int l,
+ Array4<Real> const& jz,
+ Real const * const sigjx,
+ Real const * const sigjy,
+ Real const * const sigsjz,
+ int xlo, int ylo, int zlo)
+{
+#if (AMREX_SPACEDIM == 3)
+ jz(j,k,l) = jz(j,k,l) * sigjx[j-xlo] * sigjy[k-ylo] * sigsjz[l-zlo];
+#else
+ jz(j,k,l) = jz(j,k,l) * sigjx[j-xlo] * sigsjz[k-zlo];
+#endif
+}
+
+#endif
diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90
index 380e52934..e2540350e 100644
--- a/Source/BoundaryConditions/PML_routines.F90
+++ b/Source/BoundaryConditions/PML_routines.F90
@@ -1006,4 +1006,5 @@ contains
end do
end subroutine warpx_damp_pml_f_3d
+
end module warpx_pml_module
diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp
index b0688b2c1..f5c231ddf 100644
--- a/Source/BoundaryConditions/WarpXEvolvePML.cpp
+++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp
@@ -12,6 +12,8 @@
#include <AMReX_AmrMeshInSituBridge.H>
#endif
+#include <PML_current.H>
+
using namespace amrex;
void
@@ -55,7 +57,6 @@ WarpX::DampPML (int lev, PatchType patch_type)
const Box& tbx = mfi.tilebox(Bx_nodal_flag);
const Box& tby = mfi.tilebox(By_nodal_flag);
const Box& tbz = mfi.tilebox(Bz_nodal_flag);
-
WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
tey.loVect(), tey.hiVect(),
tez.loVect(), tez.hiVect(),
@@ -79,3 +80,115 @@ WarpX::DampPML (int lev, PatchType patch_type)
}
}
}
+
+void
+WarpX::DampJPML ()
+{
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ DampJPML(lev);
+ }
+}
+
+void
+WarpX::DampJPML (int lev)
+{
+ DampJPML(lev, PatchType::fine);
+ if (lev > 0) DampJPML(lev, PatchType::coarse);
+}
+
+void
+WarpX::DampJPML (int lev, PatchType patch_type)
+{
+ if (!do_pml) return;
+ if (!do_pml_j_damping) return;
+
+ BL_PROFILE("WarpX::DampJPML()");
+
+ if (pml[lev]->ok())
+ {
+
+ const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp();
+ const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp()
+ : pml[lev]->GetMultiSigmaBox_cp();
+
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+ for ( MFIter mfi(*pml_j[0], TilingIfNotGPU()); mfi.isValid(); ++mfi )
+ {
+ auto const& pml_jxfab = pml_j[0]->array(mfi);
+ auto const& pml_jyfab = pml_j[1]->array(mfi);
+ auto const& pml_jzfab = pml_j[2]->array(mfi);
+ const Real* sigma_cumsum_fac_j_x = sigba[mfi].sigma_cumsum_fac[0].data();
+ const Real* sigma_star_cumsum_fac_j_x = sigba[mfi].sigma_star_cumsum_fac[0].data();
+#if (AMREX_SPACEDIM == 3)
+ const Real* sigma_cumsum_fac_j_y = sigba[mfi].sigma_cumsum_fac[1].data();
+ const Real* sigma_star_cumsum_fac_j_y = sigba[mfi].sigma_star_cumsum_fac[1].data();
+ const Real* sigma_cumsum_fac_j_z = sigba[mfi].sigma_cumsum_fac[2].data();
+ const Real* sigma_star_cumsum_fac_j_z = sigba[mfi].sigma_star_cumsum_fac[2].data();
+#else
+ const Real* sigma_cumsum_fac_j_y = nullptr;
+ const Real* sigma_star_cumsum_fac_j_y = nullptr;
+ const Real* sigma_cumsum_fac_j_z = sigba[mfi].sigma_cumsum_fac[1].data();
+ const Real* sigma_star_cumsum_fac_j_z = sigba[mfi].sigma_star_cumsum_fac[1].data();
+#endif
+ const Box& tjx = mfi.tilebox(jx_nodal_flag);
+ const Box& tjy = mfi.tilebox(jy_nodal_flag);
+ const Box& tjz = mfi.tilebox(jz_nodal_flag);
+
+ auto const& AMREX_RESTRICT x_lo = sigba[mfi].sigma_cumsum_fac[0].lo();
+#if (AMREX_SPACEDIM == 3)
+ auto const& AMREX_RESTRICT y_lo = sigba[mfi].sigma_cumsum_fac[1].lo();
+ auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma_cumsum_fac[2].lo();
+#else
+ int y_lo = 0;
+ auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma_cumsum_fac[1].lo();
+#endif
+
+ auto const& AMREX_RESTRICT xs_lo = sigba[mfi].sigma_star_cumsum_fac[0].lo();
+#if (AMREX_SPACEDIM == 3)
+ auto const& AMREX_RESTRICT ys_lo = sigba[mfi].sigma_star_cumsum_fac[1].lo();
+ auto const& AMREX_RESTRICT zs_lo = sigba[mfi].sigma_star_cumsum_fac[2].lo();
+#else
+ int ys_lo = 0;
+ auto const& AMREX_RESTRICT zs_lo = sigba[mfi].sigma_star_cumsum_fac[1].lo();
+#endif
+
+ amrex::ParallelFor( tjx, tjy, tjz,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ damp_jx_pml(i, j, k, pml_jxfab, sigma_star_cumsum_fac_j_x,
+ sigma_cumsum_fac_j_y, sigma_cumsum_fac_j_z,
+ xs_lo,y_lo, z_lo);
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ damp_jy_pml(i, j, k, pml_jyfab, sigma_cumsum_fac_j_x,
+ sigma_star_cumsum_fac_j_y, sigma_cumsum_fac_j_z,
+ x_lo,ys_lo, z_lo);
+ },
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ damp_jz_pml(i, j, k, pml_jzfab, sigma_cumsum_fac_j_x,
+ sigma_cumsum_fac_j_y, sigma_star_cumsum_fac_j_z,
+ x_lo,y_lo, zs_lo);
+ }
+ );
+ }
+
+ }
+}
+
+/* \brief Copy the current J from the regular grid to the PML */
+void
+WarpX::CopyJPML ()
+{
+ for (int lev = 0; lev <= finest_level; ++lev)
+ {
+ if (pml[lev]->ok()){
+ pml[lev]->CopyJtoPMLs({ current_fp[lev][0].get(),
+ current_fp[lev][1].get(),
+ current_fp[lev][2].get() },
+ { current_cp[lev][0].get(),
+ current_cp[lev][1].get(),
+ current_cp[lev][2].get() });
+ }
+ }
+}
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp
index 75643d748..594ea55ee 100644
--- a/Source/Evolve/WarpXEvolveEM.cpp
+++ b/Source/Evolve/WarpXEvolveEM.cpp
@@ -13,6 +13,7 @@
#include <AMReX_AmrMeshInSituBridge.H>
#endif
+
using namespace amrex;
void
@@ -135,7 +136,7 @@ WarpX::EvolveEM (int numsteps)
bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0);
// slice generation //
- bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0);
+ bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0);
bool do_insitu = ((step+1) >= insitu_start) &&
(insitu_int > 0) && ((step+1) % insitu_int == 0);
@@ -153,7 +154,7 @@ WarpX::EvolveEM (int numsteps)
// We might need to move j because we are going to make a plotfile.
int num_moved = MoveWindow(move_j);
-
+
if (max_level == 0) {
int num_redistribute_ghost = num_moved + 1;
mypc->RedistributeLocal(num_redistribute_ghost);
@@ -228,7 +229,7 @@ WarpX::EvolveEM (int numsteps)
// End loop on time steps
}
- bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step
+ bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step
&& (max_time_reached || istep[0] >= max_step);
bool do_insitu = (insitu_start >= istep[0]) && (insitu_int > 0) &&
@@ -255,7 +256,7 @@ WarpX::EvolveEM (int numsteps)
UpdateInSitu();
}
- if (check_int > 0 && istep[0] > last_check_file_step &&
+ if (check_int > 0 && istep[0] > last_check_file_step &&
(max_time_reached || istep[0] >= max_step)) {
WriteCheckPointFile();
}
@@ -298,6 +299,10 @@ WarpX::OneStep_nosub (Real cur_time)
SyncRho();
+ // For extended PML: copy J from regular grid to PML, and damp J in PML
+ if (do_pml && pml_has_particles) CopyJPML();
+ if (do_pml && do_pml_j_damping) DampJPML();
+
// Push E and B from {n} to {n+1}
// (And update guard cells immediately afterwards)
#ifdef WARPX_USE_PSATD
@@ -310,6 +315,7 @@ WarpX::OneStep_nosub (Real cur_time)
FillBoundaryF();
EvolveB(0.5*dt[0]); // We now have B^{n+1/2}
FillBoundaryB();
+
EvolveE(dt[0]); // We now have E^{n+1}
FillBoundaryE();
EvolveF(0.5*dt[0], DtType::SecondHalf);
@@ -319,6 +325,7 @@ WarpX::OneStep_nosub (Real cur_time)
FillBoundaryE();
}
FillBoundaryB();
+
#endif
}
@@ -545,13 +552,13 @@ WarpX::ComputeDt ()
/* \brief computes max_step for wakefield simulation in boosted frame.
* \param geom: Geometry object that contains simulation domain.
- *
- * max_step is set so that the simulation stop when the lower corner of the
+ *
+ * max_step is set so that the simulation stop when the lower corner of the
* simulation box passes input parameter zmax_plasma_to_compute_max_step.
*/
void
WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){
- // Sanity checks: can use zmax_plasma_to_compute_max_step only if
+ // Sanity checks: can use zmax_plasma_to_compute_max_step only if
// the moving window and the boost are all in z direction.
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
WarpX::moving_window_dir == AMREX_SPACEDIM-1,
@@ -565,7 +572,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){
"Can use zmax_plasma_to_compute_max_step only if " +
"warpx.boost_direction = z. TODO: all directions.");
- // Lower end of the simulation domain. All quantities are given in boosted
+ // Lower end of the simulation domain. All quantities are given in boosted
// frame except zmax_plasma_to_compute_max_step.
const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1);
// End of the plasma: Transform input argument
@@ -573,7 +580,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){
const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost;
// Plasma velocity
const Real v_plasma_boost = -beta_boost * PhysConst::c;
- // Get time at which the lower end of the simulation domain passes the
+ // Get time at which the lower end of the simulation domain passes the
// upper end of the plasma (in the z direction).
const Real interaction_time_boost = (len_plasma_boost-zmin_domain_boost)/
(moving_window_v-v_plasma_boost);
@@ -591,7 +598,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){
/* \brief Apply perfect mirror condition inside the box (not at a boundary).
* In practice, set all fields to 0 on a section of the simulation domain
- * (as for a perfect conductor with a given thickness).
+ * (as for a perfect conductor with a given thickness).
* The mirror normal direction has to be parallel to the z axis.
*/
void
@@ -608,10 +615,10 @@ WarpX::applyMirrors(Real time){
}
// Loop over levels
for(int lev=0; lev<=finest_level; lev++){
- // Make sure that the mirror contains at least
+ // Make sure that the mirror contains at least
// mirror_z_npoints[i_mirror] cells
Real dz = WarpX::CellSize(lev)[2];
- Real z_max = std::max(z_max_tmp,
+ Real z_max = std::max(z_max_tmp,
z_min+mirror_z_npoints[i_mirror]*dz);
// Get fine patch field MultiFabs
MultiFab& Ex = *Efield_fp[lev][0].get();
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 1d4491881..11d598b98 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -11,6 +11,8 @@
#include <WarpX_py.H>
#endif
+#include <PML_current.H>
+
#ifdef BL_USE_SENSEI_INSITU
#include <AMReX_AmrMeshInSituBridge.H>
#endif
@@ -445,11 +447,14 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
if (do_pml && pml[lev]->ok())
{
- if (F) pml[lev]->ExchangeF(patch_type, F);
+ if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain);
const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+ const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp();
const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+ const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp()
+ : pml[lev]->GetMultiSigmaBox_cp();
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
@@ -459,6 +464,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
const Box& tey = mfi.tilebox(Ey_nodal_flag);
const Box& tez = mfi.tilebox(Ez_nodal_flag);
+ auto const& pml_Exfab = pml_E[0]->array(mfi);
+ auto const& pml_Eyfab = pml_E[1]->array(mfi);
+ auto const& pml_Ezfab = pml_E[2]->array(mfi);
+
WRPX_PUSH_PML_EVEC(
tex.loVect(), tex.hiVect(),
tey.loVect(), tey.hiVect(),
@@ -469,7 +478,45 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- &dtsdx_c2, &dtsdy_c2, &dtsdz_c2);
+ &dtsdx_c2, &dtsdy_c2, &dtsdz_c2);
+
+ if (pml_has_particles) {
+ // Update the E field in the PML, using the current
+ // deposited by the particles in the PML
+ auto const& pml_jxfab = pml_j[0]->array(mfi);
+ auto const& pml_jyfab = pml_j[1]->array(mfi);
+ auto const& pml_jzfab = pml_j[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();
+
+ auto const& AMREX_RESTRICT x_lo = sigba[mfi].sigma[0].lo();
+#if (AMREX_SPACEDIM == 3)
+ auto const& AMREX_RESTRICT y_lo = sigba[mfi].sigma[1].lo();
+ auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma[2].lo();
+#else
+ int y_lo = 0;
+ auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma[1].lo();
+#endif
+
+ amrex::ParallelFor( tex, tey, tez,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ push_ex_pml_current(i,j,k,
+ pml_Exfab, pml_jxfab, 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,
+ pml_Eyfab, pml_jyfab, 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,
+ pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y,
+ x_lo, y_lo, mu_c2_dt);
+ }
+ );
+ }
if (pml_F)
{
@@ -481,7 +528,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
BL_TO_FORTRAN_3D((*pml_F )[mfi]),
- &dtsdx_c2, &dtsdy_c2, &dtsdz_c2,
+ &dtsdx_c2, &dtsdy_c2, &dtsdz_c2,
&WarpX::maxwell_fdtd_solver_id);
}
}
diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H
index 9db54ba18..60a9817a2 100644
--- a/Source/FortranInterface/WarpX_f.H
+++ b/Source/FortranInterface/WarpX_f.H
@@ -1,4 +1,3 @@
-
#include <AMReX_BLFort.H>
#ifdef __cplusplus
@@ -270,6 +269,7 @@ extern "C"
#endif
const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi);
+
void WRPX_SYNC_CURRENT (const int* lo, const int* hi,
BL_FORT_FAB_ARG_ANYD(crse),
const BL_FORT_FAB_ARG_ANYD(fine),
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 590c11b84..4c0dfee17 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -142,6 +142,7 @@ WarpX::InitPML ()
dt[0], nox_fft, noy_fft, noz_fft, do_nodal,
#endif
do_dive_cleaning, do_moving_window,
+ pml_has_particles, do_pml_in_domain,
do_pml_Lo_corrected, do_pml_Hi));
for (int lev = 1; lev <= finest_level; ++lev)
{
@@ -159,6 +160,7 @@ WarpX::InitPML ()
dt[lev], nox_fft, noy_fft, noz_fft, do_nodal,
#endif
do_dive_cleaning, do_moving_window,
+ pml_has_particles, do_pml_in_domain,
do_pml_Lo_MR, amrex::IntVect::TheUnitVector()));
}
}
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 9d85783b0..feb9d3564 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -18,7 +18,8 @@ WarpX::ExchangeWithPmlB (int lev)
Bfield_fp[lev][2].get() },
{ Bfield_cp[lev][0].get(),
Bfield_cp[lev][1].get(),
- Bfield_cp[lev][2].get() });
+ Bfield_cp[lev][2].get() },
+ do_pml_in_domain);
}
}
@@ -31,7 +32,8 @@ WarpX::ExchangeWithPmlE (int lev)
Efield_fp[lev][2].get() },
{ Efield_cp[lev][0].get(),
Efield_cp[lev][1].get(),
- Efield_cp[lev][2].get() });
+ Efield_cp[lev][2].get() },
+ do_pml_in_domain);
}
}
@@ -40,7 +42,8 @@ WarpX::ExchangeWithPmlF (int lev)
{
if (do_pml && pml[lev]->ok()) {
pml[lev]->ExchangeF(F_fp[lev].get(),
- F_cp[lev].get());
+ F_cp[lev].get(),
+ do_pml_in_domain);
}
}
@@ -250,9 +253,10 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type)
if (do_pml && pml[lev]->ok())
{
pml[lev]->ExchangeE(patch_type,
- { Efield_fp[lev][0].get(),
+ { Efield_fp[lev][0].get(),
Efield_fp[lev][1].get(),
- Efield_fp[lev][2].get() });
+ Efield_fp[lev][2].get() },
+ do_pml_in_domain);
pml[lev]->FillBoundaryE(patch_type);
}
@@ -265,9 +269,10 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type)
if (do_pml && pml[lev]->ok())
{
pml[lev]->ExchangeE(patch_type,
- { Efield_cp[lev][0].get(),
+ { Efield_cp[lev][0].get(),
Efield_cp[lev][1].get(),
- Efield_cp[lev][2].get() });
+ Efield_cp[lev][2].get() },
+ do_pml_in_domain);
pml[lev]->FillBoundaryE(patch_type);
}
@@ -292,9 +297,10 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type)
if (do_pml && pml[lev]->ok())
{
pml[lev]->ExchangeB(patch_type,
- { Bfield_fp[lev][0].get(),
+ { Bfield_fp[lev][0].get(),
Bfield_fp[lev][1].get(),
- Bfield_fp[lev][2].get() });
+ Bfield_fp[lev][2].get() },
+ do_pml_in_domain);
pml[lev]->FillBoundaryB(patch_type);
}
const auto& period = Geom(lev).periodicity();
@@ -306,9 +312,10 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type)
if (do_pml && pml[lev]->ok())
{
pml[lev]->ExchangeB(patch_type,
- { Bfield_cp[lev][0].get(),
- Bfield_cp[lev][1].get(),
- Bfield_cp[lev][2].get() });
+ { Bfield_cp[lev][0].get(),
+ Bfield_cp[lev][1].get(),
+ Bfield_cp[lev][2].get() },
+ do_pml_in_domain);
pml[lev]->FillBoundaryB(patch_type);
}
const auto& cperiod = Geom(lev-1).periodicity();
@@ -331,7 +338,8 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type)
{
if (do_pml && pml[lev]->ok())
{
- pml[lev]->ExchangeF(patch_type, F_fp[lev].get());
+ pml[lev]->ExchangeF(patch_type, F_fp[lev].get(),
+ do_pml_in_domain);
pml[lev]->FillBoundaryF(patch_type);
}
@@ -342,7 +350,8 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type)
{
if (do_pml && pml[lev]->ok())
{
- pml[lev]->ExchangeF(patch_type, F_cp[lev].get());
+ pml[lev]->ExchangeF(patch_type, F_cp[lev].get(),
+ do_pml_in_domain);
pml[lev]->FillBoundaryF(patch_type);
}
diff --git a/Source/WarpX.H b/Source/WarpX.H
index f1a3b1479..36d2a8f35 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -191,6 +191,12 @@ public:
void DampPML (int lev);
void DampPML (int lev, PatchType patch_type);
+ void DampJPML ();
+ void DampJPML (int lev);
+ void DampJPML (int lev, PatchType patch_type);
+
+ void CopyJPML ();
+
void PushParticlesandDepose (int lev, amrex::Real cur_time);
void PushParticlesandDepose ( amrex::Real cur_time);
@@ -499,6 +505,9 @@ private:
int do_pml = 1;
int pml_ncell = 10;
int pml_delta = 10;
+ int pml_has_particles = 0;
+ int do_pml_j_damping = 0;
+ int do_pml_in_domain = 0;
amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector();
amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector();
amrex::Vector<std::unique_ptr<PML> > pml;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 4723626ab..517fb2332 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -394,6 +394,9 @@ WarpX::ReadParameters ()
pp.query("do_pml", do_pml);
pp.query("pml_ncell", pml_ncell);
pp.query("pml_delta", pml_delta);
+ pp.query("pml_has_particles", pml_has_particles);
+ pp.query("do_pml_j_damping", do_pml_j_damping);
+ pp.query("do_pml_in_domain", do_pml_in_domain);
Vector<int> parse_do_pml_Lo(AMREX_SPACEDIM,1);
pp.queryarr("do_pml_Lo", parse_do_pml_Lo);
@@ -410,6 +413,9 @@ WarpX::ReadParameters ()
do_pml_Hi[2] = parse_do_pml_Hi[2];
#endif
+ if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){
+ amrex::Abort("J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)");
+ }
pp.query("dump_openpmd", dump_openpmd);
pp.query("openpmd_backend", openpmd_backend);
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-13 01:22:11 +0000