Implement PML for the outer RZ boundary with PSATD (#2211)

* Initial version of RZ PSATD PML BCs

* Cleaned up some bugs

* Add support of do_pml_in_domain option

* Cleaned up stuff for building

* Fix PMLPsatdAlgorithm macro

* Removed unneeded variable from SpectralSolverRZ

* Change length 3 arrays to length 2 (for 2D)

* Cleanup around DampPML

* Added more checks of pml[lev]

* Added CI test for RZ PML

* Added code to update the corner guard cells

* Further updates

* Added CI test

* Fixed EOL space

* Updated CI benchmarks, removing round off fields

* Changes to CI missed on previous commit

* Various fixes for clean up

* More fixes for clean up

* Further cleanup

* Updated benchmark

* Fixed benchmarks file

* Minor cleanup

* Added round off benchmark values

* Fixed testname in analysis_pml_psatd_rz.py

* Update comment in analysis file

* Put pml_rz code in RZ and PSATD macro blocks

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Add geometry.dims input to CI test input file, inputs_rz

* Cleanup to match recent changes

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

1 files changed, 39 insertions, 2 deletions

diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp
index d844c2317..34087ab27 100644
--- a/Source/Utils/WarpXMovingWindow.cpp
+++ b/Source/Utils/WarpXMovingWindow.cpp
@@ -9,6 +9,9 @@
 #include "WarpX.H"
 
 #include "BoundaryConditions/PML.H"
+#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
+#   include "BoundaryConditions/PML_RZ.H"
+#endif
 #include "Particles/MultiParticleContainer.H"
 #include "Parallelization/WarpXCommUtil.H"
 #include "Utils/WarpXConst.H"
@@ -179,13 +182,20 @@ WarpX::MoveWindow (const int step, bool move_j)
             if (move_j) {
                 shiftMF(*current_fp[lev][dim], geom[lev], num_shift, dir);
             }
-            if (do_pml && pml[lev]->ok()) {
+            if (pml[lev] && pml[lev]->ok()) {
                 const std::array<amrex::MultiFab*, 3>& pml_B = pml[lev]->GetB_fp();
                 const std::array<amrex::MultiFab*, 3>& pml_E = pml[lev]->GetE_fp();
                 shiftMF(*pml_B[dim], geom[lev], num_shift, dir);
                 shiftMF(*pml_E[dim], geom[lev], num_shift, dir);
             }
-
+#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
+            if (pml_rz[lev] && dim < 2) {
+                const std::array<amrex::MultiFab*, 2>& pml_rz_B = pml_rz[lev]->GetB_fp();
+                const std::array<amrex::MultiFab*, 2>& pml_rz_E = pml_rz[lev]->GetE_fp();
+                shiftMF(*pml_rz_B[dim], geom[lev], num_shift, dir);
+                shiftMF(*pml_rz_E[dim], geom[lev], num_shift, dir);
+            }
+#endif
             if (lev > 0) {
                 // coarse grid
                 shiftMF(*Bfield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, B_external_grid[dim], use_Bparser, Bfield_parser);
@@ -407,6 +417,33 @@ WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift,
             dstfab(i,j,k,n) = srcfab(i+shift.x,j+shift.y,k+shift.z,n);
         })
     }
+
+#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
+    if (WarpX::GetInstance().getPMLRZ()) {
+        // This does the exchange of data in the corner guard cells, the cells that are in the
+        // guard region both radially and longitudinally. These are the PML cells in the overlapping
+        // longitudinal region. FillBoundary normally does not update these cells.
+        // This update is needed so that the cells at the end of the FABs are updated appropriately
+        // with the data shifted from the nieghboring FAB. Without this update, the RZ PML becomes
+        // unstable with the moving grid.
+        // This code creates a temporary MultiFab using a BoxList where the radial size of all of
+        // its boxes is increased so that the radial guard cells are included in the boxes valid domain.
+        // The temporary MultiFab is setup to refer to the data of the original Multifab (this can
+        // be done since the shape of the data is all the same, just the indexing is different).
+        amrex::BoxList bl;
+        for (int i = 0, N=ba.size(); i < N; ++i) {
+            bl.push_back(amrex::grow(ba[i], 0, mf.nGrowVect()[0]));
+        }
+        amrex::BoxArray rba(std::move(bl));
+        amrex::MultiFab rmf(rba, dm, mf.nComp(), IntVect(0,mf.nGrowVect()[1]), MFInfo().SetAlloc(false));
+
+        for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) {
+            rmf.setFab(mfi, FArrayBox(mf[mfi], amrex::make_alias, 0, mf.nComp()));
+        }
+        rmf.FillBoundary(false);
+    }
+#endif
+
 }
 
 void