aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/BoundaryConditions/PML.cpp127
1 files changed, 82 insertions, 45 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 7d1906525..03d9dd9d8 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -15,7 +15,8 @@ using namespace amrex;
namespace
{
- static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum,
+ 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);
@@ -29,6 +30,7 @@ namespace
{
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;
}
@@ -36,11 +38,13 @@ namespace
{
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_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum,
+ 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);
@@ -63,7 +67,9 @@ namespace
}
}
- static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum, 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);
@@ -86,29 +92,29 @@ 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_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 [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 [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];
}
@@ -178,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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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);
@@ -188,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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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()) {
@@ -202,7 +212,8 @@ 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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], overlap);
+ 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");
@@ -216,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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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()) {
@@ -240,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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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");
}
@@ -253,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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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_cumsum[idim], sigma_star[idim], sigma_star_cumsum[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()) {
@@ -341,9 +361,9 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt)
PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
const Geometry* geom, const Geometry* cgeom,
int ncell, int delta, int ref_ratio,
- #ifdef WARPX_USE_PSATD
- Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
- #endif
+#ifdef WARPX_USE_PSATD
+ 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)
@@ -351,6 +371,11 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
m_cgeom(cgeom)
{
+ // 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)) {
@@ -365,7 +390,9 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
}
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_Lo, do_pml_Hi) : MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi);
+ const BoxArray& ba = (do_pml_in_domain)?
+ MakeBoxArray(*geom, grid_ba_reduced, ncell, do_pml_Lo, do_pml_Hi) :
+ MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi);
if (ba.size() == 0) {
m_ok = false;
@@ -460,7 +487,9 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
grid_cba.coarsen(ref_ratio);
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_Lo, do_pml_Hi) : MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi);
+ const BoxArray& cba = (do_pml_in_domain) ?
+ MakeBoxArray(*cgeom, grid_cba_reduced, ncell, do_pml_Lo, do_pml_Hi) :
+ MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi);
DistributionMapping cdm{cba};
@@ -493,8 +522,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm,
if (do_pml_in_domain){
sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta));
- }
- else {
+ } else {
sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta));
}
@@ -530,8 +558,13 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
{
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);
@@ -638,7 +671,8 @@ PML::GetF_cp ()
void
PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp,
- const std::array<amrex::MultiFab*,3>& B_cp, int do_pml_in_domain)
+ const std::array<amrex::MultiFab*,3>& B_cp,
+ int do_pml_in_domain)
{
ExchangeB(PatchType::fine, B_fp, do_pml_in_domain);
ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain);
@@ -646,7 +680,8 @@ PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp,
void
PML::ExchangeB (PatchType patch_type,
- const std::array<amrex::MultiFab*,3>& Bp, int do_pml_in_domain)
+ const std::array<amrex::MultiFab*,3>& Bp,
+ int do_pml_in_domain)
{
if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0])
{
@@ -664,7 +699,8 @@ PML::ExchangeB (PatchType patch_type,
void
PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
- const std::array<amrex::MultiFab*,3>& E_cp, int do_pml_in_domain)
+ const std::array<amrex::MultiFab*,3>& E_cp,
+ int do_pml_in_domain)
{
ExchangeE(PatchType::fine, E_fp, do_pml_in_domain);
ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain);
@@ -672,7 +708,8 @@ PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
void
PML::ExchangeE (PatchType patch_type,
- const std::array<amrex::MultiFab*,3>& Ep, int do_pml_in_domain)
+ const std::array<amrex::MultiFab*,3>& Ep,
+ int do_pml_in_domain)
{
if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0])
{
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 02:31:24 +0000