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, 15 insertions, 7 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
index e187b9d27..558d4a78a 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
@@ -6,6 +6,7 @@
  */
 #include "SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H"
 #include "SpectralAlgorithms/PsatdAlgorithmRZ.H"
+#include "SpectralAlgorithms/PMLPsatdAlgorithmRZ.H"
 #include "SpectralKSpaceRZ.H"
 #include "SpectralSolverRZ.H"
 #include "Utils/WarpXProfilerWrapper.H"
@@ -22,8 +23,7 @@
  * \param nodal    Whether the solver is applied to a nodal or staggered grid
  * \param dx       Cell size along each dimension
  * \param dt       Time step
- * \param pml      Whether the boxes in which the solver is applied are PML boxes
- *                 PML is not supported.
+ * \param with_pml Whether PML boundary will be used
  */
 SpectralSolverRZ::SpectralSolverRZ (const int lev,
                                     amrex::BoxArray const & realspace_ba,
@@ -32,6 +32,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
                                     int const norder_z, bool const nodal,
                                     const amrex::Vector<amrex::Real>& v_galilean,
                                     amrex::RealVect const dx, amrex::Real const dt,
+                                    bool const with_pml,
                                     bool const update_with_rho,
                                     const bool fft_do_time_averaging,
                                     const bool do_multi_J,
@@ -45,13 +46,16 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
     //   the spectral space corresponding to each box in `realspace_ba`,
     //   as well as the value of the corresponding k coordinates.
 
-    const bool pml = false;
+    const bool is_pml = false;
     m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging,
-                                          do_multi_J, dive_cleaning, divb_cleaning, pml);
+                                          do_multi_J, dive_cleaning, divb_cleaning, is_pml, with_pml);
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
-    //   PML is not supported.
+    if (with_pml) {
+            PML_algorithm = std::make_unique<PMLPsatdAlgorithmRZ>(
+                k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt);
+    }
     if (v_galilean[2] == 0) {
          // v_galilean is 0: use standard PSATD algorithm
         algorithm = std::make_unique<PsatdAlgorithmRZ>(
@@ -117,12 +121,16 @@ SpectralSolverRZ::BackwardTransform (const int lev,
 
 /* \brief Update the fields in spectral space, over one timestep */
 void
-SpectralSolverRZ::pushSpectralFields () {
+SpectralSolverRZ::pushSpectralFields (const bool doing_pml) {
     WARPX_PROFILE("SpectralSolverRZ::pushSpectralFields");
     // Virtual function: the actual function used here depends
     // on the sub-class of `SpectralBaseAlgorithm` that was
     // initialized in the constructor of `SpectralSolverRZ`
-    algorithm->pushSpectralFields(field_data);
+    if (doing_pml) {
+        PML_algorithm->pushSpectralFields(field_data);
+    } else {
+        algorithm->pushSpectralFields(field_data);
+    }
 }
 
 /**