Add support of do_pml_in_domain option

Source/BoundaryConditions/PML_RZ.H

@@ -37,7 +37,7 @@ class PML_RZ
 {
 public:
     PML_RZ (const int lev, const amrex::BoxArray& grid_ba, const amrex::DistributionMapping& grid_dm,
-            const amrex::Geometry* geom, const int ncell);
+            const amrex::Geometry* geom, const int ncell, const int do_pml_in_domain);
 
     void ApplyDamping(const std::array<amrex::MultiFab*,3>& E_fp,
                       const std::array<amrex::MultiFab*,3>& B_fp,
@@ -65,6 +65,7 @@ private:
     bool m_ok;
 
     const int m_ncell;
+    const int m_do_pml_in_domain;
     const amrex::Geometry* m_geom;
 
     // Only contains Er and Et, and Br and Bt

Source/BoundaryConditions/PML_RZ.cpp

@@ -49,8 +49,9 @@
 using namespace amrex;
 
 PML_RZ::PML_RZ (const int lev, const BoxArray& grid_ba, const DistributionMapping& grid_dm,
-                const Geometry* geom, const int ncell)
+                const Geometry* geom, const int ncell, const int do_pml_in_domain)
     : m_ncell(ncell),
+      m_do_pml_in_domain(do_pml_in_domain),
       m_geom(geom)
 {
 
@@ -103,16 +104,18 @@ PML_RZ::ApplyDamping(const std::array<amrex::MultiFab*,3>& E_fp,
         // They are all the same, cell centered
         amrex::Box tilebox = amrex::convert((*E_fp[1])[mfi].box(), E_fp[0]->ixType().toIntVect());
 
-        // Get upper radius of tilebox
-        const int tilebox_bigEnd_r = tilebox.bigEnd(0);
-
         // Box for the whole simulation domain
         amrex::Box const& domain = m_geom->Domain();
         const int nr_domain = domain.bigEnd(0);
 
         // Set tilebox to only include the upper radial cells
         const int nr_damp = m_ncell;
-        const int nr_damp_min = nr_domain - nr_damp;
+        int nr_damp_min;
+        if (m_do_pml_in_domain) {
+            nr_damp_min = nr_domain - nr_damp;
+        } else {
+            nr_damp_min = nr_domain;
+        }
         tilebox.setSmall(0, nr_damp_min + 1);
 
         amrex::ParallelFor( tilebox, E_fp[0]->nComp(),

Source/Initialization/WarpXInitData.cpp

@@ -242,7 +242,7 @@ WarpX::InitPML ()
 
 #ifdef WARPX_DIM_RZ
         do_pml_Lo_corrected[0] = 0; // no PML at r=0, in cylindrical geometry
-        pml_rz[0] = std::make_unique<PML_RZ>(0, boxArray(0), DistributionMap(0), &Geom(0), pml_ncell);
+        pml_rz[0] = std::make_unique<PML_RZ>(0, boxArray(0), DistributionMap(0), &Geom(0), pml_ncell, do_pml_in_domain);
 #else
         pml[0] = std::make_unique<PML>(0, boxArray(0), DistributionMap(0), &Geom(0), nullptr,
                              pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(),

Source/Parallelization/GuardCellManager.H

@@ -36,6 +36,9 @@ public:
      * \param nci_corr_stencil stencil of NCI corrector
      * \param maxwell_solver_id if of Maxwell solver
      * \param max_level max level of the simulation
+     * \param do_pml whether pml is turned on (only used by RZ PSATD)
+     * \param do_pml_in_domain whether pml is done in the domain (only used by RZ PSATD)
+     * \param pml_ncell number of cells on the pml layer (only used by RZ PSATD)
      */
     void Init(
         const amrex::Real dt,
@@ -53,7 +56,10 @@ public:
         const amrex::Array<amrex::Real,3> v_galilean,
         const amrex::Array<amrex::Real,3> v_comoving,
         const bool safe_guard_cells,
-        const int do_electrostatic);
+        const int do_electrostatic,
+        const bool do_pml,
+        const int do_pml_in_domain,
+        const int pml_ncell);
 
     // Guard cells allocated for MultiFabs E and B
     amrex::IntVect ng_alloc_EB = amrex::IntVect::TheZeroVector();

Source/Parallelization/GuardCellManager.cpp

@@ -45,7 +45,10 @@ guardCellManager::Init (
     const amrex::Array<amrex::Real,3> v_galilean,
     const amrex::Array<amrex::Real,3> v_comoving,
     const bool safe_guard_cells,
-    const int do_electrostatic)
+    const int do_electrostatic,
+    const bool do_pml,
+    const int do_pml_in_domain,
+    const int pml_ncell)
 {
     // When using subcycling, the particles on the fine level perform two pushes
     // before being redistributed ; therefore, we need one extra guard cell
@@ -174,6 +177,16 @@ guardCellManager::Init (
         IntVect ngFFT = IntVect(ngFFt_x, ngFFt_z);
 #endif
 
+#ifdef WARPX_DIM_RZ
+        if (do_pml) {
+            if (!do_pml_in_domain) {
+                ngFFT[0] = std::max(ngFFT[0], pml_ncell);
+            }
+        }
+#else
+       amrex::ignore_unused(do_pml, do_pml_in_domain, pml_ncell);
+#endif
+
         // All boxes should have the same number of guard cells, to avoid temporary parallel copies:
         // thus we take the maximum of the required number of guard cells over all available fields.
         for (int i_dim = 0; i_dim < AMREX_SPACEDIM; i_dim++) {

Source/WarpX.cpp

@@ -742,6 +742,8 @@ WarpX::ReadParameters ()
 #ifdef WARPX_DIM_RZ
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE( isAnyBoundaryPML() == false || maxwell_solver_id == MaxwellSolverAlgo::PSATD,
             "PML are not implemented in RZ geometry with FDTD; please set a different boundary condition using boundary.field_lo and boundary.field_hi.");
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( field_boundary_lo[1] != FieldBoundaryType::PML && field_boundary_hi[1] != FieldBoundaryType::PML,
+            "PML are not implemented in RZ geometry along z; please set a different boundary condition using boundary.field_lo and boundary.field_hi.");
 #endif
 
         if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){
@@ -1338,7 +1340,10 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d
         WarpX::m_v_galilean,
         WarpX::m_v_comoving,
         safe_guard_cells,
-        WarpX::do_electrostatic);
+        WarpX::do_electrostatic,
+        WarpX::isAnyBoundaryPML(),
+        WarpX::do_pml_in_domain,
+        WarpX::pml_ncell);
 
     if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) {
         n_current_deposition_buffer = 1;
@@ -1559,6 +1564,10 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         if ( fft_periodic_single_box == false ) {
             realspace_ba.grow(1, ngE[1]); // add guard cells only in z
         }
+        if (field_boundary_hi[0] == FieldBoundaryType::PML && !do_pml_in_domain) {
+            // Extend region that is solved for the PML
+            realspace_ba.growHi(0, pml_ncell);
+        }
         AllocLevelSpectralSolverRZ(spectral_solver_fp,
                                    lev,
                                    realspace_ba,