Cropping of particles at boundaries for deposition for charge conservation

Source/Particles/Deposition/CurrentDeposition.H

@@ -12,6 +12,7 @@
 #include "ablastr/parallelization/KernelTimer.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/ShapeFactors.H"
+#include "Utils/ParticleUtils.H"
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXAlgorithmSelection.H"
 #include "Utils/WarpXConst.H"
@@ -1352,6 +1353,8 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
                                            const amrex::Real dt,
                                            const amrex::XDim3 & dinv,
                                            const amrex::XDim3 & xyzmin,
+                                           const amrex::XDim3 & xyzmax,
+                                           const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & is_absorbing,
                                            const amrex::Dim3 lo,
                                            const amrex::Real q,
                                            [[maybe_unused]] const int n_rz_azimuthal_modes)
@@ -1425,13 +1428,15 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
             const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
 
             // Keep these double to avoid bug in single precision
-            double const x_new = (rp_new - xyzmin.x)*dinv.x;
+            amrex::Real const rp_new_c = ParticleUtils::crop_at_boundary(rp_new, xyzmin.x, xyzmax.x, is_absorbing[0]);
+            double const x_new = (rp_new_c - xyzmin.x)*dinv.x;
             double const x_old = (rp_old - xyzmin.x)*dinv.x;
             amrex::Real const vx = (rp_new - rp_old)/dt;
             amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
 #elif defined(WARPX_DIM_XZ)
             // Keep these double to avoid bug in single precision
-            double const x_new = (xp_np1 - xyzmin.x)*dinv.x;
+            amrex::Real const xp_np1_c = ParticleUtils::crop_at_boundary(xp_np1, xyzmin.x, xyzmax.x, is_absorbing[0]);
+            double const x_new = (xp_np1_c - xyzmin.x)*dinv.x;
             double const x_old = (xp_n[ip] - xyzmin.x)*dinv.x;
             amrex::Real const vx = (xp_np1 - xp_n[ip])/dt;
             amrex::Real const vy = uyp_nph[ip]*gaminv;
@@ -1440,16 +1445,19 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
             amrex::Real const vy = uyp_nph[ip]*gaminv;
 #elif defined(WARPX_DIM_3D)
             // Keep these double to avoid bug in single precision
-            double const x_new = (xp_np1 - xyzmin.x)*dinv.x;
+            amrex::Real const xp_np1_c = ParticleUtils::crop_at_boundary(xp_np1, xyzmin.x, xyzmax.x, is_absorbing[0]);
+            amrex::Real const yp_np1_c = ParticleUtils::crop_at_boundary(yp_np1, xyzmin.y, xyzmax.y, is_absorbing[1]);
+            double const x_new = (xp_np1_c - xyzmin.x)*dinv.x;
             double const x_old = (xp_n[ip] - xyzmin.x)*dinv.x;
-            double const y_new = (yp_np1 - xyzmin.y)*dinv.y;
+            double const y_new = (yp_np1_c - xyzmin.y)*dinv.y;
             double const y_old = (yp_n[ip] - xyzmin.y)*dinv.y;
             amrex::Real const vx = (xp_np1 - xp_n[ip])/dt;
             amrex::Real const vy = (yp_np1 - yp_n[ip])/dt;
 #endif
 
             // Keep these double to avoid bug in single precision
-            double const z_new = (zp_np1 - xyzmin.z)*dinv.z;
+            amrex::Real const zp_np1_c = ParticleUtils::crop_at_boundary(zp_np1, xyzmin.z, xyzmax.z, is_absorbing[WARPX_ZINDEX]);
+            double const z_new = (zp_np1_c - xyzmin.z)*dinv.z;
             double const z_old = (zp_n[ip] - xyzmin.z)*dinv.z;
             amrex::Real const vz = (zp_np1 - zp_n[ip])/dt;
 

Source/Particles/Gather/FieldGather.H

@@ -12,6 +12,7 @@
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/ShapeFactors.H"
 #include "Utils/WarpX_Complex.H"
+#include "Utils/ParticleUtils.H"
 
 #include <AMReX.H>
 
@@ -887,6 +888,8 @@ void doGatherPicnicShapeN (
                      [[maybe_unused]] const amrex::IndexType Bz_type,
                      const amrex::XDim3 & dinv,
                      const amrex::XDim3 & xyzmin,
+                     const amrex::XDim3 & xyzmax,
+                     const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & is_absorbing,
                      const amrex::Dim3& lo,
                      const int n_rz_azimuthal_modes)
 {
@@ -928,26 +931,30 @@ void doGatherPicnicShapeN (
         sintheta_mid = 0._rt;
     }
     const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
+    amrex::Real const rp_new_c = ParticleUtils::crop_at_boundary(rp_new, xyzmin.x, xyzmax.x, is_absorbing[0]);
     // Keep these double to avoid bug in single precision
-    double const x_new = (rp_new - xyzmin.x)*dinv.x;
+    double const x_new = (rp_new_c - xyzmin.x)*dinv.x;
     double const x_old = (rp_old - xyzmin.x)*dinv.x;
-    double const x_bar = (rp_mid - xyzmin.x)*dinv.x;
+    double const x_bar = (x_new + x_old)*0.5_rt;
 #elif !defined(WARPX_DIM_1D_Z)
+    amrex::Real const xp_np1_c = ParticleUtils::crop_at_boundary(xp_np1, xyzmin.x, xyzmax.x, is_absorbing[0]);
     // Keep these double to avoid bug in single precision
-    double const x_new = (xp_np1 - xyzmin.x)*dinv.x;
+    double const x_new = (xp_np1_c - xyzmin.x)*dinv.x;
     double const x_old = (xp_n - xyzmin.x)*dinv.x;
-    double const x_bar = (xp_nph - xyzmin.x)*dinv.x;
+    double const x_bar = (x_new + x_old)*0.5_rt;
 #endif
 #if defined(WARPX_DIM_3D)
+    amrex::Real const yp_np1_c = ParticleUtils::crop_at_boundary(yp_np1, xyzmin.y, xyzmax.y, is_absorbing[1]);
     // Keep these double to avoid bug in single precision
-    double const y_new = (yp_np1 - xyzmin.y)*dinv.y;
+    double const y_new = (yp_np1_c - xyzmin.y)*dinv.y;
     double const y_old = (yp_n   - xyzmin.y)*dinv.y;
-    double const y_bar = (yp_nph - xyzmin.y)*dinv.y;
+    double const y_bar = (y_new + y_old)*0.5_rt;
 #endif
+    amrex::Real const zp_np1_c = ParticleUtils::crop_at_boundary(zp_np1, xyzmin.z, xyzmax.z, is_absorbing[WARPX_ZINDEX]);
     // Keep these double to avoid bug in single precision
-    double const z_new = (zp_np1 - xyzmin.z)*dinv.z;
+    double const z_new = (zp_np1_c - xyzmin.z)*dinv.z;
     double const z_old = (zp_n   - xyzmin.z)*dinv.z;
-    double const z_bar = (zp_nph - xyzmin.z)*dinv.z;
+    double const z_bar = (z_new + z_old)*0.5_rt;
 
     // 1) Determine the number of segments.
     // 2) Loop over segments and gather electric field.
@@ -1711,6 +1718,8 @@ void doGatherShapeNImplicit (
                      const amrex::IndexType bz_type,
                      const amrex::XDim3 & dinv,
                      const amrex::XDim3 & xyzmin,
+                     const amrex::XDim3 & xyzmax,
+                     const amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> & is_absorbing,
                      const amrex::Dim3& lo,
                      const int n_rz_azimuthal_modes,
                      const int nox,
@@ -1749,25 +1758,25 @@ void doGatherShapeNImplicit (
                                     Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                     ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
                                     ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
-                                    dinv, xyzmin, lo, n_rz_azimuthal_modes);
+                                    dinv, xyzmin, xyzmax, is_absorbing, lo, n_rz_azimuthal_modes);
         } else if (nox == 2) {
             doGatherPicnicShapeN<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
                                     Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                     ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
                                     ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
-                                    dinv, xyzmin, lo, n_rz_azimuthal_modes);
+                                    dinv, xyzmin, xyzmax, is_absorbing, lo, n_rz_azimuthal_modes);
         } else if (nox == 3) {
             doGatherPicnicShapeN<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
                                     Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                     ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
                                     ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
-                                    dinv, xyzmin, lo, n_rz_azimuthal_modes);
+                                    dinv, xyzmin, xyzmax, is_absorbing, lo, n_rz_azimuthal_modes);
         } else if (nox == 4) {
             doGatherPicnicShapeN<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph,
                                     Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                     ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
                                     ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
-                                    dinv, xyzmin, lo, n_rz_azimuthal_modes);
+                                    dinv, xyzmin, xyzmax, is_absorbing, lo, n_rz_azimuthal_modes);
         }
     }
     else if (depos_type == CurrentDepositionAlgo::Direct) {

Source/Particles/PhysicalParticleContainer.cpp

@@ -2832,6 +2832,22 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti,
 
     // Lower corner of tile box physical domain (take into account Galilean shift)
     const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt);
+    const amrex::XDim3 xyzmax = WarpX::UpperCorner(box, gather_lev, 0._rt);
+
+    const WarpX& warpx = WarpX::GetInstance();
+    amrex::Box const& domain_box = warpx.Geom(0).Domain();
+    auto & field_boundary_lo = warpx.GetFieldBoundaryLo();
+    auto & field_boundary_hi = warpx.GetFieldBoundaryHi();
+
+    amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> is_absorbing;
+    for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
+        is_absorbing[idim][0] = (box.smallEnd(idim) <= domain_box.smallEnd(idim) &&
+                                 (field_boundary_lo[idim] == FieldBoundaryType::PEC
+                               || field_boundary_lo[idim] == FieldBoundaryType::PMC));
+        is_absorbing[idim][1] = (box.bigEnd(idim) >= domain_box.bigEnd(idim) &&
+                                 (field_boundary_hi[idim] == FieldBoundaryType::PEC
+                               || field_boundary_hi[idim] == FieldBoundaryType::PMC));
+    }
 
     const Dim3 lo = lbound(box);
 
@@ -2986,7 +3002,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti,
                 doGatherShapeNImplicit(xp_n, yp_n, zp_n, xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                        ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
                                        ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
-                                       dinv, xyzmin, lo, n_rz_azimuthal_modes, nox,
+                                       dinv, xyzmin, xyzmax, is_absorbing, lo, n_rz_azimuthal_modes, nox,
                                        depos_type );
             }
 

Source/Particles/WarpXParticleContainer.cpp

@@ -477,6 +477,21 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti,
     const Dim3 lo = lbound(tilebox);
     // Take into account Galilean shift
     const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, depos_lev, 0.5_rt*dt);
+    const amrex::XDim3 xyzmax = WarpX::UpperCorner(tilebox, depos_lev, 0.5_rt*dt);
+    amrex::Box const& domain_box = warpx.Geom(0).Domain();
+
+    auto & field_boundary_lo = warpx.GetFieldBoundaryLo();
+    auto & field_boundary_hi = warpx.GetFieldBoundaryHi();
+
+    amrex::GpuArray<amrex::GpuArray<bool,2>, AMREX_SPACEDIM> is_absorbing;
+    for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
+        is_absorbing[idim][0] = (tilebox.smallEnd(idim) <= domain_box.smallEnd(idim) &&
+                                 (field_boundary_lo[idim] == FieldBoundaryType::PEC
+                               || field_boundary_lo[idim] == FieldBoundaryType::PMC));
+        is_absorbing[idim][1] = (tilebox.bigEnd(idim) >= domain_box.bigEnd(idim) &&
+                                 (field_boundary_hi[idim] == FieldBoundaryType::PEC
+                               || field_boundary_hi[idim] == FieldBoundaryType::PMC));
+    }
 
     if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov ||
         WarpX::current_deposition_algo == CurrentDepositionAlgo::Villasenor) {
@@ -702,31 +717,31 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti,
                         GetPosition, wp.dataPtr() + offset,
                         uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset,
                         uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
-                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q,
+                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, xyzmax, is_absorbing, lo, q,
                         WarpX::n_rz_azimuthal_modes);
                 } else if (WarpX::nox == 2){
                     doVillasenorDepositionShapeNImplicit<2>(
                         xp_n_data, yp_n_data, zp_n_data,
                         GetPosition, wp.dataPtr() + offset,
                         uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset,
                         uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
-                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q,
+                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, xyzmax, is_absorbing, lo, q,
                         WarpX::n_rz_azimuthal_modes);
                 } else if (WarpX::nox == 3){
                     doVillasenorDepositionShapeNImplicit<3>(
                         xp_n_data, yp_n_data, zp_n_data,
                         GetPosition, wp.dataPtr() + offset,
                         uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset,
                         uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
-                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q,
+                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, xyzmax, is_absorbing, lo, q,
                         WarpX::n_rz_azimuthal_modes);
                 } else if (WarpX::nox == 4){
                     doVillasenorDepositionShapeNImplicit<4>(
                         xp_n_data, yp_n_data, zp_n_data,
                         GetPosition, wp.dataPtr() + offset,
                         uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset,
                         uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev,
-                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q,
+                        jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, xyzmax, is_absorbing, lo, q,
                         WarpX::n_rz_azimuthal_modes);
                 }
             }

Source/Utils/ParticleUtils.H

@@ -189,6 +189,19 @@ namespace ParticleUtils {
                           && (xlo[2] <= point.z) && (point.z <= xhi[2]));
     }
 
+    /* \brief Crops the position at the specified boundary if check is true
+     * \param[in] x The particle position
+     * \param[in] xmin The lower boundary of the domain
+     * \param[in] xmax The upper boundary of the domain
+     * \param[in] check Whether to crop at each of the boundaries
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    amrex::Real crop_at_boundary(amrex::Real x, amrex::Real xmin, amrex::Real xmax, amrex::GpuArray<bool,2> const& check) {
+        if (check[0] && x < xmin) { return xmin; }
+        if (check[1] && x > xmax) { return xmax; }
+        return x;
+    }
+
 }
 
 #endif // WARPX_PARTICLE_UTILS_H_