aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/BoundaryConditions/PML.H20
-rw-r--r--Source/BoundaryConditions/PML.cpp51
-rw-r--r--Source/Initialization/WarpXInitData.cpp2
-rw-r--r--Source/WarpX.H2
-rw-r--r--Source/WarpX.cpp7
5 files changed, 50 insertions, 32 deletions
diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H
index 20f1c3d43..d0966b312 100644
--- a/Source/BoundaryConditions/PML.H
+++ b/Source/BoundaryConditions/PML.H
@@ -43,13 +43,15 @@ struct SigmaBox
{
SigmaBox (const amrex::Box& box, const amrex::BoxArray& grids,
const amrex::Real* dx, const amrex::IntVect& ncell, const amrex::IntVect& delta,
- const amrex::Box& regdomain);
+ const amrex::Box& regdomain, const amrex::Real v_sigma);
void define_single (const amrex::Box& regdomain, const amrex::IntVect& ncell,
- const amrex::Array<amrex::Real,AMREX_SPACEDIM>& fac);
+ const amrex::Array<amrex::Real,AMREX_SPACEDIM>& fac,
+ const amrex::Real v_sigma);
void define_multiple (const amrex::Box& box, const amrex::BoxArray& grids,
const amrex::IntVect& ncell,
- const amrex::Array<amrex::Real,AMREX_SPACEDIM>& fac);
+ const amrex::Array<amrex::Real,AMREX_SPACEDIM>& fac,
+ const amrex::Real v_sigma);
void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt);
void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt);
@@ -64,6 +66,7 @@ struct SigmaBox
SigmaVect sigma_cumsum_fac;
SigmaVect sigma_star_fac;
SigmaVect sigma_star_cumsum_fac;
+ amrex::Real v_sigma;
};
@@ -73,8 +76,8 @@ class SigmaBoxFactory
public:
SigmaBoxFactory (const amrex::BoxArray& grid_ba, const amrex::Real* dx,
const amrex::IntVect& ncell, const amrex::IntVect& delta,
- const amrex::Box& regular_domain)
- : m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta), m_regdomain(regular_domain) {}
+ const amrex::Box& regular_domain, const amrex::Real v_sigma_sb)
+ : m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta), m_regdomain(regular_domain), m_v_sigma_sb(v_sigma_sb) {}
virtual ~SigmaBoxFactory () = default;
SigmaBoxFactory (const SigmaBoxFactory&) = default;
@@ -86,7 +89,7 @@ public:
virtual SigmaBox* create (const amrex::Box& box, int /*ncomps*/,
const amrex::FabInfo& /*info*/, int /*box_index*/) const final
- { return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta, m_regdomain); }
+ { return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta, m_regdomain, m_v_sigma_sb); }
virtual void destroy (SigmaBox* fab) const final {
delete fab;
}
@@ -99,6 +102,7 @@ private:
amrex::IntVect m_ncell;
amrex::IntVect m_delta;
amrex::Box m_regdomain;
+ amrex::Real m_v_sigma_sb;
};
class MultiSigmaBox
@@ -108,7 +112,7 @@ public:
MultiSigmaBox(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
const amrex::BoxArray& grid_ba, const amrex::Real* dx,
const amrex::IntVect& ncell, const amrex::IntVect& delta,
- const amrex::Box& regular_domain);
+ const amrex::Box& regular_domain, const amrex::Real v_sigma_sb);
void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt);
void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt);
private:
@@ -129,7 +133,7 @@ public:
int do_moving_window, int pml_has_particles, int do_pml_in_domain,
const bool do_multi_J,
const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning,
- int max_guard_EB,
+ int max_guard_EB, amrex::Real v_sigma_sb,
const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(),
const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector());
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 9529f1f11..67eef327b 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -58,7 +58,8 @@ namespace
{
static void FillLo (Sigma& sigma, Sigma& sigma_cumsum,
Sigma& sigma_star, Sigma& sigma_star_cumsum,
- const int olo, const int ohi, const int glo, Real fac)
+ const int olo, const int ohi, const int glo, Real fac,
+ const amrex::Real v_sigma)
{
const int slo = sigma.m_lo;
const int sslo = sigma_star.m_lo;
@@ -74,19 +75,20 @@ namespace
Real offset = static_cast<Real>(glo-i);
p_sigma[i-slo] = fac*(offset*offset);
// sigma_cumsum is the analytical integral of sigma function at same points than sigma
- p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/PhysConst::c;
+ p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
if (i <= ohi+1) {
offset = static_cast<Real>(glo-i) - 0.5_rt;
p_sigma_star[i-sslo] = fac*(offset*offset);
// sigma_star_cumsum is the analytical integral of sigma function at same points than sigma_star
- p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/PhysConst::c;
+ p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
}
});
}
static void FillHi (Sigma& sigma, Sigma& sigma_cumsum,
Sigma& sigma_star, Sigma& sigma_star_cumsum,
- const int olo, const int ohi, const int ghi, Real fac)
+ const int olo, const int ohi, const int ghi, Real fac,
+ const amrex::Real v_sigma)
{
const int slo = sigma.m_lo;
const int sslo = sigma_star.m_lo;
@@ -101,11 +103,11 @@ namespace
i += olo;
Real offset = static_cast<Real>(i-ghi-1);
p_sigma[i-slo] = fac*(offset*offset);
- p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/PhysConst::c;
+ p_sigma_cumsum[i-slo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
if (i <= ohi+1) {
offset = static_cast<Real>(i-ghi) - 0.5_rt;
p_sigma_star[i-sslo] = fac*(offset*offset);
- p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/PhysConst::c;
+ p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset*offset)/3._rt)/v_sigma;
}
});
}
@@ -139,7 +141,7 @@ namespace
SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const IntVect& ncell,
- const IntVect& delta, const amrex::Box& regdomain)
+ const IntVect& delta, const amrex::Box& regdomain, const amrex::Real v_sigma_sb)
{
BL_ASSERT(box.cellCentered());
@@ -182,14 +184,15 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const
}
if (regdomain.ok()) { // The union of the regular grids is a single box
- define_single(regdomain, ncell, fac);
+ define_single(regdomain, ncell, fac, v_sigma_sb);
} else {
- define_multiple(box, grids, ncell, fac);
+ define_multiple(box, grids, ncell, fac, v_sigma_sb);
}
}
void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell,
- const Array<Real,AMREX_SPACEDIM>& fac)
+ const Array<Real,AMREX_SPACEDIM>& fac,
+ const amrex::Real v_sigma_sb)
{
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
const int slo = sigma[idim].lo();
@@ -203,7 +206,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell,
if (ohi >= olo) {
FillLo(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
- olo, ohi, dlo, fac[idim]);
+ olo, ohi, dlo, fac[idim], v_sigma_sb);
}
#if (AMREX_SPACEDIM != 1)
@@ -223,7 +226,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell,
if (ohi >= olo) {
FillHi(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
- olo, ohi, dhi, fac[idim]);
+ olo, ohi, dhi, fac[idim], v_sigma_sb);
}
}
@@ -231,7 +234,7 @@ void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell,
}
void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const IntVect& ncell,
- const Array<Real,AMREX_SPACEDIM>& fac)
+ const Array<Real,AMREX_SPACEDIM>& fac, const amrex::Real v_sigma_sb)
{
const std::vector<std::pair<int,Box> >& isects = grids.intersections(box, false, ncell);
@@ -301,7 +304,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int
FillLo(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
looverlap.smallEnd(idim), looverlap.bigEnd(idim),
- grid_box.smallEnd(idim), fac[idim]);
+ grid_box.smallEnd(idim), fac[idim], v_sigma_sb);
}
Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]);
@@ -314,7 +317,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int
FillHi(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
hioverlap.smallEnd(idim), hioverlap.bigEnd(idim),
- grid_box.bigEnd(idim), fac[idim]);
+ grid_box.bigEnd(idim), fac[idim], v_sigma_sb);
}
if (!looverlap.ok() && !hioverlap.ok()) {
@@ -348,7 +351,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int
FillLo(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
looverlap.smallEnd(idim), looverlap.bigEnd(idim),
- grid_box.smallEnd(idim), fac[idim]);
+ grid_box.smallEnd(idim), fac[idim], v_sigma_sb);
}
Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]);
@@ -357,7 +360,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int
FillHi(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
hioverlap.smallEnd(idim), hioverlap.bigEnd(idim),
- grid_box.bigEnd(idim), fac[idim]);
+ grid_box.bigEnd(idim), fac[idim], v_sigma_sb);
}
if (!looverlap.ok() && !hioverlap.ok()) {
@@ -395,7 +398,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int
FillLo(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
looverlap.smallEnd(idim), looverlap.bigEnd(idim),
- grid_box.smallEnd(idim), fac[idim]);
+ grid_box.smallEnd(idim), fac[idim], v_sigma_sb);
}
const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]);
@@ -404,7 +407,7 @@ void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const Int
FillHi(sigma[idim], sigma_cumsum[idim],
sigma_star[idim], sigma_star_cumsum[idim],
hioverlap.smallEnd(idim), hioverlap.bigEnd(idim),
- grid_box.bigEnd(idim), fac[idim]);
+ grid_box.bigEnd(idim), fac[idim], v_sigma_sb);
}
if (!looverlap.ok() && !hioverlap.ok()) {
@@ -492,9 +495,9 @@ SigmaBox::ComputePMLFactorsE (const Real* a_dx, Real dt)
MultiSigmaBox::MultiSigmaBox (const BoxArray& ba, const DistributionMapping& dm,
const BoxArray& grid_ba, const Real* dx,
const IntVect& ncell, const IntVect& delta,
- const amrex::Box& regular_domain)
+ const amrex::Box& regular_domain, const amrex::Real v_sigma_sb)
: FabArray<SigmaBox>(ba,dm,1,0,MFInfo(),
- SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain))
+ SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain, v_sigma_sb))
{}
void
@@ -536,7 +539,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain,
const bool do_multi_J,
const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning,
- int max_guard_EB,
+ int max_guard_EB, const amrex::Real v_sigma_sb,
const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
: m_dive_cleaning(do_pml_dive_cleaning),
m_divb_cleaning(do_pml_divb_cleaning),
@@ -716,7 +719,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
Box single_domain_box = is_single_box_domain ? domain0 : Box();
// Empty box (i.e., Box()) means it's not a single box domain.
sigba_fp = std::make_unique<MultiSigmaBox>(ba, dm, grid_ba_reduced, geom->CellSize(),
- IntVect(ncell), IntVect(delta), single_domain_box);
+ IntVect(ncell), IntVect(delta), single_domain_box, v_sigma_sb);
if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
#ifndef WARPX_USE_PSATD
@@ -839,7 +842,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
single_domain_box = is_single_box_domain ? cdomain : Box();
sigba_cp = std::make_unique<MultiSigmaBox>(cba, cdm, grid_cba_reduced, cgeom->CellSize(),
- cncells, cdelta, single_domain_box);
+ cncells, cdelta, single_domain_box, v_sigma_sb);
if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
#ifndef WARPX_USE_PSATD
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 3d010bf02..b980e9403 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -475,6 +475,7 @@ WarpX::InitPML ()
do_multi_J,
do_pml_dive_cleaning, do_pml_divb_cleaning,
guard_cells.ng_FieldSolver.max(),
+ v_particle_pml,
do_pml_Lo[0], do_pml_Hi[0]);
#endif
@@ -506,6 +507,7 @@ WarpX::InitPML ()
do_moving_window, pml_has_particles, do_pml_in_domain,
do_multi_J, do_pml_dive_cleaning, do_pml_divb_cleaning,
guard_cells.ng_FieldSolver.max(),
+ v_particle_pml,
do_pml_Lo[lev], do_pml_Hi[lev]);
}
}
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 3cc018853..31d118733 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -438,6 +438,7 @@ public:
const amrex::MultiFab& getBfield_avg_cp (int lev, int direction) {return *Bfield_avg_cp[lev][direction];}
bool DoPML () const {return do_pml;}
+
#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
const PML_RZ* getPMLRZ() {return pml_rz[0].get();}
#endif
@@ -1279,6 +1280,7 @@ private:
#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
amrex::Vector<std::unique_ptr<PML_RZ> > pml_rz;
#endif
+ amrex::Real v_particle_pml;
amrex::Real moving_window_x = std::numeric_limits<amrex::Real>::max();
amrex::Real current_injection_position = 0;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 31b3bd0a9..9da534ecb 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -853,6 +853,13 @@ WarpX::ReadParameters ()
pp_warpx.query("do_pml_j_damping", do_pml_j_damping);
pp_warpx.query("do_pml_in_domain", do_pml_in_domain);
pp_warpx.query("do_similar_dm_pml", do_similar_dm_pml);
+ // Read `v_particle_pml` in units of the speed of light
+ v_particle_pml = 1._rt;
+ pp_warpx.query("v_particle_pml", v_particle_pml);
+ WARPX_ALWAYS_ASSERT_WITH_MESSAGE(0._rt < v_particle_pml && v_particle_pml <= 1._rt,
+ "Input value for the velocity warpx.v_particle_pml of the macroparticle must be in (0,1] (in units of c).");
+ // Scale by the speed of light
+ v_particle_pml = v_particle_pml * PhysConst::c;
// Default values of WarpX::do_pml_dive_cleaning and WarpX::do_pml_divb_cleaning:
// false for FDTD solver, true for PSATD solver.
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 11:14:42 +0000