Accessing the particles that hit the boundaries
WarpX can automatically save the particles that hit the boundaries
(see save_particles_at_xlo/ylo/zlo, save_particles_at_xhi/yhi/zhi,
and save_particles_at_eb in Inputs: Parameter List).
This data can be accessed in Python via the ParticleBoundaryBufferWrapper object,
which can is initialized as shown below.
from pywarpx import particle_containers
buffer = particle_containers.ParticleBoundaryBufferWrapper()
The ParticleBoundaryBufferWrapper object provides the following methods to access the particle boundary buffer data:
- class pywarpx.particle_containers.ParticleBoundaryBufferWrapper
Wrapper around particle boundary buffer containers. This provides a convenient way to query data in the particle boundary buffer containers.
- clear_buffer()
Clear the buffer that holds the particles lost at the boundaries.
- get_particle_boundary_buffer(species_name, boundary, comp_name, level)
This returns a list of numpy or cupy arrays containing the particle array data for a species that has been scraped by a specific simulation boundary.
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 an example of a simple case of particle-boundary interaction (reflection).
- Parameters:
species_name (str) – The species name that the data will be returned for.
boundary (str) – The boundary from which to get the scraped particle data in the form x/y/z_hi/lo or eb.
comp_name (str) – The component of the array data that will be returned. “x”, “y”, “z”, “ux”, “uy”, “uz”, “w” “stepScraped”,”deltaTimeScraped”, if boundary=’eb’: “nx”, “ny”, “nz”
level (int) – Which AMR level to retrieve scraped particle data from.
- get_particle_boundary_buffer_size(species_name, boundary, local=False)
This returns the number of particles that have been scraped so far in the simulation from the specified boundary and of the specified species.
- Parameters:
species_name (str) – Return the number of scraped particles of this species
boundary (str) – The boundary from which to get the scraped particle data in the form x/y/z_hi/lo
local (bool) – Whether to only return the number of particles in the current processor’s buffer
- get_particle_scraped_this_step(species_name, boundary, comp_name, level)
This returns a list of numpy or cupy arrays containing the particle array data for particles that have been scraped at the current timestep, for a specific species and simulation boundary.
The data for the arrays is a view of the underlying boundary buffer in WarpX ; writing to these arrays will therefore also modify the underlying boundary buffer.
- Parameters:
species_name (str) – The species name that the data will be returned for.
boundary (str) – The boundary from which to get the scraped particle data in the form x/y/z_hi/lo or eb.
comp_name (str) – The component of the array data that will be returned. “x”, “y”, “z”, “ux”, “uy”, “uz”, “w” “stepScraped”,”deltaTimeScraped”, if boundary=’eb’: “nx”, “ny”, “nz”
level (int) – Which AMR level to retrieve scraped particle data from.
This can be used to implement custom processes that occur at the boundaries (e.g., secondary emission), as in the example below.
Full example
Examples/Tests/secondary_ion_emission/inputs_test_rz_secondary_ion_emission_picmi.py.#!/usr/bin/env python3
# This is the script that tests secondary ion emission when ions hit an embedded boundary
# with a specified secondary emission yield of delta_H = 0.4. Specifically, a callback
# function at each time step ensures that the correct number of secondary electrons is
# emitted when ions impact the embedded boundary, following the given secondary emission
# model defined in sigma_nescap function. This distribution depends on the ion's energy and
# suggests that for an ion incident with 1 keV energy, an average of 0.4 secondary
# electrons will be emitted.
# Simulation is initialized with four ions with i_dist distribution and spherical
# embedded boundary given by implicit function.
import numpy as np
from scipy.constants import e, elementary_charge, m_e, proton_mass
from pywarpx import callbacks, particle_containers, picmi
from pywarpx.LoadThirdParty import load_cupy
##########################
# numerics parameters
##########################
dt = 0.000000075
# --- Nb time steps
Te = 0.0259 # in eV
dist_th = np.sqrt(Te * elementary_charge / m_e)
max_steps = 3
diagnostic_interval = 1
# --- grid
nr = 64
nz = 64
rmin = 0.0
rmax = 2
zmin = -2
zmax = 2
delta_H = 0.4
E_HMax = 250
##########################
# 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=["none", "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
##########################
i_dist = picmi.ParticleListDistribution(
x=[
0.025,
0.0,
-0.1,
-0.14,
],
y=[0.0, 0.0, 0.0, 0],
z=[-0.26, -0.29, -0.25, -0.23],
ux=[0.18e6, 0.1e6, 0.15e6, 0.21e6],
uy=[0.0, 0.0, 0.0, 0.0],
uz=[8.00e5, 7.20e5, 6.40e5, 5.60e5],
weight=[1, 1, 1, 1],
)
electrons = picmi.Species(
particle_type="electron", # Specify the particle type
name="electrons", # Name of the species
)
ions = picmi.Species(
name="ions",
particle_type="proton",
charge=e,
initial_distribution=i_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"],
warpx_format="openpmd",
)
part_diag = picmi.ParticleDiagnostic(
name="diag1",
period=diagnostic_interval,
species=[ions, electrons],
warpx_format="openpmd",
)
##########################
# 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=[0, 0, 0], grid=grid),
)
sim.add_species(
ions,
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
##########################
xp, _ = load_cupy()
xp.random.seed(10025015)
def concat(list_of_arrays):
if len(list_of_arrays) == 0:
# Return a 1d array of size 0
return xp.empty(0)
else:
return xp.concatenate(list_of_arrays)
def sigma_nascap(energy_kEv, delta_H, E_HMax):
"""
Compute sigma_nascap for each element in the energy array using a loop.
Parameters:
- energy_kEv: ndarray or list, energy values in KeV
- delta_H: float, parameter for the formula
- E_HMax: float, parameter for the formula in KeV
Returns:
- numpy or cupy array, computed probability sigma_nascap
"""
energy_kEv = xp.asarray(energy_kEv)
# Loop through each energy value
sigma_nascap = xp.where(
energy_kEv > 0.0,
delta_H
* (E_HMax + 1.0)
/ (E_HMax * 1.0 + energy_kEv)
* xp.sqrt(energy_kEv / 1.0),
0.0,
)
return sigma_nascap
def to_numpy(arr):
if hasattr(arr, "get"):
return arr.get()
else:
return arr
def secondary_emission():
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)
n = buffer.get_particle_boundary_buffer_size("ions", "eb")
electrons = sim.particles.get("electrons")
if n != 0:
r = concat(buffer.get_particle_scraped_this_step("ions", "eb", "r", lev))
theta = concat(
buffer.get_particle_scraped_this_step("ions", "eb", "theta", lev)
)
z = concat(buffer.get_particle_scraped_this_step("ions", "eb", "z", lev))
x = r * xp.cos(theta) # from RZ coordinates to 3D coordinates
y = r * xp.sin(theta)
ux = concat(buffer.get_particle_scraped_this_step("ions", "eb", "ux", lev))
uy = concat(buffer.get_particle_scraped_this_step("ions", "eb", "uy", lev))
uz = concat(buffer.get_particle_scraped_this_step("ions", "eb", "uz", lev))
w = concat(buffer.get_particle_scraped_this_step("ions", "eb", "w", lev))
nx = concat(buffer.get_particle_scraped_this_step("ions", "eb", "nx", lev))
ny = concat(buffer.get_particle_scraped_this_step("ions", "eb", "ny", lev))
nz = concat(buffer.get_particle_scraped_this_step("ions", "eb", "nz", lev))
delta_t = concat(
buffer.get_particle_scraped_this_step("ions", "eb", "deltaTimeScraped", lev)
)
energy_ions = 0.5 * proton_mass * w * (ux**2 + uy**2 + uz**2)
energy_ions_in_kEv = energy_ions / (e * 1000)
sigma_nascap_ions = sigma_nascap(energy_ions_in_kEv, delta_H, E_HMax)
x = to_numpy(x)
y = to_numpy(y)
z = to_numpy(z)
w = to_numpy(w)
nx = to_numpy(nx)
ny = to_numpy(ny)
nz = to_numpy(nz)
delta_t = to_numpy(delta_t)
sigma_nascap_ions = to_numpy(sigma_nascap_ions)
xe = np.array([])
ye = np.array([])
ze = np.array([])
we = np.array([])
delta_te = np.array([])
uxe = np.array([])
uye = np.array([])
uze = np.array([])
# Loop over all ions that have been scraped in the last timestep
for i in range(0, len(w)):
sigma = sigma_nascap_ions[i]
# Ne_sec is number of the secondary electrons to be emitted
Ne_sec = int(sigma + np.random.uniform())
for _ in range(Ne_sec):
# Random thermal momenta distribution
ux_th = np.random.normal(0, dist_th)
uy_th = np.random.normal(0, dist_th)
uz_th = np.random.normal(0, dist_th)
un_th = nx[i] * ux_th + ny[i] * uy_th + nz[i] * uz_th
if un_th < 0:
ux_th_reflect = (
-2 * un_th * nx[i] + ux_th
) # for a "mirror reflection" u(sym)=-2(u.n)n+u
uy_th_reflect = -2 * un_th * ny[i] + uy_th
uz_th_reflect = -2 * un_th * nz[i] + uz_th
uxe = np.append(uxe, ux_th_reflect)
uye = np.append(uye, uy_th_reflect)
uze = np.append(uze, uz_th_reflect)
else:
uxe = np.append(uxe, ux_th)
uye = np.append(uye, uy_th)
uze = np.append(uze, uz_th)
# Also convert the position and weight arrays
xe = np.append(xe, x[i])
ye = np.append(ye, y[i])
ze = np.append(ze, z[i])
we = np.append(we, w[i])
delta_te = np.append(delta_te, delta_t[i])
electrons.add_particles(
x=xe + (dt - delta_te) * uxe,
y=ye + (dt - delta_te) * uye,
z=ze + (dt - delta_te) * uze,
ux=uxe,
uy=uye,
uz=uze,
w=we,
)
# using the new particle container modified at the last step
callbacks.installafterstep(secondary_emission)
##########################
# simulation run
##########################
sim.step(max_steps) # the whole process is done "max_steps" times