[WIP] Save the normal components in the particle buffer when there is EB

Docs/source/usage/parameters.rst

@@ -2752,8 +2752,9 @@ The data collected at each boundary is written out to a subdirectory of the diag
 By default, all of the collected particle data is written out at the end of the simulation. Optionally, the ``<diag_name>.intervals`` parameter can be given to specify writing out the data more often.
 This can be important if a large number of particles are lost, avoiding filling up memory with the accumulated lost particle data.
 
-In addition to their usual attributes, the saved particles have an integer attribute ``timestamp``, which
-indicates the PIC iteration at which each particle was absorbed at the boundary.
+In addition to their usual attributes, the saved particles have real attributes
+``timestamp``, which indicates the exact time at which each particle was absorbed at the boundary ;
+And ``nx``, ``ny``, ``nz``, which indicate the components of the normal vector to the EB on the exact position.
 
 ``BoundaryScrapingDiagnostics`` can be used with ``<diag_name>.<species>.random_fraction``, ``<diag_name>.<species>.uniform_stride``, and ``<diag_name>.<species>.plot_filter_function``, which have the same behavior as for ``FullDiagnostics``. For ``BoundaryScrapingDiagnostics``, these filters are applied at the time the data is written to file. An implication of this is that more particles may initially be accumulated in memory than are ultimately written. ``t`` in ``plot_filter_function`` refers to the time the diagnostic is written rather than the time the particle crossed the boundary.
 

Examples/Tests/particle_boundary_process/PICMI_inputs_reflection.py

@@ -131,12 +131,12 @@
 print("Number of electrons in lower buffer:", n)
 assert n == 67
 
-scraped_steps = buffer.get_particle_boundary_buffer("electrons", 'z_hi', 'step_scraped', 0)
+scraped_steps = buffer.get_particle_boundary_buffer("electrons", 'z_hi', 'timestamp', 0)
 for arr in scraped_steps:
     # print(arr)
     assert all(arr == 4)
 
-scraped_steps = buffer.get_particle_boundary_buffer("electrons", 'z_lo', 'step_scraped', 0)
+scraped_steps = buffer.get_particle_boundary_buffer("electrons", 'z_lo', 'timestamp', 0)
 for arr in scraped_steps:
     # print(arr)
     assert all(arr == 8)

Examples/Tests/particle_boundary_scrape/PICMI_inputs_scrape.py

@@ -133,7 +133,7 @@
 print(f"Number of electrons in buffer (proc #{my_id}): {n}")
 assert n == 612
 
-scraped_steps = particle_buffer.get_particle_boundary_buffer("electrons", 'eb', 'step_scraped', 0)
+scraped_steps = particle_buffer.get_particle_boundary_buffer("electrons", 'eb', 'timestamp', 0)
 for arr in scraped_steps:
     assert all(np.array(arr, copy=False) > 40)
 

Examples/Tests/point_of_contact_EB/analysis.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+
+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_scraping = OpenPMDTimeSeries('./diags/diag2/particles_at_eb/')
+
+
+it=ts_scraping.iterations
+timestamp,x,y,z,nx,ny,nz=ts_scraping.get_particle( ['timestamp','x','y','z','nx','ny','nz'], species='electron', iteration=it )
+
+timestamp_analytic=3.58
+x_analytic=-0.1983
+y_analytic=0.02584
+z_analytic=0.0000
+nx_analytic=-0.99
+ny_analytic=0.13
+nz_analytic=0.0
+
+timestamp=timestamp[0]*1e10
+print('\n')
+print('NUMERICAL coordinates of the point of contact:')
+print('time_stamp=%5.4f e-10, x=%5.4f, y=%5.4f, z=%5.4f, nx=%5.4f, ny=%5.4f, nz=%5.4f' % (timestamp, x[0], y[0], z[0], nx[0], ny[0], nz[0]))
+print('\n')
+print('ANALYTICAL coordinates of the point of contact:')
+print('time_stamp=%5.2f, x=%5.4f, y=%5.4f, z=%5.4f, nx=%5.4f, ny=%5.4f, nz=%5.4f' % (timestamp_analytic, x_analytic, y_analytic, z_analytic, nx_analytic, ny_analytic, nz_analytic))
+
+tolerance=0.001
+tolerance_t=0.003
+tolerance_n=0.01
+print('\n')
+print("tolerance = "+ str(tolerance *100) + '%')
+print("tolerance for the time_stamp = "+ str(tolerance_t *100) + '%')
+print("tolerance for the normal components = "+ str(tolerance_n *100) + '%')
+diff_timestamp=np.abs((timestamp-timestamp_analytic)/timestamp_analytic)
+diff_x=np.abs((x[0]-x_analytic)/x_analytic)
+diff_y=np.abs((y[0]-y_analytic)/y_analytic)
+diff_nx=np.abs((nx[0]-nx_analytic)/nx_analytic)
+diff_ny=np.abs((ny[0]-ny_analytic)/ny_analytic)
+
+print('\n')
+print("percentage error for timestamp = %5.4f %%" %(diff_timestamp *100))
+print("percentage error for x = %5.4f %%" %(diff_x *100))
+print("percentage error for y = %5.4f %%" %(diff_y *100))
+print("percentage error for nx = %5.2f %%" %(diff_nx *100))
+print("percentage error for ny = %5.2f %%" %(diff_ny *100))
+print("nz = %5.2f " %(nz[0]))
+
+
+assert (diff_timestamp < tolerance_t) and (diff_x < tolerance) and (diff_y < tolerance) and (diff_nx < tolerance_n) and (diff_ny < tolerance_n) and (np.abs(nz) < tolerance_n), 'Test point_of_contact did not pass'

Examples/Tests/point_of_contact_EB/analysis_rz.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+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_scraping = OpenPMDTimeSeries('./diags/diag2/particles_at_eb/')
+
+
+it=ts_scraping.iterations
+r,theta,z=ts_scraping.get_particle( ['x','theta','z'], species='electron', iteration=it )
+
+r_analytic=0.2
+theta_analytic=3.0
+z_analytic=0.0000
+
+print('NUMERICAL coordinates of the point of contact:')
+print('r=%5.4f, theta=%5.4f, z=%5.4f' % (r[0], theta[0], z[0]))
+print('\n')
+print('ANALYTICAL coordinates of the point of contact:')
+print('x=%5.4f, y=%5.4f, z=%5.4f' % (r_analytic, theta_analytic, z_analytic))
+
+tolerance=0.001
+print("tolerance = "+ str(tolerance *100) + '%')
+
+diff_x=np.abs((x[0]-x_analytic)/x_analytic)
+diff_y=np.abs((y[0]-y_analytic)/y_analytic)
+
+print("percentage error for x = %5.4f %%" %(diff_x *100))
+print("percentage error for y = %5.4f %%" %(diff_y *100))
+
+assert (diff_x < tolerance) and (diff_y < tolerance) and (np.abs(z[0]) < tolerance), 'Test point_of_contact did not pass'

Examples/Tests/point_of_contact_EB/inputs_3d

@@ -0,0 +1,51 @@
+amr.max_level = 0
+
+max_step = 3
+
+
+amr.n_cell = 64 64 64
+amr.blocking_factor = 8
+amr.max_grid_size = 128
+
+geometry.dims = 3
+geometry.prob_lo     = -0.26 -0.26 -0.26
+geometry.prob_hi     = 0.26 0.26 0.26
+
+boundary.field_lo = pec pec pec
+boundary.field_hi = pec pec pec
+boundary.potential_lo_x = 0
+boundary.potential_hi_x = 0
+boundary.potential_lo_y = 0
+boundary.potential_hi_y = 0
+boundary.potential_lo_z = 0
+boundary.potential_hi_z = 0
+
+warpx.const_dt = 1.0e-10
+warpx.eb_implicit_function = "-(x**2+y**2+z**2-0.2**2)"
+
+# Do not evolve the E and B fields
+algo.maxwell_solver = none
+algo.field_gathering = momentum-conserving
+algo.particle_shape = 1
+
+particles.species_names = electron
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "SingleParticle"
+electron.single_particle_pos = -0.25 0 0
+electron.single_particle_u = 1 0.5 0.0
+electron.single_particle_weight = 1
+
+
+
+electron.save_particles_at_eb = 1
+
+diagnostics.diags_names = diag1 diag2
+
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex
+diag1.format = openpmd
+
+diag2.diag_type = BoundaryScraping
+diag2.format = openpmd

Examples/Tests/point_of_contact_EB/inputs_rz

@@ -0,0 +1,53 @@
+amr.max_level = 1
+warpx.fine_tag_lo = 0.0  -0.5
+warpx.fine_tag_hi = 0.25  0.0
+
+max_step = 3
+
+
+amr.n_cell = 64 64
+amr.blocking_factor = 8
+amr.max_grid_size = 128
+
+geometry.dims = RZ
+geometry.prob_lo     = 0.0 -0.5
+geometry.prob_hi     = 0.5  0.5
+
+boundary.field_lo = none periodic
+boundary.field_hi = pec periodic
+boundary.potential_lo_x = 0
+boundary.potential_hi_x = 0
+boundary.potential_lo_y = 0
+boundary.potential_hi_y = 0
+boundary.potential_lo_z = 0
+boundary.potential_hi_z = 0
+
+warpx.const_dt = 1.0e-10
+warpx.eb_implicit_function = "-(x**2 -0.2**2)"
+
+# Do not evolve the E and B fields
+algo.maxwell_solver = none
+algo.field_gathering = momentum-conserving
+algo.particle_shape = 1
+
+particles.species_names = electron
+electron.charge = -q_e
+electron.mass = m_e
+electron.injection_style = "SingleParticle"
+electron.single_particle_pos = -0.25 0 0
+electron.single_particle_u = 1 0.5 0.0
+electron.single_particle_weight = 1
+
+
+
+electron.save_particles_at_eb = 1
+
+diagnostics.diags_names = diag1 diag2
+
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Er
+diag1.format = openpmd
+
+diag2.diag_type = BoundaryScraping
+diag2.format = openpmd

Python/pywarpx/particle_containers.py

@@ -772,9 +772,9 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level)
                 form x/y/z_hi/lo or eb.
 
             comp_name      : str
-                The component of the array data that will be returned. If
-                "step_scraped" the special attribute holding the timestep at
-                which a particle was scraped will be returned.
+                The component of the array data that will be returned.
+                "x", "y", "z", "ux", "uy", "uz", "w"
+                "timestamp", "nx", "ny", "nz"
 
             level          : int
                 Which AMR level to retrieve scraped particle data from.
@@ -784,21 +784,13 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level)
         part_container = self.particle_buffer.get_particle_container(
             species_name, self._get_boundary_number(boundary)
         )
+        comp_idx = part_container.get_comp_index(comp_name)
         data_array = []
-        if comp_name == 'step_scraped':
-            # the step scraped is always the final integer component
-            comp_idx = part_container.num_int_comps - 1
-            for ii, pti in enumerate(libwarpx.libwarpx_so.BoundaryBufferParIter(part_container, level)):
-                soa = pti.soa()
-                data_array.append(xp.array(soa.GetIntData(comp_idx), copy=False))
-        else:
-            container_wrapper = ParticleContainerWrapper(species_name)
-            comp_idx = container_wrapper.particle_container.get_comp_index(comp_name)
-            for ii, pti in enumerate(libwarpx.libwarpx_so.BoundaryBufferParIter(part_container, level)):
-                soa = pti.soa()
-                data_array.append(xp.array(soa.GetRealData(comp_idx), copy=False))
-        return data_array
 
+        for ii, pti in enumerate(libwarpx.libwarpx_so.BoundaryBufferParIter(part_container, level)):
+            soa = pti.soa()
+            data_array.append(xp.array(soa.GetRealData(comp_idx), copy=False))
+        return data_array
 
     def clear_buffer(self):
         '''

Regression/Checksum/benchmarks_json/Point_of_contact_EB_3d.json

@@ -0,0 +1,14 @@
+{
+  "electron": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 0.0,
+    "particle_position_y": 0.0,
+    "particle_position_z": 0.0,
+    "particle_weight": 0.0
+  },
+  "lev=0": {
+    "Ex": 0.0
+  }
+}

Regression/Checksum/benchmarks_json/Point_of_contact_EB_rz.json

@@ -0,0 +1,19 @@
+{
+  "lev=0": {
+    "Er": 0.0
+  },
+  "lev=1": {
+    "Er": 0.0
+  },
+  "electron": {
+    "particle_position/theta": 0.0,
+    "particle_positionOffset/y": 0.0,
+    "particle_theta": 0.0,
+    "particle_position_x": 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/WarpX-tests.ini

@@ -4498,6 +4498,45 @@ doVis = 0
 outputFile = BeamBeamCollision_plt
 analysisRoutine = Examples/analysis_default_openpmd_regression.py
 
+[Point_of_contact_EB_3d]
+buildDir = .
+inputFile = Examples/Tests/point_of_contact_EB/inputs_3d
+runtime_params =
+dim = 3
+addToCompileString = USE_EB=TRUE
+cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_EB=ON
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 0
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+particleTypes = electrons
+outputFile = Point_of_contact_EB_3d_plt
+analysisRoutine = Examples/Tests/point_of_contact_EB/analysis.py
+
+[Point_of_contact_EB_rz]
+buildDir = .
+inputFile = Examples/Tests/point_of_contact_EB/inputs_rz
+runtime_params =
+dim = 2
+addToCompileString = USE_RZ=TRUE USE_EB=TRUE
+cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_EB=ON
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 0
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+particleTypes = electrons
+outputFile = Point_of_contact_EB_rz_plt
+analysisRoutine = Examples/Tests/point_of_contact_EB/analysis_rz.py
+
+
 [ImplicitPicard_1d]
 buildDir = .
 inputFile = Examples/Tests/Implicit/inputs_1d

Source/Diagnostics/WarpXOpenPMD.cpp

@@ -231,7 +231,7 @@ namespace detail
             //       and I/O in Cartesian.
             //       Other attributes like momentum are consequently
             //       stored in x,y,z internally.
-            positionComponents = vs{"x", "y", "z"};
+            positionComponents = vs{"x", "theta", "z"};
 #elif defined(WARPX_DIM_3D)
             positionComponents = vs{"x", "y", "z"};
 #else