Correct particle positions outside refined injection (ECP-WarpX#3463)

* Correct particle positions outside refined injection * Improve test * Update benchmark * Update Examples/Physics_applications/laser_acceleration/analysis_refined_injection.py
dpgrote · Nov 22, 2022 · cf011aa · cf011aa
1 parent 0e3760a
commit cf011aa
Show file tree

Hide file tree

Showing 3 changed files with 47 additions and 33 deletions.
diff --git a/Examples/Physics_applications/laser_acceleration/analysis_refined_injection.py b/Examples/Physics_applications/laser_acceleration/analysis_refined_injection.py
@@ -31,8 +31,7 @@
 np = ad['electrons', 'particle_id'].size
 
 # the number of coarse particle streams
-# (odd because one that exactly hits the fine patch is not refined)
-n_coarse = 11
+n_coarse = 10
 
 # the number of fine particle streams
 n_fine = 64
@@ -50,5 +49,15 @@
 
 assert( np == np_expected )
 
+# Test uniformity of rho, by taking a slice of rho that
+# crosses the edge of the refined injection region
+# (but is ahead of the mesh refinement patch)
+ds.force_periodicity()
+ad = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+rho = ad['rho'].to_ndarray().squeeze()
+rho_slice = rho[13:51, 475]
+# Test uniformity up to 0.5% relative variation
+assert( rho_slice.std() < 0.005*abs(rho_slice.mean()) )
+
 test_name = os.path.split(os.getcwd())[1]
 checksumAPI.evaluate_checksum(test_name, fn)
diff --git a/Regression/Checksum/benchmarks_json/RefinedInjection.json b/Regression/Checksum/benchmarks_json/RefinedInjection.json
@@ -8,35 +8,35 @@
     "particle_weight": 12483018148921.525
   },
   "electrons": {
-    "particle_momentum_x": 1.7019663108287842e-19,
-    "particle_momentum_y": 7.804841508613444e-19,
-    "particle_momentum_z": 3.6023951025729006e-19,
-    "particle_position_x": 0.13262784674591108,
-    "particle_position_y": 0.28121954338547195,
-    "particle_weight": 122071289062500.27
+    "particle_momentum_x": 1.701962385014116e-19,
+    "particle_momentum_y": 7.804483053552035e-19,
+    "particle_momentum_z": 3.6023937355618695e-19,
+    "particle_position_x": 0.13035674941360564,
+    "particle_position_y": 0.2792572644186751,
+    "particle_weight": 119232421875000.25
   },
   "lev=0": {
-    "Bx": 24121596.196462017,
-    "By": 78274.29859712059,
-    "Bz": 601058.4634075701,
-    "Ex": 22866959486331.44,
-    "Ey": 6508550078474765.0,
-    "Ez": 32117037637504.49,
-    "jx": 1125098206038016.2,
-    "jy": 2.3251392468384067e+18,
-    "jz": 6796866375752497.0,
-    "rho": 156567078.90542063
+    "Bx": 24121596.46539713,
+    "By": 78274.29859708338,
+    "Bz": 601059.1001251702,
+    "Ex": 22866961334698.695,
+    "Ey": 6508550266827824.0,
+    "Ez": 32117002621363.535,
+    "jx": 1125083620993366.5,
+    "jy": 2.325138134178515e+18,
+    "jz": 6796844511773390.0,
+    "rho": 152886703.9390533
   },
   "lev=1": {
-    "Bx": 53303287.754619256,
-    "By": 213834.35319447558,
-    "Bz": 1321619.3400202503,
-    "Ex": 47849984748795.125,
-    "Ey": 1.4551014572946544e+16,
-    "Ez": 65308104831571.74,
-    "jx": 43901994433620.09,
-    "jy": 9.274835473220729e+18,
-    "jz": 1.9263572931827268e+16,
-    "rho": 67240256.28814039
+    "Bx": 53303287.79524416,
+    "By": 213834.37094616296,
+    "Bz": 1321619.3483672105,
+    "Ex": 47849992572011.23,
+    "Ey": 1.4551014574617832e+16,
+    "Ez": 65308101458272.28,
+    "jx": 43895764938937.25,
+    "jy": 9.274835323597818e+18,
+    "jz": 1.926356834571391e+16,
+    "rho": 67208930.86201674
   }
 }
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
@@ -1130,9 +1130,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
                 ParticleType& p = pp[ip];
                 p.id() = pid+ip;
                 p.cpu() = cpuid;
-
-                const XDim3 r =
-                    inj_pos->getPositionUnitBox(i_part, lrrfac, engine);
+                const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ?
+                  // In the refined injection region: use refinement ratio `lrrfac`
+                  inj_pos->getPositionUnitBox(i_part, lrrfac, engine) :
+                  // Otherwise: use 1 as the refinement ratio
+                  inj_pos->getPositionUnitBox(i_part, 1, engine);
                 auto pos = getCellCoords(overlap_corner, dx, r, iv);
 
 #if defined(WARPX_DIM_3D)
@@ -1645,8 +1647,11 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt)
                 p.cpu() = cpuid;
 
                 // This assumes the inj_pos is of type InjectorPositionRandomPlane
-                const XDim3 r =
-                    inj_pos->getPositionUnitBox(i_part, lrrfac, engine);
+                const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ?
+                  // In the refined injection region: use refinement ratio `lrrfac`
+                  inj_pos->getPositionUnitBox(i_part, lrrfac, engine) :
+                  // Otherwise: use 1 as the refinement ratio
+                  inj_pos->getPositionUnitBox(i_part, 1, engine);
                 auto pos = getCellCoords(overlap_corner, dx, r, iv);
                 auto ppos = PDim3(pos);