[WIP] Add automated test for particle absorption on new stair-case approximation

Examples/Tests/CMakeLists.txt

@@ -15,6 +15,7 @@ add_subdirectory(electrostatic_sphere)
 add_subdirectory(electrostatic_sphere_eb)
 add_subdirectory(embedded_boundary_cube)
 add_subdirectory(embedded_boundary_diffraction)
+add_subdirectory(embedded_boundary_em_particle_absorption)
 add_subdirectory(embedded_boundary_python_api)
 add_subdirectory(embedded_boundary_rotated_cube)
 add_subdirectory(embedded_circle)

Examples/Tests/embedded_boundary_cube/inputs_base_3d

@@ -13,8 +13,8 @@ boundary.field_lo = pec pec pec
 boundary.field_hi = pec pec pec
 
 eb2.geom_type = box
-eb2.box_lo = -0.5 -0.5 -0.5
-eb2.box_hi = 0.5 0.5 0.5
+eb2.box_lo = -0.501 -0.501 -0.501 # Ensures that the stair-case EB is exactly at -0.5
+eb2.box_hi =  0.501  0.501  0.501 # Ensures that the stair-case EB is exactly at  0.5
 eb2.box_has_fluid_inside = true
 # Alternatively one could use parser to build EB
 # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular.

Examples/Tests/embedded_boundary_cube/inputs_test_2d_embedded_boundary_cube

@@ -12,10 +12,10 @@ warpx.abort_on_warning_threshold = medium
 boundary.field_lo = pec pec
 boundary.field_hi = pec pec
 
-my_constants.xmin = -0.5
-my_constants.zmin = -0.5
-my_constants.xmax = 0.5
-my_constants.zmax = 0.5
+my_constants.xmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5
+my_constants.zmin = -0.501 # Ensures that the stair-case EB is exactly at -0.5
+my_constants.xmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5
+my_constants.zmax = 0.501 # Ensures that the stair-case EB is exactly at 0.5
 # Note that for amrex EB implicit function, >0 is covered, =0 is boundary and <0 is regular.
 warpx.eb_implicit_function = "max(max(x+xmin,-(x+xmax)), max(z+zmin,-(z+zmax)))"
 

Examples/Tests/embedded_boundary_em_particle_absorption/CMakeLists.txt

@@ -0,0 +1,39 @@
+# Add tests (alphabetical order) ##############################################
+#
+
+if(WarpX_EB)
+    add_warpx_test(
+        test_3d_embedded_boundary_em_particle_absorption_sh_factor_1  # name
+        3  # dims
+        1  # nprocs
+        inputs_test_3d_embedded_boundary_em_particle_absorption_sh_factor_1  # inputs
+        "analysis.py"  # analysis
+        "analysis_default_regression.py --path diags/diag1"  # checksum
+        OFF  # dependency
+    )
+endif()
+
+
+if(WarpX_EB)
+    add_warpx_test(
+        test_2d_embedded_boundary_em_particle_absorption_sh_factor_1  # name
+        2  # dims
+        1  # nprocs
+        inputs_test_2d_embedded_boundary_em_particle_absorption_sh_factor_1  # inputs
+        "analysis.py"  # analysis
+        "analysis_default_regression.py --path diags/diag1"  # checksum
+        OFF  # dependency
+    )
+endif()
+
+if(WarpX_EB)
+    add_warpx_test(
+        test_rz_embedded_boundary_em_particle_absorption_sh_factor_1  # name
+        RZ  # dims
+        1  # nprocs
+        inputs_test_rz_embedded_boundary_em_particle_absorption_sh_factor_1  # inputs
+        "analysis.py"  # analysis
+        "analysis_default_regression.py --path diags/diag1"  # checksum
+        OFF  # dependency
+    )
+endif()

Examples/Tests/embedded_boundary_em_particle_absorption/analysis.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+
+"""
+This analysis script checks that there is no spurious charge build-up when a particle is absorbed by an embedded boundary.
+
+More specifically, this test simulates two particles of oppposite charge that are initialized at
+the same position and then move in opposite directions. The particles are surrounded by a cylindrical
+embedded boundary, and are absorbed when their trajectory intersects this boundary. With an
+electromagnetic solver, this can lead to spurious charge build-up (i.e., div(E)!= rho/epsion_0)
+that remains at the position where particle was absorbed.
+
+Note that, in this test, there will also be a (non-spurious) component of div(E) that propagates
+along the embedded boundary, due to electromagnetic waves reflecting on this boundary.
+When checking for static, spurious charge build-up, we average div(E) in time to remove this component.
+
+The test is performed in 2D, 3D and RZ.
+(In 2D, the cylindrical embeded boundary becomes two parallel plates)
+"""
+
+from openpmd_viewer import OpenPMDTimeSeries
+
+ts = OpenPMDTimeSeries("./diags/diag1/")
+
+divE_stacked = ts.iterate(
+    lambda iteration: ts.get_field("divE", iteration=iteration)[0]
+)
+start_avg_iter = 25
+end_avg_iter = 100
+divE_avg = divE_stacked[start_avg_iter:end_avg_iter].mean(axis=0)
+
+# Adjust the tolerance so that the remaining error due to the propagating
+# div(E) (after averaging) is below this tolerance, but so that any typical
+# spurious charge build-up is above this tolerance. This is dimension-dependent.
+dim = ts.fields_metadata["divE"]["geometry"]
+if dim == "3dcartesian":
+    tolerance = 7e-11
+elif dim == "2dcartesian":
+    tolerance = 3.5e-10
+elif dim == "thetaMode":
+    # In RZ: there are issues with divE on axis
+    # Set the few cells around the axis to 0 for this test
+    divE_avg[13:19] = 0
+    tolerance = 4e-12
+
+
+def check_tolerance(array, tolerance):
+    assert abs(array).max() <= tolerance, (
+        f"Test did not pass: the max error {abs(array).max()} exceeded the tolerance of {tolerance}."
+    )
+    print("All elements of are within the tolerance.")
+
+
+check_tolerance(divE_avg, tolerance)

Examples/Tests/embedded_boundary_em_particle_absorption/analysis_default_regression.py

@@ -0,0 +1 @@
+../../analysis_default_regression.py

Examples/Tests/embedded_boundary_em_particle_absorption/inputs_base

@@ -0,0 +1,35 @@
+max_step = 100
+amr.max_level = 0
+amr.blocking_factor = 8
+amr.max_grid_size = 256
+
+algo.charge_deposition = standard
+algo.field_gathering = energy-conserving
+warpx.const_dt = 1.17957283598e-09
+warpx.use_filter = 0
+
+my_constants.R = 6.35
+warpx.eb_implicit_function = "(x**2 + y**2 - R**2)"
+
+particles.species_names = electron positron
+
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "SingleParticle"
+electron.single_particle_pos = 0.0  0.0 0.0
+electron.single_particle_u = 1.e20  0.0 0.4843221e20  # gamma*beta
+electron.single_particle_weight = 1.0
+
+positron.charge = q_e
+positron.mass = m_e
+positron.injection_style = "SingleParticle"
+positron.single_particle_pos = 0.0   0.0 0.0
+positron.single_particle_u = -1.e20  0.0 -0.4843221e20  # gamma*beta
+positron.single_particle_weight = 1.0
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = divE rho
+diag1.format = openpmd

Examples/Tests/embedded_boundary_em_particle_absorption/inputs_test_2d_embedded_boundary_em_particle_absorption_sh_factor_1

@@ -0,0 +1,11 @@
+# base input parameters
+FILE = inputs_base
+
+geometry.dims = 2
+amr.n_cell = 32 32
+geometry.prob_lo     = -10 -10
+geometry.prob_hi     =  10  10
+boundary.field_lo = pec absorbing_silver_mueller
+boundary.field_hi = pec absorbing_silver_mueller
+
+algo.particle_shape = 1

Examples/Tests/embedded_boundary_em_particle_absorption/inputs_test_3d_embedded_boundary_em_particle_absorption_sh_factor_1

@@ -0,0 +1,11 @@
+# base input parameters
+FILE = inputs_base
+
+geometry.dims = 3
+amr.n_cell = 32 32 32
+geometry.prob_lo     = -10 -10 -10
+geometry.prob_hi     =  10  10  10
+boundary.field_lo = pec pec absorbing_silver_mueller
+boundary.field_hi = pec pec absorbing_silver_mueller
+
+algo.particle_shape = 1

Examples/Tests/embedded_boundary_em_particle_absorption/inputs_test_rz_embedded_boundary_em_particle_absorption_sh_factor_1

@@ -0,0 +1,11 @@
+# base input parameters
+FILE = inputs_base
+
+geometry.dims = RZ
+amr.n_cell = 16 32
+geometry.prob_lo     = 0  -10
+geometry.prob_hi     = 10  10
+boundary.field_lo = none absorbing_silver_mueller
+boundary.field_hi = pec absorbing_silver_mueller
+
+algo.particle_shape = 1

Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_cube.json

@@ -3,8 +3,8 @@
     "Bx": 9.263694545408503e-05,
     "By": 0.00031905198933489145,
     "Bz": 7.328424783762594e-05,
-    "Ex": 8553.906698053046,
+    "Ex": 8553.90669811286,
     "Ey": 60867.04830538045,
-    "Ez": 8.439422682267567e-07
+    "Ez": 4.223902107031194e-06
   }
 }

Regression/Checksum/benchmarks_json/test_2d_embedded_boundary_em_particle_absorption_sh_factor_1.json

@@ -0,0 +1,24 @@
+{
+  "lev=0": {
+    "divE": 3.059581906777539e-08,
+    "rho": 0.0
+  },
+  "electron": {
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_weight": 0.0
+  },
+  "positron": {
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_weight": 0.0
+  }
+}

Regression/Checksum/benchmarks_json/test_2d_field_probe.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
     "Bx": 0.0,
-    "By": 126826.78487921853,
+    "By": 123510.69657444415,
     "Bz": 0.0,
-    "Ex": 32517064310550.266,
+    "Ex": 31206368949280.34,
     "Ey": 0.0,
-    "Ez": 17321323003697.61
+    "Ez": 16921005306450.537
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_eb_picmi.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Bx": 148673.005859208,
-    "By": 148673.00585920806,
-    "Bz": 3371.758117878558,
-    "Ex": 55378581103426.71,
-    "Ey": 55378581103426.72,
-    "Ez": 68412803445328.25
+    "Bx": 144495.08082507108,
+    "By": 144495.08082507114,
+    "Bz": 8481.958724628861,
+    "Ex": 54500496182517.92,
+    "Ey": 54500496182517.91,
+    "Ez": 70231240245509.39
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Bx": 4.060477854092961e-18,
+    "Bx": 4.166971025838921e-18,
     "By": 0.006628374119786834,
     "Bz": 0.006628374119786834,
-    "Ex": 5102618.4711524295,
-    "Ey": 6.323754160591239e-05,
-    "Ez": 6.323754160591239e-05
+    "Ex": 5102618.471153786,
+    "Ey": 1.4283859321773714e-05,
+    "Ez": 1.4283859321773714e-05
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_cube_macroscopic.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Bx": 4.20930075273562e-18,
+    "Bx": 4.228863291892693e-18,
     "By": 0.005101824310293573,
     "Bz": 0.005101824310293573,
-    "Ex": 4414725.184731115,
-    "Ey": 6.32375413967707e-05,
-    "Ez": 6.32375413967707e-05
+    "Ex": 4414725.184732471,
+    "Ey": 1.4283895626502055e-05,
+    "Ez": 1.4283895626502055e-05
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_em_particle_absorption_sh_factor_1.json

@@ -0,0 +1,24 @@
+{
+  "lev=0": {
+    "divE": 4.928354322096152e-07,
+    "rho": 0.0
+  },
+  "electron": {
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_weight": 0.0
+  },
+  "positron": {
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.0,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_weight": 0.0
+  }
+}

Regression/Checksum/benchmarks_json/test_3d_embedded_boundary_rotated_cube.json

@@ -1,11 +1,10 @@
 {
   "lev=0": {
-    "Bx": 1.280747509243305e-05,
-    "By": 2.473900144296397e-02,
-    "Bz": 2.473890786894079e-02,
-    "Ex": 1.025322901921306e+07,
-    "Ey": 1.042254197269831e+04,
-    "Ez": 1.040011664019071e+04
+    "Bx": 1.252616939910365e-05,
+    "By": 0.02473895628331097,
+    "Bz": 0.024738956316621142,
+    "Ex": 10253221.850298548,
+    "Ey": 10387.334582977643,
+    "Ez": 10387.532806510022
   }
-}
-
+}

Regression/Checksum/benchmarks_json/test_3d_particle_absorption.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Bx": 202106.71291347666,
-    "By": 202106.71291347663,
-    "Bz": 3371.897999274175,
-    "Ex": 38304043178806.11,
-    "Ey": 38304043178806.11,
-    "Ez": 83057027925874.84
+    "Bx": 198610.0530604908,
+    "By": 198610.0530604909,
+    "Bz": 8482.656173586969,
+    "Ex": 37232105734622.53,
+    "Ey": 37232105734622.54,
+    "Ez": 85094015810307.19
   }
-}
+}

Regression/Checksum/benchmarks_json/test_3d_particle_scrape.json

@@ -1,10 +1,10 @@
 {
   "lev=0": {
-    "Bx": 148673.005859208,
-    "By": 148673.00585920803,
-    "Bz": 3371.758117878557,
-    "Ex": 55378581103426.695,
-    "Ey": 55378581103426.7,
-    "Ez": 68412803445328.25
+    "Bx": 144495.08082507108,
+    "By": 144495.0808250711,
+    "Bz": 8481.95872462886,
+    "Ex": 54500496182517.914,
+    "Ey": 54500496182517.914,
+    "Ez": 70231240245509.4
   }
-}
+}