Add automated test for particle absorption on new stair-case approximation (#5562)

This introduces automated tests in 2D, 3D and RZ, as a follow up to
#5534, to check that a particle does not leave a spurious charge behind,
when absorbed by the embedded boundary, and when using an EM solver.

In this test, the embedded boundary is a cylinder aligned with the z
axis. (In 2D, this reduces to two parallel plates.)

The tests fail on the `development` branch, before #5534 is merged. They
pass after #5534 is merged.

- `development` branch, before #5534 is merged:


![movie](https://github.com/user-attachments/assets/7d00f181-0744-47bd-ab30-4cc0b9e62956)


- after #5534 is merged:

![movie](https://github.com/user-attachments/assets/dc4af640-502b-4c42-a51d-f30007228061)

Note that these tests fails for higher-order shapes. This will be fixed
in #5209

---------

Co-authored-by: Roelof Groenewald <[email protected]>

Author: RemiLehe

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_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 embedded 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_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_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
+  }
+}