Rigid injection: Center field scaling around the v push (#5389)

RemiLehe · web-flow · commit 20a79544daea · 2024-11-01T19:21:25.000Z
In the rigid injection, the fields where scale by the fraction of time
spent between `n*dt` and `(n+1)*dt` to the right of the injection plane.
However, to be consistent with the leap-frog velocity update, this needs
to be between `(n-1/2)*dt` and `(n+1/2)*dt` instead.

As a side-effect of this PR, saving and re-setting `u` and
`optical_depth` to their original value is not needed anymore since the
scaling factor for E and B is 0 for particles to the left of the plane.
diff --git a/Regression/Checksum/benchmarks_json/test_1d_fel.json b/Regression/Checksum/benchmarks_json/test_1d_fel.json
@@ -1,31 +1,31 @@
 {
   "lev=0": {
     "Bx": 0.0,
-    "By": 514.5044890273722,
+    "By": 473.98537926589177,
     "Bz": 0.0,
-    "Ex": 154245109024.33972,
+    "Ex": 142097845843.78326,
     "Ey": 0.0,
     "Ez": 0.0,
-    "jx": 1161126105.5594487,
+    "jx": 1260205974.7220135,
     "jy": 0.0,
     "jz": 0.0
   },
   "electrons": {
     "particle_position_x": 0.0,
     "particle_position_y": 0.0,
-    "particle_position_z": 13607.569953355982,
-    "particle_momentum_x": 3.095483353687591e-19,
+    "particle_position_z": 13607.572916093213,
+    "particle_momentum_x": 3.2646797960476606e-19,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 1.5419514460764825e-16,
+    "particle_momentum_z": 1.541338620507345e-16,
     "particle_weight": 1349823909946836.0
   },
   "positrons": {
     "particle_position_x": 0.0,
     "particle_position_y": 0.0,
-    "particle_position_z": 13607.569953355982,
-    "particle_momentum_x": 3.095483353687591e-19,
+    "particle_position_z": 13607.572916093213,
+    "particle_momentum_x": 3.2646797960476606e-19,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 1.5419514460764825e-16,
+    "particle_momentum_z": 1.541338620507345e-16,
     "particle_weight": 1349823909946836.0
   }
 }
diff --git a/Regression/Checksum/benchmarks_json/test_2d_comoving_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/test_2d_comoving_psatd_hybrid.json
@@ -1,38 +1,38 @@
 {
   "lev=0": {
-    "Bx": 1118808.3686978193,
-    "By": 3248970.5506422943,
-    "Bz": 280612.7921641442,
-    "Ex": 975536649649286.1,
-    "Ey": 402861835403418.1,
-    "Ez": 159049265640492.28,
-    "jx": 2.9996888133195436e+16,
-    "jy": 8.866654944519546e+16,
-    "jz": 3.164008885453435e+17,
-    "rho": 1059988299.6088305
-  },
-  "ions": {
-    "particle_momentum_x": 1.6150513873065298e-18,
-    "particle_momentum_y": 2.233426695677123e-18,
-    "particle_momentum_z": 4.279249529993671e-13,
-    "particle_position_x": 1.4883816864183497,
-    "particle_position_y": 16.452386504127254,
-    "particle_weight": 1.234867369440658e+18
+    "Bx": 1118808.3708538802,
+    "By": 3248949.0437452313,
+    "Bz": 280612.7768961371,
+    "Ex": 975530336896144.1,
+    "Ey": 402861838033488.6,
+    "Ez": 159049784131625.12,
+    "jx": 2.9997142632475216e+16,
+    "jy": 8.866655055001146e+16,
+    "jz": 3.163953981093208e+17,
+    "rho": 1059970922.1974506
   },
   "electrons": {
-    "particle_momentum_x": 7.058167362825288e-19,
-    "particle_momentum_y": 2.204239326446281e-18,
-    "particle_momentum_z": 2.530521998715408e-16,
-    "particle_position_x": 1.5006581263609764,
-    "particle_position_y": 16.454388313398017,
+    "particle_momentum_x": 7.058252826278211e-19,
+    "particle_momentum_y": 2.204239315713169e-18,
+    "particle_momentum_z": 2.530521235191952e-16,
+    "particle_position_x": 1.5006579649318788,
+    "particle_position_y": 16.454388304724286,
     "particle_weight": 1.234867020725368e+18
   },
   "beam": {
-    "particle_momentum_x": 6.869222298759882e-19,
-    "particle_momentum_y": 4.374719809060106e-19,
-    "particle_momentum_z": 6.4523206583503136e-18,
-    "particle_position_x": 0.001290816359726098,
-    "particle_position_y": 0.3586691102823157,
+    "particle_momentum_x": 6.88879318082965e-19,
+    "particle_momentum_y": 4.37466174746362e-19,
+    "particle_momentum_z": 6.4299296650127095e-18,
+    "particle_position_x": 0.0012936414423443238,
+    "particle_position_y": 0.3587414953163842,
     "particle_weight": 3120754537230.3823
+  },
+  "ions": {
+    "particle_momentum_x": 1.6150618501530563e-18,
+    "particle_momentum_y": 2.2334266731098355e-18,
+    "particle_momentum_z": 4.279249530957972e-13,
+    "particle_position_x": 1.488381686539698,
+    "particle_position_y": 16.4523865041322,
+    "particle_weight": 1.234867369440658e+18
   }
 }
diff --git a/Regression/Checksum/benchmarks_json/test_2d_galilean_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/test_2d_galilean_psatd_hybrid.json
@@ -1,38 +1,38 @@
 {
   "lev=0": {
-    "Bx": 1086729.9718613266,
-    "By": 2886554.482275311,
-    "Bz": 264259.55093734514,
-    "Ex": 867387781289915.2,
-    "Ey": 392666724461952.5,
-    "Ez": 146897592531660.03,
-    "jx": 2.702866174672266e+16,
-    "jy": 8.615938361747776e+16,
-    "jz": 2.7329155817806224e+17,
-    "rho": 915945723.7934376
+    "Bx": 1086729.9879225595,
+    "By": 2886531.8361456757,
+    "Bz": 264259.55266959703,
+    "Ex": 867381192933999.0,
+    "Ey": 392666738858258.7,
+    "Ez": 146898030091111.84,
+    "jx": 2.702892158065604e+16,
+    "jy": 8.615938867870698e+16,
+    "jz": 2.7328506574305683e+17,
+    "rho": 915924567.6956444
+  },
+  "beam": {
+    "particle_momentum_x": 7.006049777955171e-19,
+    "particle_momentum_y": 4.374916846096741e-19,
+    "particle_momentum_z": 6.173292885825711e-18,
+    "particle_position_x": 0.0016046573777589298,
+    "particle_position_y": 0.35899824939059793,
+    "particle_weight": 3120754537230.3823
   },
   "ions": {
-    "particle_momentum_x": 1.4394902513923003e-18,
-    "particle_momentum_y": 1.5967629157922875e-18,
-    "particle_momentum_z": 4.287340658051679e-13,
-    "particle_position_x": 1.4911814217142487,
-    "particle_position_y": 16.521964978771,
+    "particle_momentum_x": 1.4395010524514718e-18,
+    "particle_momentum_y": 1.596762923413923e-18,
+    "particle_momentum_z": 4.2873406589510426e-13,
+    "particle_position_x": 1.4911814218338595,
+    "particle_position_y": 16.52196497877563,
     "particle_weight": 1.2372405194129536e+18
   },
   "electrons": {
-    "particle_momentum_x": 6.240933687389075e-19,
-    "particle_momentum_y": 1.5790611427694247e-18,
-    "particle_momentum_z": 2.5064357834741096e-16,
-    "particle_position_x": 1.501413766926399,
-    "particle_position_y": 16.523781713952324,
+    "particle_momentum_x": 6.241019323257125e-19,
+    "particle_momentum_y": 1.5790611706036782e-18,
+    "particle_momentum_z": 2.5064350576384073e-16,
+    "particle_position_x": 1.501413593465263,
+    "particle_position_y": 16.52378170448397,
     "particle_weight": 1.2372401466086835e+18
-  },
-  "beam": {
-    "particle_momentum_x": 7.000932845220306e-19,
-    "particle_momentum_y": 4.374936866729326e-19,
-    "particle_momentum_z": 6.194468548032543e-18,
-    "particle_position_x": 0.0016030835496557787,
-    "particle_position_y": 0.3589262705964349,
-    "particle_weight": 3120754537230.3823
   }
 }
diff --git a/Regression/Checksum/benchmarks_json/test_2d_laser_acceleration_boosted.json b/Regression/Checksum/benchmarks_json/test_2d_laser_acceleration_boosted.json
@@ -1,38 +1,38 @@
 {
   "lev=0": {
-    "Bx": 4818955.480797943,
-    "By": 1752.8025638791275,
-    "Bz": 14516.212782554387,
-    "Ex": 2366115503598.9224,
-    "Ey": 1446112025635674.2,
-    "Ez": 21864189507357.867,
-    "jx": 1996366349775593.5,
-    "jy": 5.312583827155926e+16,
-    "jz": 2.0491352624508764e+16,
-    "rho": 68443961.71852128
+    "Bx": 4818955.485307051,
+    "By": 1752.8020185365554,
+    "Bz": 14516.212849649737,
+    "Ex": 2366115529014.2324,
+    "Ey": 1446112026998942.5,
+    "Ez": 21864189485739.55,
+    "jx": 1996366372981548.5,
+    "jy": 5.312583836344946e+16,
+    "jz": 2.049135259966133e+16,
+    "rho": 68443961.64027263
+  },
+  "beam": {
+    "particle_momentum_x": 3.535736052190267e-19,
+    "particle_momentum_y": 4.363217976210739e-19,
+    "particle_momentum_z": 5.658515465395611e-17,
+    "particle_position_x": 0.008314855161869274,
+    "particle_position_y": 1.170433573157185,
+    "particle_weight": 62415090744.60765
   },
   "electrons": {
-    "particle_momentum_x": 2.2135945391319113e-23,
-    "particle_momentum_y": 2.8224559499558413e-22,
-    "particle_momentum_z": 5.260626010214114e-22,
-    "particle_position_x": 0.010800577787628052,
-    "particle_position_y": 0.2111506062831815,
+    "particle_momentum_x": 2.213594541883545e-23,
+    "particle_momentum_y": 2.8224559261549207e-22,
+    "particle_momentum_z": 5.260626007410037e-22,
+    "particle_position_x": 0.010800577787636243,
+    "particle_position_y": 0.2111506062831794,
     "particle_weight": 4.121554826246186e+16
   },
   "ions": {
-    "particle_momentum_x": 6.248472277246885e-23,
-    "particle_momentum_y": 4.449097689427654e-22,
-    "particle_momentum_z": 5.768168724998047e-22,
+    "particle_momentum_x": 6.24847229412907e-23,
+    "particle_momentum_y": 4.449097671673176e-22,
+    "particle_momentum_z": 5.768168722032957e-22,
     "particle_position_x": 0.010800001678510512,
     "particle_position_y": 0.21114947608115425,
     "particle_weight": 4.121554826246186e+16
-  },
-  "beam": {
-    "particle_momentum_x": 3.5357456351701565e-19,
-    "particle_momentum_y": 4.363391839372122e-19,
-    "particle_momentum_z": 5.658606416951653e-17,
-    "particle_position_x": 0.008314723025211468,
-    "particle_position_y": 1.1704335743854242,
-    "particle_weight": 62415090744.60765
   }
 }
diff --git a/Source/Particles/Gather/ScaleFields.H b/Source/Particles/Gather/ScaleFields.H
@@ -47,8 +47,21 @@ struct ScaleFields
         // This only approximates what should be happening. The particles
         // should by advanced a fraction of a time step instead.
         // Scaling the fields is much easier and may be good enough.
-        const amrex::Real dtscale = 1._rt - (m_z_plane_previous - zp)/(m_vz_ave_boosted + m_v_boost)/m_dt;
-        if (0._rt < dtscale && dtscale < 1._rt)
+
+        // The scaling factor corresponds to the fraction of time that
+        // the particles spends to the right of the injection plane,
+        // between (n-1/2)*dt and (n+1/2)*dt, which is the interval
+        // over which the velocity is updated, in the leap-frog velocity push.
+        // (Note that here, `zp` is the particle position at time n*dt)
+        amrex::Real dtscale = 0.5_rt - (m_z_plane_previous - zp)/(m_vz_ave_boosted + m_v_boost)/m_dt;
+        // If the particle stays to the left of the plane during the
+        // whole push, simply set the scaling factor to 0, and thus
+        // the velocity push leaves the velocity unchanged.
+        if (dtscale < 0._rt) {
+            dtscale = 0;
+        }
+        // Scale the fields.
+        if (dtscale < 1._rt)
         {
             Exp *= dtscale;
             Eyp *= dtscale;
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -176,10 +176,6 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
 
     // Save the position, momentum and optical depth, making copies
     amrex::Gpu::DeviceVector<ParticleReal> xp_save, yp_save, zp_save;
-    amrex::Gpu::DeviceVector<ParticleReal> uxp_save, uyp_save, uzp_save;
-#ifdef WARPX_QED
-    amrex::Gpu::DeviceVector<ParticleReal> optical_depth_save;
-#endif
 
     const auto GetPosition = GetParticlePosition<PIdx>(pti, offset);
           auto SetPosition = SetParticlePosition<PIdx>(pti, offset);
@@ -188,54 +184,24 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
     amrex::ParticleReal* const AMREX_RESTRICT uy = uyp.dataPtr() + offset;
     amrex::ParticleReal* const AMREX_RESTRICT uz = uzp.dataPtr() + offset;
 
-#ifdef WARPX_QED
-    const bool loc_has_quantum_sync = has_quantum_sync();
-    amrex::ParticleReal* AMREX_RESTRICT p_optical_depth = nullptr;
-    amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_save = nullptr;
-#endif
-
     if (!done_injecting_lev)
     {
         // If the old values are not already saved, create copies here.
         xp_save.resize(np_to_push);
         yp_save.resize(np_to_push);
         zp_save.resize(np_to_push);
 
-        uxp_save.resize(np_to_push);
-        uyp_save.resize(np_to_push);
-        uzp_save.resize(np_to_push);
-
         amrex::ParticleReal* const AMREX_RESTRICT xp_save_ptr = xp_save.dataPtr();
         amrex::ParticleReal* const AMREX_RESTRICT yp_save_ptr = yp_save.dataPtr();
         amrex::ParticleReal* const AMREX_RESTRICT zp_save_ptr = zp_save.dataPtr();
 
-        amrex::ParticleReal* const AMREX_RESTRICT uxp_save_ptr = uxp_save.dataPtr();
-        amrex::ParticleReal* const AMREX_RESTRICT uyp_save_ptr = uyp_save.dataPtr();
-        amrex::ParticleReal* const AMREX_RESTRICT uzp_save_ptr = uzp_save.dataPtr();
-
-#ifdef WARPX_QED
-        if(loc_has_quantum_sync){
-            p_optical_depth = pti.GetAttribs(particle_comps["opticalDepthQSR"]).dataPtr()
-                              + offset;
-            optical_depth_save.resize(np_to_push);
-            p_optical_depth_save = optical_depth_save.dataPtr();
-        }
-#endif
-
         amrex::ParallelFor( np_to_push,
                             [=] AMREX_GPU_DEVICE (long i) {
                                 amrex::ParticleReal xp, yp, zp;
                                 GetPosition(i, xp, yp, zp);
                                 xp_save_ptr[i] = xp;
                                 yp_save_ptr[i] = yp;
                                 zp_save_ptr[i] = zp;
-                                uxp_save_ptr[i] = ux[i];
-                                uyp_save_ptr[i] = uy[i];
-                                uzp_save_ptr[i] = uz[i];
-#ifdef WARPX_QED
-                                if(loc_has_quantum_sync){
-                                    p_optical_depth_save[i] = p_optical_depth[i];}
-#endif
                             });
     }
 
@@ -252,9 +218,6 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
         amrex::ParticleReal* AMREX_RESTRICT x_save = xp_save.dataPtr();
         amrex::ParticleReal* AMREX_RESTRICT y_save = yp_save.dataPtr();
         amrex::ParticleReal* AMREX_RESTRICT z_save = zp_save.dataPtr();
-        amrex::ParticleReal* AMREX_RESTRICT ux_save = uxp_save.dataPtr();
-        amrex::ParticleReal* AMREX_RESTRICT uy_save = uyp_save.dataPtr();
-        amrex::ParticleReal* AMREX_RESTRICT uz_save = uzp_save.dataPtr();
 
         // Undo the push for particles not injected yet.
         // The zp are advanced a fixed amount.
@@ -267,9 +230,6 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
                                 amrex::ParticleReal xp, yp, zp;
                                 GetPosition(i, xp, yp, zp);
                                 if (zp <= z_plane_lev) {
-                                    ux[i] = ux_save[i];
-                                    uy[i] = uy_save[i];
-                                    uz[i] = uz_save[i];
                                     xp = x_save[i];
                                     yp = y_save[i];
                                     if (rigid) {
@@ -281,10 +241,6 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
                                         zp = z_save[i] + dt*uz[i]*gi;
                                     }
                                     SetPosition(i, xp, yp, zp);
-#ifdef WARPX_QED
-                                    if(loc_has_quantum_sync){
-                                        p_optical_depth[i] = p_optical_depth_save[i];}
-#endif
                                 }
                             });
     }
