Fixes to RZ PSATD (#1945)

* For RZ, changed the sign of the density corrections near the axis

* Further fixes for deposition correction near axis

* Yet one more sign fix for charge density

* For RZ spectral solver, filled in the guard cells below the radial axis

* Fix white space at end of line

* In RZ spectral backtransform, ensure box is valid

* For RZ inverse volume scaling, fixed use of nGrow to use nGrowVect

* Temporary fix adding damped cells in the domain interior

* Bug fix for RZ PSATD scalar backward transform

* Fixes for damping of the fields in the z-guards

* Bug fix in DampFieldsInGuards

* Bug fix in DampFieldsInGuards (for tiling)

* Added warpx_amr_check_input input parameter

* Removed unneeded damp and zero_in_domain input

* Removed damping related code from picmi

* Improved some comments in code copying field to the radial guard cells

* Update Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp

Simplify the expression for the sign

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Updated benchmarks

* Updated tolerance for Langmuir analysis script

* Updated CI test galilean_rz_psatd_current_correction

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

1 files changed, 41 insertions, 30 deletions

diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 1e375f0eb..15ce940d6 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -751,7 +751,7 @@ WarpX::DampFieldsInGuards(std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield
 void
 WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, MultiFab* Jz, int lev)
 {
-    const long ngJ = Jx->nGrow();
+    const amrex::IntVect ngJ = Jx->nGrowVect();
     const std::array<Real,3>& dx = WarpX::CellSize(lev);
     const Real dr = dx[0];
 
@@ -793,16 +793,16 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
         // Grow the tileboxes to include the guard cells, except for the
         // guard cells at negative radius.
         if (rmin > 0.) {
-           tbr.growLo(0, ngJ);
-           tbt.growLo(0, ngJ);
-           tbz.growLo(0, ngJ);
+           tbr.growLo(0, ngJ[0]);
+           tbt.growLo(0, ngJ[0]);
+           tbz.growLo(0, ngJ[0]);
         }
-        tbr.growHi(0, ngJ);
-        tbt.growHi(0, ngJ);
-        tbz.growHi(0, ngJ);
-        tbr.grow(1, ngJ);
-        tbt.grow(1, ngJ);
-        tbz.grow(1, ngJ);
+        tbr.growHi(0, ngJ[0]);
+        tbt.growHi(0, ngJ[0]);
+        tbz.growHi(0, ngJ[0]);
+        tbr.grow(1, ngJ[1]);
+        tbt.grow(1, ngJ[1]);
+        tbz.grow(1, ngJ[1]);
 
         // Rescale current in r-z mode since the inverse volume factor was not
         // included in the current deposition.
@@ -812,7 +812,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
             // Wrap the current density deposited in the guard cells around
             // to the cells above the axis.
             // Note that Jr(i==0) is at 1/2 dr.
-            if (rmin == 0. && 0 <= i && i < ngJ) {
+            if (rmin == 0. && 0 <= i && i < ngJ[0]) {
                 Jr_arr(i,j,0,0) -= Jr_arr(-1-i,j,0,0);
             }
             // Apply the inverse volume scaling
@@ -825,9 +825,9 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
                 // Wrap the current density deposited in the guard cells around
                 // to the cells above the axis.
                 // Note that Jr(i==0) is at 1/2 dr.
-                if (rmin == 0. && 0 <= i && i < ngJ) {
-                    Jr_arr(i,j,0,2*imode-1) -= Jr_arr(-1-i,j,0,2*imode-1);
-                    Jr_arr(i,j,0,2*imode) -= Jr_arr(-1-i,j,0,2*imode);
+                if (rmin == 0. && 0 <= i && i < ngJ[0]) {
+                    Jr_arr(i,j,0,2*imode-1) += std::pow(-1, imode+1)*Jr_arr(-1-i,j,0,2*imode-1);
+                    Jr_arr(i,j,0,2*imode) += std::pow(-1, imode+1)*Jr_arr(-1-i,j,0,2*imode);
                 }
                 // Apply the inverse volume scaling
                 // Since Jr is never node centered in r, no need for distinction
@@ -842,7 +842,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
             // to the cells above the axis.
             // If Jt is node centered, Jt[0] is located on the boundary.
             // If Jt is cell centered, Jt[0] is at 1/2 dr.
-            if (rmin == 0. && 1-ishift_t <= i && i <= ngJ-ishift_t) {
+            if (rmin == 0. && 1-ishift_t <= i && i <= ngJ[0]-ishift_t) {
                 Jt_arr(i,j,0,0) -= Jt_arr(-ishift_t-i,j,0,0);
             }
 
@@ -858,9 +858,9 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
             for (int imode=1 ; imode < nmodes ; imode++) {
                 // Wrap the current density deposited in the guard cells around
                 // to the cells above the axis.
-                if (rmin == 0. && 1-ishift_t <= i && i <= ngJ-ishift_t) {
-                    Jt_arr(i,j,0,2*imode-1) -= Jt_arr(-ishift_t-i,j,0,2*imode-1);
-                    Jt_arr(i,j,0,2*imode) -= Jt_arr(-ishift_t-i,j,0,2*imode);
+                if (rmin == 0. && 1-ishift_t <= i && i <= ngJ[0]-ishift_t) {
+                    Jt_arr(i,j,0,2*imode-1) += std::pow(-1, imode+1)*Jt_arr(-ishift_t-i,j,0,2*imode-1);
+                    Jt_arr(i,j,0,2*imode) += std::pow(-1, imode+1)*Jt_arr(-ishift_t-i,j,0,2*imode);
                 }
 
                 // Apply the inverse volume scaling
@@ -880,8 +880,8 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
             // to the cells above the axis.
             // If Jz is node centered, Jt[0] is located on the boundary.
             // If Jz is cell centered, Jt[0] is at 1/2 dr.
-            if (rmin == 0. && 1-ishift_z <= i && i <= ngJ-ishift_z) {
-                Jz_arr(i,j,0,0) -= Jz_arr(-ishift_z-i,j,0,0);
+            if (rmin == 0. && 1-ishift_z <= i && i <= ngJ[0]-ishift_z) {
+                Jz_arr(i,j,0,0) += Jz_arr(-ishift_z-i,j,0,0);
             }
 
             // Apply the inverse volume scaling
@@ -896,9 +896,9 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
             for (int imode=1 ; imode < nmodes ; imode++) {
                 // Wrap the current density deposited in the guard cells around
                 // to the cells above the axis.
-                if (rmin == 0. && 1-ishift_z <= i && i <= ngJ-ishift_z) {
-                    Jz_arr(i,j,0,2*imode-1) -= Jz_arr(-ishift_z-i,j,0,2*imode-1);
-                    Jz_arr(i,j,0,2*imode) -= Jz_arr(-ishift_z-i,j,0,2*imode);
+                if (rmin == 0. && 1-ishift_z <= i && i <= ngJ[0]-ishift_z) {
+                    Jz_arr(i,j,0,2*imode-1) -= std::pow(-1, imode+1)*Jz_arr(-ishift_z-i,j,0,2*imode-1);
+                    Jz_arr(i,j,0,2*imode) -= std::pow(-1, imode+1)*Jz_arr(-ishift_z-i,j,0,2*imode);
                 }
 
                 // Apply the inverse volume scaling
@@ -919,7 +919,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu
 void
 WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev)
 {
-    const long ngRho = Rho->nGrow();
+    const amrex::IntVect ngRho = Rho->nGrowVect();
     const std::array<Real,3>& dx = WarpX::CellSize(lev);
     const Real dr = dx[0];
 
@@ -949,23 +949,34 @@ WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev)
         // Grow the tilebox to include the guard cells, except for the
         // guard cells at negative radius.
         if (rmin > 0.) {
-           tb.growLo(0, ngRho);
+           tb.growLo(0, ngRho[0]);
         }
-        tb.growHi(0, ngRho);
-        tb.grow(1, ngRho);
+        tb.growHi(0, ngRho[0]);
+        tb.grow(1, ngRho[1]);
 
         // Rescale charge in r-z mode since the inverse volume factor was not
         // included in the charge deposition.
         // Note that the loop is also over ncomps, which takes care of the RZ modes,
         // as well as the old and new rho.
-        amrex::ParallelFor(tb, Rho->nComp(),
+        int const ncomp = Rho->nComp();
+        amrex::ParallelFor(tb, ncomp,
         [=] AMREX_GPU_DEVICE (int i, int j, int /*k*/, int icomp)
         {
             // Wrap the charge density deposited in the guard cells around
             // to the cells above the axis.
             // Rho is located on the boundary
-            if (rmin == 0. && 1-ishift <= i && i <= ngRho-ishift) {
-                Rho_arr(i,j,0,icomp) -= Rho_arr(-ishift-i,j,0,icomp);
+            if (rmin == 0. && 1-ishift <= i && i <= ngRho[0]-ishift) {
+                int imode;
+                if (icomp == 0 || icomp == ncomp/2) {
+                    imode = 0;
+                }
+                else if (icomp < ncomp/2) {
+                    imode = (icomp+1)/2;
+                }
+                else {
+                    imode = (icomp - ncomp/2 + 1)/2;
+                }
+                Rho_arr(i,j,0,icomp) -= std::pow(-1, imode+1)*Rho_arr(-ishift-i,j,0,icomp);
             }
 
             // Apply the inverse volume scaling