Cropping of particles at boundaries for deposition for charge conservation

Examples/Tests/implicit/analysis_vandb_jfnk_2d.py

@@ -14,7 +14,13 @@
 
 import numpy as np
 import yt
-from scipy.constants import e, epsilon_0
+from scipy.constants import e
+from scipy.constants import epsilon_0 epsilon_0_old
+
+print(f'epsilon_0_old = {epsilon_0_old}')
+
+# Use the same value of epsilon_0 as in Source/ablastr/constant.H
+epsilon_0 = 8.8541878128e-12
 
 field_energy = np.loadtxt("diags/reducedfiles/field_energy.txt", skiprows=1)
 particle_energy = np.loadtxt("diags/reducedfiles/particle_energy.txt", skiprows=1)

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)
@@ -1403,73 +1406,45 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
 #endif
             ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n[ip];
 
-            // computes current and old position in grid units
-#if defined(WARPX_DIM_RZ)
-            amrex::Real const xp_new = xp_np1;
-            amrex::Real const yp_new = yp_np1;
-            amrex::Real const xp_mid = xp_nph;
-            amrex::Real const yp_mid = yp_nph;
-            amrex::Real const xp_old = xp_n[ip];
-            amrex::Real const yp_old = yp_n[ip];
-            amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
-            amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
-            amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
-            amrex::Real costheta_mid, sintheta_mid;
-            if (rp_mid > 0._rt) {
-                costheta_mid = xp_mid/rp_mid;
-                sintheta_mid = yp_mid/rp_mid;
-            } else {
-                costheta_mid = 1._rt;
-                sintheta_mid = 0._rt;
-            }
-            const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
+            // 1) Determine the number of segments.
+            // 2) Loop over segments and deposit current.
+
+            // cell crossings are defined at cell edges if depos_order is odd
+            // cell crossings are defined at cell centers if depos_order is even
+
+            int num_segments = 1;
+            double shift = 0.0;
+            if ( (depos_order % 2) == 0 ) { shift = 0.5; }
+
+#if defined(WARPX_DIM_3D)
+
+            amrex::Real xp_np1_c = xp_np1;
+            amrex::Real yp_np1_c = yp_np1;
+            amrex::Real zp_np1_c = zp_np1;
+            ParticleUtils::crop_at_boundary(xp_n[ip], yp_n[ip], zp_n[ip], xp_np1_c, yp_np1_c, zp_np1_c, xyzmin.x, xyzmax.x, is_absorbing[0]);
+            ParticleUtils::crop_at_boundary(yp_n[ip], zp_n[ip], xp_n[ip], yp_np1_c, zp_np1_c, xp_np1_c, xyzmin.y, xyzmax.y, is_absorbing[1]);
+            ParticleUtils::crop_at_boundary(zp_n[ip], xp_n[ip], yp_n[ip], zp_np1_c, xp_np1_c, yp_np1_c, xyzmin.z, xyzmax.z, is_absorbing[2]);
 
             // Keep these double to avoid bug in single precision
-            double const x_new = (rp_new - 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;
-            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;
-#elif defined(WARPX_DIM_1D_Z)
-            amrex::Real const vx = uxp_nph[ip]*gaminv;
-            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;
+            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;
+            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 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;
-            double const z_old = (zp_n[ip] - xyzmin.z)*dinv.z;
             amrex::Real const vz = (zp_np1 - zp_n[ip])/dt;
 
+            double const dxp = (xp_np1 - xp_n[ip])*dinv.x;
+            double const dyp = (yp_np1 - yp_n[ip])*dinv.y;
+            double const dzp = (zp_np1 - zp_n[ip])*dinv.z;
+
             // Define velocity kernals to deposit
             amrex::Real const wqx = wq*vx*invvol;
             amrex::Real const wqy = wq*vy*invvol;
             amrex::Real const wqz = wq*vz*invvol;
 
-            // 1) Determine the number of segments.
-            // 2) Loop over segments and deposit current.
-
-            // cell crossings are defined at cell edges if depos_order is odd
-            // cell crossings are defined at cell centers if depos_order is even
-
-            int num_segments = 1;
-            double shift = 0.0;
-            if ( (depos_order % 2) == 0 ) { shift = 0.5; }
-
-#if defined(WARPX_DIM_3D)
-
             // compute cell crossings in X-direction
             const auto i_old = static_cast<int>(x_old-shift);
             const auto i_new = static_cast<int>(x_new-shift);
@@ -1494,9 +1469,6 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
 
             // compute total change in particle position and the initial cell
             // locations in each direction used to find the position at cell crossings.
-            const double dxp = x_new - x_old;
-            const double dyp = y_new - y_old;
-            const double dzp = z_new - z_old;
             const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
             const auto dirY_sign = static_cast<double>(dyp < 0. ? -1. : 1.);
             const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
@@ -1661,6 +1633,56 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
 
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
 
+#if defined(WARPX_DIM_RZ)
+            amrex::Real const xp_new = xp_np1;
+            amrex::Real const yp_new = yp_np1;
+            amrex::Real const xp_mid = xp_nph;
+            amrex::Real const yp_mid = yp_nph;
+            amrex::Real const xp_old = xp_n[ip];
+            amrex::Real const yp_old = yp_n[ip];
+            amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new);
+            amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old);
+            amrex::Real const rp_mid = (rp_new + rp_old)/2._rt;
+            amrex::Real costheta_mid, sintheta_mid;
+            if (rp_mid > 0._rt) {
+                costheta_mid = xp_mid/rp_mid;
+                sintheta_mid = yp_mid/rp_mid;
+            } else {
+                costheta_mid = 1._rt;
+                sintheta_mid = 0._rt;
+            }
+            const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid};
+
+            amrex::Real const xp_n_c = rp_old;
+            amrex::Real xp_np1_c = rp_new;
+            amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv;
+#elif defined(WARPX_DIM_XZ)
+            amrex::Real const xp_n_c = xp_n[ip];
+            amrex::Real xp_np1_c = xp_np1;
+            amrex::Real const vy = uyp_nph[ip]*gaminv;
+#endif
+
+            amrex::Real const vx = (xp_np1_c - xp_n_c)/dt;
+            amrex::Real const vz = (zp_np1   - zp_n[ip])/dt;
+
+            amrex::Real zp_np1_c = zp_np1;
+            ParticleUtils::crop_at_boundary(xp_n_c, zp_n[ip], xp_np1_c, zp_np1_c, xyzmin.x, xyzmax.x, is_absorbing[0]);
+            ParticleUtils::crop_at_boundary(zp_n[ip], xp_n_c, zp_np1_c, xp_np1_c, xyzmin.z, xyzmax.z, is_absorbing[1]);
+
+            // Keep these double to avoid bug in single precision
+            double const x_new = (xp_np1_c - xyzmin.x)*dinv.x;
+            double const x_old = (xp_n_c   - xyzmin.x)*dinv.x;
+            double const z_new = (zp_np1_c - xyzmin.z)*dinv.z;
+            double const z_old = (zp_n[ip] - xyzmin.z)*dinv.z;
+
+            double const dxp = (xp_np1 - xp_n[ip])*dinv.x;
+            double const dzp = (zp_np1 - zp_n[ip])*dinv.z;
+
+            // Define velocity kernals to deposit
+            amrex::Real const wqx = wq*vx*invvol;
+            amrex::Real const wqy = wq*vy*invvol;
+            amrex::Real const wqz = wq*vz*invvol;
+
             // compute cell crossings in X-direction
             const auto i_old = static_cast<int>(x_old-shift);
             const auto i_new = static_cast<int>(x_new-shift);
@@ -1679,8 +1701,6 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
 
             // compute total change in particle position and the initial cell
             // locations in each direction used to find the position at cell crossings.
-            const double dxp = x_new - x_old;
-            const double dzp = z_new - z_old;
             const auto dirX_sign = static_cast<double>(dxp < 0. ? -1. : 1.);
             const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
             double Xcell = 0., Zcell = 0.;
@@ -1834,6 +1854,23 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
 
 #elif defined(WARPX_DIM_1D_Z)
 
+            amrex::Real zp_np1_c = zp_np1;
+            ParticleUtils::crop_at_boundary(zp_np1_c, xyzmin.z, xyzmax.z, is_absorbing[0]);
+
+            // Keep these double to avoid bug in single precision
+            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 vx = uxp_nph[ip]*gaminv;
+            amrex::Real const vy = uyp_nph[ip]*gaminv;
+            amrex::Real const vz = (zp_np1 - zp_n[ip])/dt;
+
+            double const dzp = (zp_np1 - zp_n[ip])*dinv.z;
+
+            // Define velocity kernals to deposit
+            amrex::Real const wqx = wq*vx*invvol;
+            amrex::Real const wqy = wq*vy*invvol;
+            amrex::Real const wqz = wq*vz*invvol;
+
             // compute cell crossings in Z-direction
             const auto k_old = static_cast<int>(z_old-shift);
             const auto k_new = static_cast<int>(z_new-shift);
@@ -1845,7 +1882,6 @@ void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::Particle
             // e.g., if (num_segments > 3) ...
 
             // compute dzp and the initial cell location used to find the cell crossings.
-            double const dzp = z_new - z_old;
             const auto dirZ_sign = static_cast<double>(dzp < 0. ? -1. : 1.);
             double Zcell = static_cast<double>(k_old) + shift + 0.5*(1.-dirZ_sign);