New Python test: Particle-Boundary interaction (#4729)

* enable the diagnostic of ParticleScraping in Python

* Update picmi.py

* Update picmi.py

* new test

* python update

* modification of the script

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update PICMI_inputs_rz.py

* json

* update

* Update PICMI_inputs_rz.py

* Update particle_boundary_interaction.json

* Update PICMI_inputs_rz.py

* Update PICMI_inputs_rz.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update PICMI_inputs_rz.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Fix hanging script in parallel

* Make the test executable

* Update analysis script

* Update particle_containers.py

* Update PICMI_inputs_rz.py

* Update analysis.py

* Update analysis.py

* Update particle_containers.py

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py

@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+# @Eya Dammak supervised by @Remi Lehe, 2024
+# --- Input file for particle-boundary interaction testing in RZ.
+# --- This input is a simple case of reflection
+# --- of one electron on the surface of a sphere.
+import numpy as np
+
+from pywarpx import callbacks, particle_containers, picmi
+
+##########################
+# numerics parameters
+##########################
+
+dt = 1.0e-11
+
+# --- Nb time steps
+
+max_steps = 23
+diagnostic_interval = 1
+
+# --- grid
+
+nr = 64
+nz= 64
+
+rmin = 0.0
+rmax = 2
+zmin = -2
+zmax = 2
+
+##########################
+# numerics components
+##########################
+
+grid = picmi.CylindricalGrid(
+    number_of_cells = [nr, nz],
+    n_azimuthal_modes = 1,
+    lower_bound = [rmin, zmin],
+    upper_bound = [rmax, zmax],
+    lower_boundary_conditions = ['none', 'dirichlet'],
+    upper_boundary_conditions =  ['dirichlet', 'dirichlet'],
+    lower_boundary_conditions_particles = ['absorbing', 'reflecting'],
+    upper_boundary_conditions_particles =  ['absorbing', 'reflecting']
+)
+
+
+solver = picmi.ElectrostaticSolver(
+    grid=grid, method='Multigrid',
+    warpx_absolute_tolerance=1e-7
+)
+
+embedded_boundary = picmi.EmbeddedBoundary(
+    implicit_function="-(x**2+y**2+z**2-radius**2)",
+    radius = 0.2
+)
+
+##########################
+# physics components
+##########################
+
+#one particle
+e_dist=picmi.ParticleListDistribution(x=0.0, y=0.0, z=-0.25, ux=0.5e10, uy=0.0, uz=1.0e10, weight=1)
+
+electrons = picmi.Species(
+    particle_type='electron', name='electrons', initial_distribution=e_dist, warpx_save_particles_at_eb=1
+)
+
+##########################
+# diagnostics
+##########################
+
+field_diag = picmi.FieldDiagnostic(
+    name = 'diag1',
+    grid = grid,
+    period = diagnostic_interval,
+    data_list = ['Er', 'Ez', 'phi', 'rho','rho_electrons'],
+    warpx_format = 'openpmd',
+    write_dir = '.',
+    warpx_file_prefix = 'particle_boundary_interaction_plt'
+)
+
+part_diag = picmi.ParticleDiagnostic(name = 'diag1',
+    period = diagnostic_interval,
+    species = [electrons],
+    warpx_format = 'openpmd',
+    write_dir = '.',
+    warpx_file_prefix = 'particle_boundary_interaction_plt'
+)
+
+##########################
+# simulation setup
+##########################
+
+sim = picmi.Simulation(
+    solver=solver,
+    time_step_size = dt,
+    max_steps = max_steps,
+    warpx_embedded_boundary=embedded_boundary,
+    warpx_amrex_the_arena_is_managed=1,
+)
+
+sim.add_species(
+    electrons,
+    layout = picmi.GriddedLayout(
+        n_macroparticle_per_cell=[10, 1, 1], grid=grid
+    )
+)
+sim.add_diagnostic(part_diag)
+sim.add_diagnostic(field_diag)
+
+sim.initialize_inputs()
+sim.initialize_warpx()
+
+##########################
+# python particle data access
+##########################
+
+def concat( list_of_arrays ):
+    if len(list_of_arrays) == 0:
+        # Return a 1d array of size 0
+        return np.empty(0)
+    else:
+        return np.concatenate( list_of_arrays )
+
+
+def mirror_reflection():
+    buffer = particle_containers.ParticleBoundaryBufferWrapper() #boundary buffer
+
+    #STEP 1: extract the different parameters of the boundary buffer (normal, time, position)
+    lev = 0 # level 0 (no mesh refinement here)
+    delta_t = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'deltaTimeScraped', lev))
+    r = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'x', lev))
+    theta = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'theta', lev))
+    z = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'z', lev))
+    x= r*np.cos(theta) #from RZ coordinates to 3D coordinates
+    y= r*np.sin(theta)
+    ux = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'ux', lev))
+    uy = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'uy', lev))
+    uz = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'uz', lev))
+    w = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'w', lev))
+    nx = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'nx', lev))
+    ny = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'ny', lev))
+    nz = concat(buffer.get_particle_boundary_buffer("electrons", 'eb', 'nz', lev))
+
+    #STEP 2: use these parameters to inject particle from the same position in the plasma
+    elect_pc = particle_containers.ParticleContainerWrapper('electrons') #general particle container
+
+    ####this part is specific to the case of simple reflection.
+    un=ux*nx+uy*ny+uz*nz
+    ux_reflect=-2*un*nx+ux #for a "mirror reflection" u(sym)=-2(u.n)n+u
+    uy_reflect=-2*un*ny+uy
+    uz_reflect=-2*un*nz+uz
+    elect_pc.add_particles(
+        x=x + (dt-delta_t)*ux_reflect, y=y + (dt-delta_t)*uy_reflect, z=z + (dt-delta_t)*uz_reflect,
+        ux=ux_reflect, uy=uy_reflect, uz=uz_reflect,
+        w=w
+        ) #adds the particle in the general particle container at the next step
+        #### Can be modified depending on the model of interaction.
+
+    buffer.clear_buffer() #reinitialise the boundary buffer
+
+callbacks.installafterstep(mirror_reflection) #mirror_reflection is called at the next step
+                                              # using the new particle container modified at the last step
+
+##########################
+# simulation run
+##########################
+
+sim.step(max_steps) #the whole process is done "max_steps" times

Examples/Tests/particle_boundary_interaction/analysis.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+# @Eya Dammak supervised by @Remi Lehe, 2024
+"""
+This script tests the last coordinate after adding an electron.
+The sphere is centered on O and has a radius of 0.2 (EB)
+The electron is initially at: (0,0,-0.25) and moves with a velocity:
+(0.5e10,0,1.0e10) with a time step of 1e-11.
+An input file PICMI_inputs_rz.py is used.
+"""
+import os
+import sys
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+import yt
+
+yt.funcs.mylog.setLevel(0)
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# Open plotfile specified in command line
+filename = sys.argv[1]
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename, output_format='openpmd')
+
+ts = OpenPMDTimeSeries('./particle_boundary_interaction_plt')
+
+it=ts.iterations
+x,y,z=ts.get_particle(['x','y','z'], species='electrons', iteration=it[-1])
+
+# Analytical results calculated
+x_analytic=0.03532
+y_analytic=0.00000
+z_analytic=-0.20531
+
+print('NUMERICAL coordinates of the point of contact:')
+print('x=%5.5f, y=%5.5f, z=%5.5f' % (x[0], y[0], z[0]))
+print('\n')
+print('ANALYTICAL coordinates of the point of contact:')
+print('x=%5.5f, y=%5.5f, z=%5.5f' % (x_analytic, y_analytic, z_analytic))
+
+tolerance=1e-5
+
+diff_x=np.abs((x[0]-x_analytic)/x_analytic)
+diff_z=np.abs((z[0]-z_analytic)/z_analytic)
+
+print('\n')
+print("percentage error for x = %5.4f %%" %(diff_x *100))
+print("percentage error for z = %5.4f %%" %(diff_z *100))
+
+assert (diff_x < tolerance) and (y[0] < 1e-8) and (diff_z < tolerance), 'Test particle_boundary_interaction did not pass'

Python/pywarpx/particle_containers.py

@@ -761,6 +761,9 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level)
         The data for the arrays are not copied, but share the underlying
         memory buffer with WarpX. The arrays are fully writeable.
 
+        You can find `here https://github.com/ECP-WarpX/WarpX/blob/319e55b10ad4f7c71b84a4fb21afbafe1f5b65c2/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py`
+        an example of a simple case of particle-boundary interaction (reflection).
+
         Parameters
         ----------
 

Regression/Checksum/benchmarks_json/particle_boundary_interaction.json

@@ -0,0 +1,18 @@
+{
+  "electrons": {
+    "particle_momentum_x": 7.082238783412152e-21,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 7.319015172219235e-21,
+    "particle_position_x": 0.03532003602713794,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.20531131195378569,
+    "particle_weight": 1.0
+  },
+  "lev=0": {
+    "Er": 1.3287679891086577e-06,
+    "Ez": 1.791737740110229e-06,
+    "phi": 3.3983929077087634e-07,
+    "rho": 7.811529217958718e-16,
+    "rho_electrons": 7.811529217958718e-16
+  }
+}

Regression/WarpX-tests.ini

@@ -4681,3 +4681,24 @@ doVis = 0
 compareParticles = 0
 doComparison = 0
 analysisRoutine = Examples/Tests/gaussian_beam/analysis_focusing_beam.py
+
+[particle_boundary_interaction]
+buildDir = .
+inputFile = Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py
+runtime_params =
+customRunCmd = python3 PICMI_inputs_rz.py
+dim = 2
+addToCompileString = USE_PYTHON_MAIN=TRUE USE_RZ=TRUE
+cmakeSetupOpts = -DWarpX_DIMS="RZ" -DWarpX_EB=ON -DWarpX_PYTHON=ON
+target = pip_install
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 1
+particleTypes = electrons
+outputFile = particle_boundary_interaction_plt
+analysisRoutine = Examples/Tests/particle_boundary_interaction/analysis.py