Skip to content

Commit c3aa425

Browse files
RemiLeheRemiLehe
authored
Handle impact ionization for cases where the non-target incident species is not in the products (#5852)
# Original issue The structure `m_num_products_host` was being initialized in a complicated way: - First filled up as ``` if (it != product_species.end()) { m_num_product_species = 3; m_num_products_host.push_back(2); // the non-target species m_num_products_host.push_back(1); // the target species m_num_products_host.push_back(1); // corresponds to whichever ionization product species1 is not (ion or electron) } else { m_num_product_species = 4; m_num_products_host.push_back(1); // the non-target species m_num_products_host.push_back(1); // the target species m_num_products_host.push_back(1); // first product species m_num_products_host.push_back(1); // second product species } ``` - then updated with: ``` if ((i < 2) & (num_products == 1)) { num_products = 0; } ``` This combination ends up begin incorrect in the case of impact ionization reactions involving 4 species. (e.g. $H + H_2 \rightarrow H_2 + H^+ + e^{-}$), since in this case, we need to allocate space in the particle array for the non-ionized species ($H_2$) to have a scattered macroparticle - for conservation of energy. This results in a `SegFault` in practice, when using impact ionization reactions involving 4 species (which was not automatically tested, before this PR). # Changes in this PR Since `m_num_products_host` is only used for product-producing DSMC reactions, this code can be simplified to: ``` if (it != product_species.end()) { m_num_product_species = 3; m_num_products_host.push_back(2); // the non-target species m_num_products_host.push_back(0); // the target species m_num_products_host.push_back(1); // corresponds to whichever ionization product species1 is not (ion or electron) } else { m_num_product_species = 4; m_num_products_host.push_back(1); // the non-target species m_num_products_host.push_back(0); // the target species m_num_products_host.push_back(1); // first product species m_num_products_host.push_back(1); // second product species } ``` and then we don't need: ``` if ((i < 2) & (num_products == 1)) { num_products = 0; } ``` # Automated test I added a test featuring the reaction $H + H_2 \rightarrow H_2 + H^+ + e^{-}$, where a beam of fast $H^+$ goes through a background gas of $H_2$.
1 parent 47af36e commit c3aa425

File tree

6 files changed

+220
-16
lines changed

6 files changed

+220
-16
lines changed

‎Examples/Tests/ionization_dsmc/CMakeLists.txt‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,3 +20,13 @@ add_warpx_test(
2020
"analysis_default_regression.py --path diags/diag1000250" # checksum
2121
OFF # dependency
2222
)
23+
24+
add_warpx_test(
25+
test_1d_ionization_neutral_dsmc # name
26+
1 # dims
27+
2 # nprocs
28+
inputs_test_1d_ionization_neutral_dsmc # inputs
29+
"analysis_ionization_dsmc_1d.py" # analysis
30+
"analysis_default_regression.py --path diags/diag" # checksum
31+
OFF # dependency
32+
)

‎Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_1d.py‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
#!/usr/bin/env python3
2+
3+
import matplotlib.pyplot as plt
4+
import numpy as np
5+
from openpmd_viewer import OpenPMDTimeSeries
6+
from scipy.constants import c
7+
8+
# Read in the cross-section data for comparison with theory
9+
cross_section_data = np.loadtxt(
10+
"../../../../warpx-data/MCC_cross_sections/H/H_on_H2_ionization.dat"
11+
)
12+
E_eV = 12e3
13+
E_com_eV = E_eV * 2.0 / 3
14+
# Find the cross-section at the center of mass energy
15+
cross_section = np.interp(E_com_eV, cross_section_data[:, 0], cross_section_data[:, 1])
16+
print(f"Cross-section at {E_com_eV} eV: {cross_section} m^2")
17+
18+
# Read in the data
19+
ts = OpenPMDTimeSeries("./diags/diag")
20+
21+
# Compute the beam flux for ions and neutral, as a function of z, and compute with theory
22+
iteration = 100
23+
sim_flux = {"Hneutral": [], "Hplus": []}
24+
theory_flux = {"Hneutral": [], "Hplus": []}
25+
26+
# Compute the simulation flux
27+
# Loop over species
28+
for species in ["Hneutral", "Hplus"]:
29+
z, w, uz = ts.get_particle(
30+
["z", "w", "uz"],
31+
species=species,
32+
iteration=iteration,
33+
select={"uz": [1e-3, None]},
34+
)
35+
w_binned, bins = np.histogram(z, bins=25, weights=w * uz)
36+
# Convert from number of particles per bin to beam current
37+
dz_bin = bins[1] - bins[0]
38+
sim_flux[species] = w_binned / dz_bin * c
39+
40+
# Compute the theoretical flux
41+
n = 1e21
42+
flux = 1e20
43+
zmax = 0.2
44+
z_th = bins[:-1]
45+
theory_flux["Hneutral"] = flux * np.exp(-z_th * n * cross_section)
46+
theory_flux["Hplus"] = flux * (1 - np.exp(-z_th * n * cross_section))
47+
48+
# Compare the fluxes
49+
assert np.allclose(sim_flux["Hneutral"], theory_flux["Hneutral"], atol=5e-2 * flux)
50+
assert np.allclose(sim_flux["Hplus"], theory_flux["Hplus"], atol=5e-2 * flux)
51+
52+
# Plot the computed fluxes
53+
plt.figure()
54+
for species, color in zip(["Hneutral", "Hplus"], ["b", "r"]):
55+
plt.plot(bins[:-1], sim_flux[species], color=color, label=species)
56+
plt.plot(z_th, theory_flux[species], color=color, ls=":")
57+
plt.legend(loc=0)
58+
plt.ylabel("Beam flux")
59+
plt.xlabel("z [m]")
60+
plt.savefig("Beam_fluxes.png")

‎Examples/Tests/ionization_dsmc/inputs_test_1d_ionization_neutral_dsmc‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
# In this test, a beam of fast H neutral goes through a background gas of H2
2+
# Due to the ionization reaction H2 + Hneutral -> H2 + Hplus + electrons,
3+
# the beam of fast neutral is progressively converted into a beam of fast H+ ions,
4+
# with a rate that depends on the cross section and background gas density.
5+
# In the analysis script, we then check that the population of particles in the
6+
# beam is consistent with this rate.
7+
8+
#############################
9+
# Define parameter of the
10+
my_constants.beam_energy_ev = 12e3 # [eV]
11+
my_constants.beam_beta = sqrt(2*beam_energy_ev*q_e/m_p)/clight # [-]
12+
my_constants.beam_flux = 1e20 # [m^-2 s^-1]
13+
my_constants.gas_density = 1e21 # [m^-3]
14+
my_constants.gas_length = 0.2 # m
15+
16+
#############################
17+
# Simulation time
18+
max_step = 100
19+
warpx.const_dt = gas_length/beam_beta/clight/32
20+
21+
#############################
22+
# Grid
23+
amr.n_cell = 32
24+
amr.max_level = 0
25+
amr.max_grid_size = 256
26+
geometry.dims = 1
27+
28+
geometry.prob_lo = 0.0
29+
geometry.prob_hi = gas_length
30+
31+
boundary.particle_lo = absorbing
32+
boundary.particle_hi = absorbing
33+
34+
#############################
35+
# Algorithms:
36+
algo.particle_shape = 1
37+
38+
# - Without space charge
39+
algo.maxwell_solver = none
40+
boundary.field_lo = pec
41+
boundary.field_hi = pec
42+
43+
#############################
44+
# Create species
45+
particles.species_names = Hplus Hneutral electrons H2
46+
47+
Hneutral.do_continuous_injection = 1
48+
Hneutral.injection_style = NFluxPerCell
49+
Hneutral.mass = m_p
50+
Hneutral.charge = 0
51+
Hneutral.surface_flux_pos = 0
52+
Hneutral.flux_normal_axis = z
53+
Hneutral.flux_direction = 1
54+
Hneutral.num_particles_per_cell = 800
55+
Hneutral.flux_profile = constant
56+
Hneutral.flux = beam_flux
57+
Hneutral.momentum_distribution_type = constant
58+
Hneutral.uz = beam_beta
59+
60+
Hplus.injection_style = none
61+
Hplus.mass = m_p
62+
Hplus.charge = 0
63+
64+
H2.injection_style = NRandomPerCell
65+
H2.mass = 2*m_p
66+
H2.charge = 0
67+
H2.num_particles_per_cell = 1600
68+
H2.profile = constant
69+
H2.density = gas_density
70+
H2.momentum_distribution_type = constant
71+
H2.ux = 0
72+
H2.uy = 0
73+
H2.uz = 0
74+
H2.zmin = 0
75+
H2.zmax = gas_length
76+
77+
electrons.injection_style = none
78+
electrons.mass = m_e
79+
electrons.charge = -q_e
80+
81+
#############################
82+
# Add diagnostics
83+
diagnostics.diags_names = diag
84+
diag.intervals = 50
85+
diag.fields_to_plot = rho
86+
diag.diag_type = Full
87+
diag.format = openpmd
88+
89+
#############################
90+
91+
# Setup charge exchange
92+
collisions.collision_names = neutral_to_ion
93+
94+
neutral_to_ion.type = dsmc
95+
neutral_to_ion.scattering_processes = ionization
96+
neutral_to_ion.species = H2 Hneutral
97+
neutral_to_ion.ionization_target_species = Hneutral
98+
neutral_to_ion.product_species = Hplus electrons
99+
neutral_to_ion.ionization_cross_section = ../../../../warpx-data/MCC_cross_sections/H/H_on_H2_ionization.dat
100+
neutral_to_ion.ionization_energy = 13.6 # eV

‎Regression/Checksum/benchmarks_json/test_1d_ionization_neutral_dsmc.json‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
{
2+
"lev=0": {
3+
"rho": 3.0120123448636258e-05
4+
},
5+
"H2": {
6+
"particle_position_x": 0.0,
7+
"particle_position_y": 0.0,
8+
"particle_position_z": 10206.150677880647,
9+
"particle_momentum_x": 1.004262286923891e-18,
10+
"particle_momentum_y": 1.0070028443208195e-18,
11+
"particle_momentum_z": 6.86406551132128e-17,
12+
"particle_weight": 1.999999908139092e+20
13+
},
14+
"Hneutral": {
15+
"particle_position_x": 0.0,
16+
"particle_position_y": 0.0,
17+
"particle_position_z": 889.8320165555774,
18+
"particle_momentum_x": 0.0,
19+
"particle_momentum_y": 0.0,
20+
"particle_momentum_z": 3.1109861365015515e-17,
21+
"particle_weight": 6320718815591.17
22+
},
23+
"Hplus": {
24+
"particle_position_x": 0.0,
25+
"particle_position_y": 0.0,
26+
"particle_position_z": 5077.795115967319,
27+
"particle_momentum_x": 9.69219028133845e-19,
28+
"particle_momentum_y": 9.698576498693891e-19,
29+
"particle_momentum_z": 3.426639687689534e-17,
30+
"particle_weight": 20939216697480.555
31+
},
32+
"electrons": {
33+
"particle_position_x": 0.0,
34+
"particle_position_y": 0.0,
35+
"particle_position_z": 240.78671120953402,
36+
"particle_momentum_x": 7.189843092211492e-20,
37+
"particle_momentum_y": 7.214563802138814e-20,
38+
"particle_momentum_z": 8.578072077537713e-21,
39+
"particle_weight": 1227354057123.8418
40+
}
41+
}

‎Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H‎

Lines changed: 3 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -115,16 +115,8 @@ public:
115115

116116
for (int i = 0; i < m_num_product_species; i++)
117117
{
118-
// Add the number of product producing events to the species involved
119-
// in those processes. For the two colliding particles, if either is set to
120-
// have just 1 copy in the products it indicates that that species is not a
121-
// product of the product producing reaction (instead it is just tracked as
122-
// an outgoing particle in non-product producing reactions), and therefore
123-
// it does not count in the products.
118+
// Add the number of product producing events to the species involved in those processes.
124119
int num_products = m_num_products_host[i];
125-
if ((i < 2) & (num_products == 1)) {
126-
num_products = 0;
127-
}
128120
const index_type num_added = with_product_total * num_products;
129121
num_added_vec[i] += static_cast<int>(num_added);
130122
}
@@ -200,6 +192,7 @@ public:
200192
ioniz_product2_offset = products_np[2];
201193
}
202194
}
195+
203196
// Grab the masses of ionization products
204197
amrex::ParticleReal m_ioniz_product1 = 0;
205198
amrex::ParticleReal m_ioniz_product2 = 0;
@@ -331,7 +324,7 @@ public:
331324
const auto species1_index = products_np_data[0] + no_product_total + with_product_p_offsets[i];
332325
// Make a copy of the particle from species 1
333326
copy_species1[0](soa_products_data[0], soa_1, static_cast<int>(p_pair_indices_1[i]),
334-
static_cast<int>(species1_index), engine);
327+
static_cast<int>(species1_index), engine);
335328
// Set the weight of the new particles to p_pair_reaction_weight[i]
336329
soa_products_data[0].m_rdata[PIdx::w][species1_index] = p_pair_reaction_weight[i];
337330

‎Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -47,13 +47,13 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
4747

4848
if (it != product_species.end()) {
4949
m_num_product_species = 3;
50-
m_num_products_host.push_back(2); // the non-target species
51-
m_num_products_host.push_back(1); // the target species
52-
m_num_products_host.push_back(1); // corresponds to whichever ionization product species1 is not (ion or electron)
50+
m_num_products_host.push_back(2); // the non-target species: one particle for the product ; one particle for the fact that the incident particle loses energy
51+
m_num_products_host.push_back(0); // the target species
52+
m_num_products_host.push_back(1); // the other product of the reaction
5353
} else {
5454
m_num_product_species = 4;
5555
m_num_products_host.push_back(1); // the non-target species
56-
m_num_products_host.push_back(1); // the target species
56+
m_num_products_host.push_back(0); // the target species
5757
m_num_products_host.push_back(1); // first product species
5858
m_num_products_host.push_back(1); // second product species
5959
}
@@ -63,8 +63,8 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
6363

6464
} else {
6565
m_num_product_species = 2;
66-
m_num_products_host.push_back(1);
67-
m_num_products_host.push_back(1);
66+
m_num_products_host.push_back(0);
67+
m_num_products_host.push_back(0);
6868
}
6969
}
7070
else

0 commit comments

Comments
 (0)