Skip to content
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

Adding support for pml in domain for multiple patches #4632

Merged
merged 16 commits into from
Feb 28, 2024
Merged

Adding support for pml in domain for multiple patches #4632

Changes from 3 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
d6c8909
this works only on single proc, but lifts the assumption of single pa…
RevathiJambunathan Jan 24, 2024
96cc7fa
fix cbl
RevathiJambunathan Jan 24, 2024
af19b69
building reduced ba and cba
RevathiJambunathan Jan 24, 2024
00c136e
Removing print statements
RevathiJambunathan Jan 24, 2024
4dbfaf0
removing comment
RevathiJambunathan Jan 24, 2024
3cdacfe
assert that box length must be greater than size of pml when its gro…
RevathiJambunathan Jan 24, 2024
528fc6f
same assertion for box length for coarse grid
RevathiJambunathan Jan 24, 2024
82225d6
simplified_list and no need for assert in coarse, but if fine works, …
RevathiJambunathan Jan 24, 2024
d72fd6e
boxlist instead of simplified list
RevathiJambunathan Jan 24, 2024
be9fe48
comments
RevathiJambunathan Feb 26, 2024
3762427
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 26, 2024
1a6b723
CI test for double patch with particles entering and leaving pml insi…
RevathiJambunathan Feb 27, 2024
9dcf20a
checksum for new particles in pml 2d mr test with double patch
RevathiJambunathan Feb 27, 2024
ea0babc
renaming plt to diag1 again
RevathiJambunathan Feb 28, 2024
3afca6d
include level 0 maxfield to capture no j damping effect
RevathiJambunathan Feb 28, 2024
74ce961
removing incorect comment about heavy paritcle in input
RevathiJambunathan Feb 28, 2024
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
71 changes: 42 additions & 29 deletions Source/BoundaryConditions/PML.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -566,34 +566,41 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
// 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)
// minimalBox provides the bounding box around grid_ba for level, lev.
// 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();

BoxArray grid_ba_reduced = grid_ba;
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);
BoxList bl = grid_ba.boxList();
Copy link
Member

Choose a reason for hiding this comment

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

We could also use simplified_list() instead of boxList here. That way, it's less likely to hit the assert failure below.

Copy link
Member Author

Choose a reason for hiding this comment

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

Thank you Weiqun for the suggestions!
Done

for (auto& b : bl) {
Copy link
Member

Choose a reason for hiding this comment

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

Could you add a comment to explain what this for loop does?
(In my understanding, it loops over the box at this level, and whenever a box touches a PML, it reduces the box by ncell, since this is where the surrounding "in-domain" PML will be created.)

Copy link
Member Author

Choose a reason for hiding this comment

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

added comments to clarify

for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (do_pml_Lo[idim]) {
Box const& bb = amrex::adjCellLo(b, idim);
Copy link
Member

@RemiLehe RemiLehe Feb 12, 2024 •
edited
Loading

Choose a reason for hiding this comment

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

Could you add a comment to explains what this code (the call to adjCellLo and the subsequent intersects) does?

Copy link
Member Author

Choose a reason for hiding this comment

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

added comments to clarify

if ( ! grid_ba.intersects(bb) ) {
b.growLo(idim, -ncell);
}
}
if (do_pml_Hi[idim]) {
Box const& bb = amrex::adjCellHi(b, idim);
if ( ! grid_ba.intersects(bb) ) {
b.growHi(idim, -ncell);
}
}
}
}
grid_ba_reduced = BoxArray(std::move(bl));
}
const BoxArray grid_ba_reduced = (do_pml_in_domain) ?
BoxArray(grid_ba.boxList().intersect(domain0)) : grid_ba;

Box const domain0 = grid_ba_reduced.minimalBox();
const bool is_single_box_domain = domain0.numPts() == grid_ba_reduced.numPts();
const BoxArray& ba = MakeBoxArray(is_single_box_domain, domain0, *geom, grid_ba_reduced,
IntVect(ncell), do_pml_in_domain, do_pml_Lo, do_pml_Hi);


if (ba.empty()) {
m_ok = false;
return;
} else {
m_ok = true;
}

// Define the number of guard cells in each direction, for E, B, and F
// Define the number of guard cells in each di;rection, for E, B, and F
auto nge = IntVect(AMREX_D_DECL(2, 2, 2));
auto ngb = IntVect(AMREX_D_DECL(2, 2, 2));
int ngf_int = 0;
Expand Down Expand Up @@ -760,39 +767,43 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri
BoxArray grid_cba = grid_ba;
grid_cba.coarsen(ref_ratio);

// assuming that the bounding box around grid_cba is a single patch, and not disjoint patches, similar to fine patch.
amrex::Box cdomain = grid_cba.minimalBox();
BoxArray grid_cba_reduced = grid_cba;
Copy link
Member

Choose a reason for hiding this comment

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

The code below duplicates the above code.
Maybe it would be good, in a separate follow-up PR, to create a separate function for this.

Copy link
Member Author

Choose a reason for hiding this comment

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

yes. I ll do this in a separate PR

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
cdomain.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
cdomain.growHi(idim, -ncell/ref_ratio[idim]);
BoxList bl = grid_cba.boxList();
Copy link
Member

Choose a reason for hiding this comment

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

Here we can use simplified_list as well.

Copy link
Member Author

Choose a reason for hiding this comment

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

done

Copy link
Member Author

Choose a reason for hiding this comment

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

discussed offline, switched back to boxList due to failing CI

for (auto& b : bl) {
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
if (do_pml_Lo[idim]) {
Box const& bb = amrex::adjCellLo(b, idim);
if ( ! grid_cba.intersects(bb) ) {
b.growLo(idim, -ncell/ref_ratio[idim]);
}
}
if (do_pml_Hi[idim]) {
Box const& bb = amrex::adjCellHi(b, idim);
if ( ! grid_cba.intersects(bb) ) {
b.growHi(idim, -ncell/ref_ratio[idim]);
}
}
}
}
grid_cba_reduced = BoxArray(std::move(bl));
}
const BoxArray grid_cba_reduced = (do_pml_in_domain) ?
BoxArray(grid_cba.boxList().intersect(cdomain)) : grid_cba;
Box const cdomain = grid_cba_reduced.minimalBox();

const IntVect cncells = IntVect(ncell)/ref_ratio;
const IntVect cdelta = IntVect(delta)/ref_ratio;

// Assuming that refinement ratio is equal in all dimensions
const BoxArray& cba = MakeBoxArray(is_single_box_domain, cdomain, *cgeom, grid_cba_reduced,
cncells, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
amrex::Print() << " cba : " << cba << "\n";
DistributionMapping cdm;
if (WarpX::do_similar_dm_pml) {
auto ng_sim = amrex::elemwiseMax(amrex::elemwiseMax(nge, ngb), ngf);
cdm = amrex::MakeSimilarDM(cba, grid_cba_reduced, grid_dm, ng_sim);
} else {
cdm.define(cba);
}

const amrex::BoxArray cba_Ex = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,0).ixType().toIntVect());
const amrex::BoxArray cba_Ey = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,1).ixType().toIntVect());
const amrex::BoxArray cba_Ez = amrex::convert(cba, WarpX::GetInstance().getEfield_cp(1,2).ixType().toIntVect());
Expand Down Expand Up @@ -867,8 +878,10 @@ PML::MakeBoxArray (bool is_single_box_domain, const amrex::Box& regular_domain,
const amrex::IntVect& do_pml_Lo, const amrex::IntVect& do_pml_Hi)
{
if (is_single_box_domain) {
amrex::Print() << " calling MakeBoxArraySingle \n";
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
amrex::Print() << " calling MakeBoxArrayMultiple \n";
return MakeBoxArray_multiple(geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi);
}
}
Expand Down
Loading