Skip to content

More efficient PML BoxArray #2631

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 32 commits into from
Jan 11, 2022
Merged

More efficient PML BoxArray #2631

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
32 commits
Select commit Hold shift + click to select a range
a4ce15c
More efficient PML BoxArray
WeiqunZhang Dec 6, 2021
4b6752e
Update Source/BoundaryConditions/PML.cpp
WeiqunZhang Dec 8, 2021
c1b511a
Apply suggestions from code review
WeiqunZhang Dec 8, 2021
070f35a
Merge branch 'development' into pml_grids
WeiqunZhang Dec 9, 2021
2999304
reset Python_wrappers benchmark
WeiqunZhang Dec 9, 2021
1063dab
Merge branch 'development' into pml_grids
WeiqunZhang Dec 10, 2021
9bf3f6c
Merge branch 'development' into pml_grids
WeiqunZhang Dec 13, 2021
95abac1
Merge branch 'development' into pml_grids
WeiqunZhang Dec 15, 2021
d156028
fix the computation of sigmas for the new BoxArray
WeiqunZhang Dec 17, 2021
0b1939c
Merge branch 'development' into pml_grids
WeiqunZhang Dec 17, 2021
fb2c715
Revert "reset Python_wrappers benchmark"
WeiqunZhang Dec 17, 2021
a0335eb
fix warning
WeiqunZhang Dec 17, 2021
046130b
fix 1d
WeiqunZhang Dec 17, 2021
13747eb
initialize to quiet NaN
WeiqunZhang Dec 18, 2021
5b67892
Merge branch 'development' into pml_grids
WeiqunZhang Dec 18, 2021
425cc51
Merge branch 'development' into pml_grids
WeiqunZhang Dec 22, 2021
0e5cd17
Merge branch 'development' into weiqun_pml_grids
EZoni Jan 5, 2022
4cfacf9
Reset Benchmark: pml_x_psatd
EZoni Jan 6, 2022
1b31579
Reset Benchmark: LaserAccelerationMR
EZoni Jan 6, 2022
0018393
Reset Benchmark: LaserOnFine
EZoni Jan 6, 2022
6926240
Reset Benchmark: PlasmaAccelerationMR
EZoni Jan 6, 2022
b74be74
Reset Benchmark: RefinedInjection
EZoni Jan 7, 2022
8d83686
Reset Benchmark: momentum-conserving-gather
EZoni Jan 7, 2022
794585a
Reset Benchmark: subcyclingMR
EZoni Jan 7, 2022
5db73e9
Reset Benchmark: Langmuir_multi_2d_MR
EZoni Jan 7, 2022
abc271a
Reset Benchmark: Langmuir_multi_2d_MR_psatd
EZoni Jan 7, 2022
eb8c0ec
Reset Benchmark: Python_LaserAccelerationMR
EZoni Jan 7, 2022
fcde666
Reset Benchmark: Python_wrappers
EZoni Jan 7, 2022
b828680
Reset Benchmark: pml_psatd_dive_divb_cleaning
EZoni Jan 7, 2022
1cb1410
Remove an assertion. We will fix it later
WeiqunZhang Jan 11, 2022
d70b5f2
Reset Benchmark: Langmuir_multi_2d_MR_anisotropic
EZoni Jan 11, 2022
a67483a
Reset Benchmark: PEC_field_mr
EZoni Jan 11, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 15 additions & 3 deletions Source/BoundaryConditions/PML.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,11 +220,23 @@ private:
std::unique_ptr<SpectralSolver> spectral_solver_cp;
#endif

static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom,
static amrex::BoxArray MakeBoxArray (const amrex::Box& regular_domain,
const amrex::Geometry& geom,
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());
const amrex::IntVect do_pml_Lo,
const amrex::IntVect do_pml_Hi);

static amrex::BoxArray MakeBoxArray_single (const amrex::Box& regular_domain,
const amrex::BoxArray& grid_ba,
int ncell, const amrex::IntVect do_pml_Lo,
const amrex::IntVect do_pml_Hi);

static amrex::BoxArray MakeBoxArray_multiple (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);

static void CopyToPML (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom);
};
Expand Down
105 changes: 80 additions & 25 deletions Source/BoundaryConditions/PML.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -491,19 +491,22 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g
// Note that this is okay to build pml inside domain for a single patch, or joint patches
// with same [min,max]. But it does not support multiple disjoint refinement patches.
Box domain0 = grid_ba.minimalBox();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (do_pml_Lo[idim]){
domain0.growLo(idim, -ncell);
}
if (do_pml_Hi[idim]){
domain0.growHi(idim, -ncell);
if (do_pml_in_domain) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (do_pml_Lo[idim]){
domain0.growLo(idim, -ncell);
}
if (do_pml_Hi[idim]){
domain0.growHi(idim, -ncell);
}
Comment on lines +554 to +560
Copy link
Member

@EZoni EZoni Dec 8, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One question to make sure I understand what is being done here (it wasn't added in this PR, but still I'd like to ask). My understanding is that this is the case where we want the PML to overlap with the last ncell of the regular domain. Shouldn't we then add ncell (i.e. grow by +ncell, instead of -ncell) to the low index and subtract ncell from the high index (i.e. grow by -ncell, as already done here)? If ncell = 10, the low index would have to be 0 instead of -10 (so -10 + (+ncell) = -10 + (+10) = 0), such that the PML overlaps with the cells 0:9 of the regular domain rather than spanning the cells -10:-1 outside of the domain. While when I read domain0.growLo(idim, -ncell) it seems like we are lowering the low index even more. I guess I'm missing either the starting point or what growLo does precisely. Could you clarify this for me? Thanks!

Copy link
Member Author

@WeiqunZhang WeiqunZhang Dec 8, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

growLo means growing in the direction that points to the lo (i.e., -inf) direction.

}
}
const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0));
const BoxArray grid_ba_reduced = (do_pml_in_domain) ?
BoxArray(grid_ba.boxList().intersect(domain0)) : grid_ba;

const BoxArray& ba = MakeBoxArray(domain0, *geom, grid_ba_reduced, ncell,
do_pml_in_domain, do_pml_Lo, do_pml_Hi);

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.empty()) {
m_ok = false;
return;
Expand Down Expand Up @@ -657,24 +660,27 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g

// assuming that the bounding box around grid_cba is a single patch, and not disjoint patches, similar to fine patch.
amrex::Box domain1 = grid_cba.minimalBox();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (do_pml_Lo[idim]){
// ncell is divided by refinement ratio to ensure that the
// physical width of the PML region is equal in fine and coarse patch
domain1.growLo(idim, -ncell/ref_ratio[idim]);
}
if (do_pml_Hi[idim]){
// ncell is divided by refinement ratio to ensure that the
// physical width of the PML region is equal in fine and coarse patch
domain1.growHi(idim, -ncell/ref_ratio[idim]);
if (do_pml_in_domain) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (do_pml_Lo[idim]){
// ncell is divided by refinement ratio to ensure that the
// physical width of the PML region is equal in fine and coarse patch
domain1.growLo(idim, -ncell/ref_ratio[idim]);
}
if (do_pml_Hi[idim]){
// ncell is divided by refinement ratio to ensure that the
// physical width of the PML region is equal in fine and coarse patch
domain1.growHi(idim, -ncell/ref_ratio[idim]);
}
}
}
const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain1));
const BoxArray grid_cba_reduced = (do_pml_in_domain) ?
BoxArray(grid_cba.boxList().intersect(domain1)) : grid_cba;

// Assuming that refinement ratio is equal in all dimensions
const BoxArray& cba = (do_pml_in_domain) ?
MakeBoxArray(*cgeom, grid_cba_reduced, ncell/ref_ratio[0], 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);
const BoxArray& cba = MakeBoxArray(domain1, *cgeom, grid_cba_reduced,
ncell/ref_ratio[0], do_pml_in_domain, do_pml_Lo,
do_pml_Hi);
DistributionMapping cdm{cba};

pml_E_cp[0] = std::make_unique<MultiFab>(amrex::convert( cba,
Expand Down Expand Up @@ -759,9 +765,58 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g
}

BoxArray
PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba,
PML::MakeBoxArray (const amrex::Box& regular_domain,
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)
{
if (regular_domain.numPts() == grid_ba.numPts()) {
return MakeBoxArray_single(regular_domain, grid_ba, ncell, do_pml_Lo, do_pml_Hi);
} else { // the union of the regular grids is *not* a single rectangular domain
return MakeBoxArray_multiple(geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
}
}

BoxArray
PML::MakeBoxArray_single (const amrex::Box& regular_domain, const amrex::BoxArray& grid_ba,
int ncell, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
{
BoxList bl;
for (int i = 0, N = grid_ba.size(); i < N; ++i) {
Box const& b = grid_ba[i];
for (OrientationIter oit; oit.isValid(); ++oit) { // Iteration over all faces
Orientation ori = oit();
const int idim = ori.coordDir();
bool pml_bndry = false;
if (ori.isLow() && do_pml_Lo[idim]) {
pml_bndry = b.smallEnd(idim) == regular_domain.smallEnd(idim);
} else if (ori.isHigh() && do_pml_Hi[idim]) {
pml_bndry = b.bigEnd(idim) == regular_domain.bigEnd(idim);
}
if (pml_bndry) {
Box bbox = amrex::adjCell(b, ori, ncell);
for (int jdim = 0; jdim < idim; ++jdim) {
if (do_pml_Lo[jdim] &&
bbox.smallEnd(jdim) == regular_domain.smallEnd(jdim)) {
bbox.growLo(jdim, ncell);
}
if (do_pml_Hi[jdim] &&
bbox.bigEnd(jdim) == regular_domain.bigEnd(jdim)) {
bbox.growHi(jdim, ncell);
}
}
bl.push_back(bbox);
}
}
}

return BoxArray(std::move(bl));
}

BoxArray
PML::MakeBoxArray_multiple (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();
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
Expand Down