Skip to content

Commit

Permalink
add pml parameters for ref patch : abcinpml, cubicsigma, dampingstrength
Browse files Browse the repository at this point in the history
  • Loading branch information
@RevathiJambunathan
RevathiJambunathan committed Mar 2, 2024
1 parent 1efba85 commit 1f04d31
Show file tree
Hide file tree
Showing 7 changed files with 220 additions and 57 deletions.
22 changes: 15 additions & 7 deletions Source/BoundaryConditions/PML.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,15 +44,18 @@ 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, amrex::Real v_sigma);
const amrex::Box& regdomain, amrex::Real v_sigma,
const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength);

void define_single (const amrex::Box& regdomain, const amrex::IntVect& ncell,
const amrex::Array<amrex::Real,AMREX_SPACEDIM>& fac,
amrex::Real v_sigma);
amrex::Real v_sigma,
const int do_cubic_sigma_pml);
void define_multiple (const amrex::Box& box, const amrex::BoxArray& grids,
const amrex::IntVect& ncell,
const amrex::Array<amrex::Real,AMREX_SPACEDIM>& fac,
amrex::Real v_sigma);
amrex::Real v_sigma,
const int do_cubic_sigma_pml);

void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt);
void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt);
Expand All @@ -79,8 +82,9 @@ 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, 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) {}
const amrex::Box& regular_domain, const amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength)
: m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta), m_regdomain(regular_domain), m_v_sigma_sb(v_sigma_sb), m_do_cubic_sigma_pml(do_cubic_sigma_pml), m_pml_damping_strength(pml_damping_strength) {}
~SigmaBoxFactory () override = default;

SigmaBoxFactory (const SigmaBoxFactory&) = default;
Expand All @@ -93,7 +97,7 @@ public:
[[nodiscard]] 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, m_v_sigma_sb);
return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta, m_regdomain, m_v_sigma_sb, m_do_cubic_sigma_pml, m_pml_damping_strength);
}

void destroy (SigmaBox* fab) const final
Expand All @@ -114,6 +118,8 @@ private:
amrex::IntVect m_delta;
amrex::Box m_regdomain;
amrex::Real m_v_sigma_sb;
int m_do_cubic_sigma_pml;
amrex::Real m_pml_damping_strength;
};

class MultiSigmaBox
Expand All @@ -123,7 +129,8 @@ 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, amrex::Real v_sigma_sb);
const amrex::Box& regular_domain, amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength);
void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt);
void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt);
private:
Expand All @@ -147,6 +154,7 @@ public:
const amrex::IntVect& fill_guards_fields,
const amrex::IntVect& fill_guards_current,
int max_guard_EB, amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength,
amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(),
amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector());

Expand Down
86 changes: 56 additions & 30 deletions Source/BoundaryConditions/PML.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,11 +61,11 @@ namespace
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 amrex::Real v_sigma)
const amrex::Real v_sigma, const int do_cubic_sigma_pml)
{
const int slo = sigma.m_lo;
const int sslo = sigma_star.m_lo;

int sigma_exp = (do_cubic_sigma_pml==1) ? 3 : 2;
const int N = ohi+1-olo+1;
Real* p_sigma = sigma.data();
Real* p_sigma_cumsum = sigma_cumsum.data();
Expand All @@ -75,25 +75,34 @@ namespace
{
i += olo;
Real offset = static_cast<Real>(glo-i);
p_sigma[i-slo] = fac*(offset*offset);
amrex::Real offset_exp = 1._rt;
for (int iexp = 0; iexp < sigma_exp; ++iexp){
offset_exp *= offset;
}
p_sigma[i-slo] = fac*(offset_exp);
// 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)/v_sigma;
p_sigma_cumsum[i-slo] = (fac*(offset*offset_exp)/(sigma_exp+1._rt))/v_sigma;
if (i <= ohi+1) {
offset = static_cast<Real>(glo-i) - 0.5_rt;
p_sigma_star[i-sslo] = fac*(offset*offset);
offset_exp = 1._rt;
for (int iexp = 0; iexp < sigma_exp; ++iexp){
offset_exp *= offset;
}
p_sigma_star[i-sslo] = fac*(offset_exp);
// 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)/v_sigma;
p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset_exp)/(sigma_exp + 1._rt))/v_sigma;
}
});
}

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 amrex::Real v_sigma)
const amrex::Real v_sigma, const int do_cubic_sigma_pml)
{
const int slo = sigma.m_lo;
const int sslo = sigma_star.m_lo;
int sigma_exp = (do_cubic_sigma_pml==1) ? 3 : 2;

const int N = ohi+1-olo+1;
Real* p_sigma = sigma.data();
Expand All @@ -104,12 +113,20 @@ 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)/v_sigma;
amrex::Real offset_exp = 1._rt;
for (int iexp = 0; iexp < sigma_exp; ++iexp){
offset_exp *= offset;
}
p_sigma[i-slo] = fac*(offset_exp);
p_sigma_cumsum[i-slo] = (fac*(offset*offset_exp)/(sigma_exp+1._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)/v_sigma;
offset_exp = 1._rt;
for (int iexp = 0; iexp < sigma_exp; ++iexp){
offset_exp *= offset;
}
p_sigma_star[i-sslo] = fac*(offset_exp);
p_sigma_star_cumsum[i-sslo] = (fac*(offset*offset_exp)/(sigma_exp + 1._rt))/v_sigma;
}
});
}
Expand Down Expand Up @@ -143,7 +160,8 @@ namespace


SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const IntVect& ncell,
const IntVect& delta, const amrex::Box& regdomain, const amrex::Real v_sigma_sb)
const IntVect& delta, const amrex::Box& regdomain, const amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength)
{
BL_ASSERT(box.cellCentered());

Expand Down Expand Up @@ -179,22 +197,27 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, const
sigma_star_cumsum_fac[idim].m_lo = lo[idim];
sigma_star_cumsum_fac[idim].m_hi = hi[idim]+1;
}

int sigma_exp = (do_cubic_sigma_pml==1) ? 3 : 2;
Array<Real,AMREX_SPACEDIM> fac;
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
fac[idim] = 4.0_rt*PhysConst::c/(dx[idim]*static_cast<Real>(delta[idim]*delta[idim]));
amrex::Real delta_exp = 1._rt;
for (int iexp = 0; iexp < sigma_exp; ++iexp) {
delta_exp *= static_cast<amrex::Real>(delta[idim]);
}
fac[idim] = pml_damping_strength*PhysConst::c/(dx[idim]*delta_exp);
}

if (regdomain.ok()) { // The union of the regular grids is a single box
define_single(regdomain, ncell, fac, v_sigma_sb);
define_single(regdomain, ncell, fac, v_sigma_sb, do_cubic_sigma_pml);
} else {
define_multiple(box, grids, ncell, fac, v_sigma_sb);
define_multiple(box, grids, ncell, fac, v_sigma_sb, do_cubic_sigma_pml);
}
}

void SigmaBox::define_single (const Box& regdomain, const IntVect& ncell,
const Array<Real,AMREX_SPACEDIM>& fac,
const amrex::Real v_sigma_sb)
const amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml)
{
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
const int slo = sigma[idim].lo();
Expand All @@ -208,7 +231,8 @@ 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], v_sigma_sb);
olo, ohi, dlo, fac[idim], v_sigma_sb,
do_cubic_sigma_pml);
}

#if (AMREX_SPACEDIM != 1)
Expand All @@ -228,15 +252,16 @@ 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], v_sigma_sb);
olo, ohi, dhi, fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}
}

amrex::Gpu::streamSynchronize();
}

void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, const IntVect& ncell,
const Array<Real,AMREX_SPACEDIM>& fac, const amrex::Real v_sigma_sb)
const Array<Real,AMREX_SPACEDIM>& fac, const amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml)
{
const std::vector<std::pair<int,Box> >& isects = grids.intersections(box, false, ncell);

Expand Down Expand Up @@ -306,7 +331,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], v_sigma_sb);
grid_box.smallEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}

Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]);
Expand All @@ -319,7 +344,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], v_sigma_sb);
grid_box.bigEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}

WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
Expand Down Expand Up @@ -355,7 +380,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], v_sigma_sb);
grid_box.smallEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}

Box hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]);
Expand All @@ -364,7 +389,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], v_sigma_sb);
grid_box.bigEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}

WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
Expand Down Expand Up @@ -405,7 +430,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], v_sigma_sb);
grid_box.smallEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}

const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell[idim]);
Expand All @@ -414,7 +439,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], v_sigma_sb);
grid_box.bigEnd(idim), fac[idim], v_sigma_sb, do_cubic_sigma_pml);
}

WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
Expand Down Expand Up @@ -505,9 +530,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::Real v_sigma_sb)
const amrex::Box& regular_domain, const amrex::Real v_sigma_sb, const int do_cubic_sigma_pml, const amrex::Real pml_damping_strength)
: FabArray<SigmaBox>(ba,dm,1,0,MFInfo(),
SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain, v_sigma_sb))
SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain, v_sigma_sb,do_cubic_sigma_pml,pml_damping_strength))
{}

void
Expand Down Expand Up @@ -552,6 +577,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
const amrex::IntVect& fill_guards_fields,
const amrex::IntVect& fill_guards_current,
int max_guard_EB, const amrex::Real v_sigma_sb,
const int do_cubic_sigma_pml, amrex::Real pml_damping_strength,
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),
Expand Down Expand Up @@ -737,7 +763,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, v_sigma_sb);
IntVect(ncell), IntVect(delta), single_domain_box, v_sigma_sb, do_cubic_sigma_pml, pml_damping_strength);

if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
#ifndef WARPX_USE_PSATD
Expand Down Expand Up @@ -851,7 +877,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, v_sigma_sb);
cncells, cdelta, single_domain_box, v_sigma_sb, do_cubic_sigma_pml, pml_damping_strength);

if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
#ifndef WARPX_USE_PSATD
Expand Down
14 changes: 7 additions & 7 deletions Source/BoundaryConditions/WarpXEvolvePML.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ WarpX::DampPML_Cartesian (const int lev, PatchType patch_type)
{
const bool dive_cleaning = WarpX::do_pml_dive_cleaning;
const bool divb_cleaning = WarpX::do_pml_divb_cleaning;

const bool abc_in_pml = WarpX::do_abc_in_pml;
if (pml[lev]->ok())
{
const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
Expand Down Expand Up @@ -154,39 +154,39 @@ WarpX::DampPML_Cartesian (const int lev, PatchType patch_type)

warpx_damp_pml_ex(i, j, k, pml_Exfab, Ex_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
dive_cleaning);
dive_cleaning, abc_in_pml);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {

warpx_damp_pml_ey(i, j, k, pml_Eyfab, Ey_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
dive_cleaning);
dive_cleaning, abc_in_pml);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {

warpx_damp_pml_ez(i, j, k, pml_Ezfab, Ez_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
dive_cleaning);
dive_cleaning, abc_in_pml);
});

amrex::ParallelFor(tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {

warpx_damp_pml_bx(i, j, k, pml_Bxfab, Bx_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
divb_cleaning);
divb_cleaning, abc_in_pml);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {

warpx_damp_pml_by(i, j, k, pml_Byfab, By_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
divb_cleaning);
divb_cleaning, abc_in_pml);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {

warpx_damp_pml_bz(i, j, k, pml_Bzfab, Bz_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
divb_cleaning);
divb_cleaning, abc_in_pml);
});

// For warpx_damp_pml_F(), mfi.nodaltilebox is used in the ParallelFor loop and here we
Expand Down
Loading

0 comments on commit 1f04d31

Please sign in to comment.