diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 253f9ca0071..a8a9dd9d141 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -2070,7 +2070,7 @@ Details about the collision models can be found in the :ref:`theory section .product_species`` (`strings`) - Only for ``nuclearfusion``. The name(s) of the species in which to add + Only for ``dsmc`` and ``nuclearfusion``. The name(s) of the species in which to add the new macroparticles created by the reaction. * ``.ndt`` (`int`) optional @@ -2185,7 +2185,7 @@ Details about the collision models can be found in the :ref:`theory section ._energy`` (`float`) - Only for ``background_mcc``. If the scattering process is either + Only for ``dsmc`` and ``background_mcc``. If the scattering process is either ``excitationX`` or ``ionization`` the energy cost of that process must be given in eV. * ``.ionization_species`` (`float`) @@ -2193,6 +2193,10 @@ Details about the collision models can be found in the :ref:`theory section .ionization_target_species`` (`string`) + Only for ``dsmc`` with impact ionization. This specifies which one of the + colliding particles is ionized. + .. _running-cpp-parameters-numerics: Numerics and algorithms diff --git a/Examples/Physics_applications/capacitive_discharge/analysis_dsmc.py b/Examples/Physics_applications/capacitive_discharge/analysis_dsmc.py index cdaa6bed58f..e0ffd794946 100755 --- a/Examples/Physics_applications/capacitive_discharge/analysis_dsmc.py +++ b/Examples/Physics_applications/capacitive_discharge/analysis_dsmc.py @@ -6,39 +6,39 @@ # fmt: off ref_density = np.array([ - 1.27942709e+14, 2.23579371e+14, 2.55384387e+14, 2.55660663e+14, - 2.55830911e+14, 2.55814337e+14, 2.55798906e+14, 2.55744891e+14, - 2.55915585e+14, 2.56083194e+14, 2.55942354e+14, 2.55833026e+14, - 2.56036175e+14, 2.56234141e+14, 2.56196179e+14, 2.56146141e+14, - 2.56168022e+14, 2.56216909e+14, 2.56119961e+14, 2.56065167e+14, - 2.56194764e+14, 2.56416398e+14, 2.56465239e+14, 2.56234337e+14, - 2.56234503e+14, 2.56316003e+14, 2.56175023e+14, 2.56030269e+14, - 2.56189301e+14, 2.56286379e+14, 2.56130396e+14, 2.56295225e+14, - 2.56474082e+14, 2.56340375e+14, 2.56350864e+14, 2.56462330e+14, - 2.56469391e+14, 2.56412726e+14, 2.56241788e+14, 2.56355650e+14, - 2.56650599e+14, 2.56674748e+14, 2.56642480e+14, 2.56823508e+14, - 2.57025029e+14, 2.57110614e+14, 2.57042364e+14, 2.56950884e+14, - 2.57051822e+14, 2.56952148e+14, 2.56684016e+14, 2.56481130e+14, - 2.56277073e+14, 2.56065774e+14, 2.56190033e+14, 2.56411074e+14, - 2.56202418e+14, 2.56128368e+14, 2.56227002e+14, 2.56083004e+14, - 2.56056768e+14, 2.56343831e+14, 2.56443659e+14, 2.56280541e+14, - 2.56191572e+14, 2.56147304e+14, 2.56342794e+14, 2.56735473e+14, - 2.56994680e+14, 2.56901500e+14, 2.56527131e+14, 2.56490824e+14, - 2.56614730e+14, 2.56382744e+14, 2.56588214e+14, 2.57160270e+14, - 2.57230435e+14, 2.57116530e+14, 2.57065771e+14, 2.57236507e+14, - 2.57112865e+14, 2.56540177e+14, 2.56416828e+14, 2.56648954e+14, - 2.56625594e+14, 2.56411003e+14, 2.56523754e+14, 2.56841108e+14, - 2.56856368e+14, 2.56757912e+14, 2.56895134e+14, 2.57144419e+14, - 2.57001944e+14, 2.56371759e+14, 2.56179404e+14, 2.56541905e+14, - 2.56715727e+14, 2.56851681e+14, 2.57114458e+14, 2.57001739e+14, - 2.56825690e+14, 2.56879682e+14, 2.56699673e+14, 2.56532841e+14, - 2.56479582e+14, 2.56630989e+14, 2.56885996e+14, 2.56694637e+14, - 2.56250819e+14, 2.56045278e+14, 2.56366075e+14, 2.56693733e+14, - 2.56618530e+14, 2.56580918e+14, 2.56812781e+14, 2.56754216e+14, - 2.56444736e+14, 2.56473391e+14, 2.56538398e+14, 2.56626551e+14, - 2.56471950e+14, 2.56274969e+14, 2.56489423e+14, 2.56645266e+14, - 2.56611124e+14, 2.56344324e+14, 2.56244156e+14, 2.24183727e+14, - 1.27909856e+14 + 1.27939695e+14, 2.23589080e+14, 2.55400046e+14, 2.55652603e+14, + 2.55810704e+14, 2.55816145e+14, 2.55810457e+14, 2.55743643e+14, + 2.55908052e+14, 2.56076623e+14, 2.55948081e+14, 2.55841574e+14, + 2.56029524e+14, 2.56320511e+14, 2.56608595e+14, 2.56755504e+14, + 2.56699377e+14, 2.56700767e+14, 2.56497253e+14, 2.56481560e+14, + 2.56832303e+14, 2.57064841e+14, 2.57023000e+14, 2.56614315e+14, + 2.56368670e+14, 2.56370666e+14, 2.56227710e+14, 2.56240281e+14, + 2.56673842e+14, 2.56837209e+14, 2.56625623e+14, 2.56729845e+14, + 2.56975973e+14, 2.56801701e+14, 2.56491181e+14, 2.56516559e+14, + 2.56468688e+14, 2.56251727e+14, 2.56243466e+14, 2.56484137e+14, + 2.56637978e+14, 2.56448971e+14, 2.56140684e+14, 2.56117358e+14, + 2.56274706e+14, 2.56233588e+14, 2.56047578e+14, 2.56087060e+14, + 2.56365128e+14, 2.56357745e+14, 2.56269776e+14, 2.56419914e+14, + 2.56392856e+14, 2.56202826e+14, 2.56363244e+14, 2.56572545e+14, + 2.56351695e+14, 2.56393353e+14, 2.56759784e+14, 2.56767115e+14, + 2.56700246e+14, 2.56618056e+14, 2.56234915e+14, 2.56237788e+14, + 2.56606031e+14, 2.56520133e+14, 2.56316818e+14, 2.56184858e+14, + 2.56246807e+14, 2.56626394e+14, 2.56747253e+14, 2.56630112e+14, + 2.56518940e+14, 2.56358089e+14, 2.56249884e+14, 2.56271535e+14, + 2.56420396e+14, 2.56704340e+14, 2.56912250e+14, 2.56823163e+14, + 2.56694985e+14, 2.56822690e+14, 2.56736406e+14, 2.56438911e+14, + 2.56359312e+14, 2.56356028e+14, 2.56415261e+14, 2.56408702e+14, + 2.56267048e+14, 2.56274807e+14, 2.56494202e+14, 2.56789842e+14, + 2.56939719e+14, 2.56875327e+14, 2.56831776e+14, 2.56827482e+14, + 2.56698383e+14, 2.56712727e+14, 2.56879409e+14, 2.56629297e+14, + 2.56322165e+14, 2.56377317e+14, 2.56277894e+14, 2.56112364e+14, + 2.56171697e+14, 2.56370929e+14, 2.56855124e+14, 2.57621107e+14, + 2.57656000e+14, 2.56760729e+14, 2.56449741e+14, 2.56716250e+14, + 2.56721224e+14, 2.56506121e+14, 2.56236691e+14, 2.56270200e+14, + 2.56745053e+14, 2.56940581e+14, 2.56539958e+14, 2.56403313e+14, + 2.56600509e+14, 2.56776206e+14, 2.56884434e+14, 2.56755321e+14, + 2.56558818e+14, 2.56400159e+14, 2.56223931e+14, 2.23879043e+14, + 1.27601051e+14 ]) # fmt: on diff --git a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py index a03cf1954ad..0d1b9755af8 100644 --- a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py +++ b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py @@ -268,32 +268,55 @@ def setup_run(self): ####################################################################### cross_sec_direc = "../../../../warpx-data/MCC_cross_sections/He/" - electron_colls = picmi.MCCCollisions( - name="coll_elec", - species=self.electrons, - background_density=self.gas_density, - background_temperature=self.gas_temp, - background_mass=self.ions.mass, - ndt=self.mcc_subcycling_steps, - scattering_processes={ - "elastic": { - "cross_section": cross_sec_direc + "electron_scattering.dat" - }, - "excitation1": { - "cross_section": cross_sec_direc + "excitation_1.dat", - "energy": 19.82, - }, - "excitation2": { - "cross_section": cross_sec_direc + "excitation_2.dat", - "energy": 20.61, - }, - "ionization": { - "cross_section": cross_sec_direc + "ionization.dat", - "energy": 24.55, - "species": self.ions, - }, + + electron_scattering_processes = { + "elastic": {"cross_section": cross_sec_direc + "electron_scattering.dat"}, + "excitation1": { + "cross_section": cross_sec_direc + "excitation_1.dat", + "energy": 19.82, }, - ) + "excitation2": { + "cross_section": cross_sec_direc + "excitation_2.dat", + "energy": 20.61, + }, + "ionization": { + "cross_section": cross_sec_direc + "ionization.dat", + "energy": 24.55, + "species": self.ions, + }, + } + if self.dsmc: + ionization = {"ionization": electron_scattering_processes.pop("ionization")} + ionization["ionization"]["target_species"] = self.neutrals + ionization["ionization"].pop("species") + electron_colls_dsmc = picmi.DSMCCollisions( + name="coll_elec_dsmc", + species=[self.electrons, self.neutrals], + product_species=[self.ions, self.electrons], + ndt=4, + scattering_processes=ionization, + ) + electron_colls_mcc = picmi.MCCCollisions( + name="coll_elec", + species=self.electrons, + background_density=self.gas_density, + background_temperature=self.gas_temp, + background_mass=self.ions.mass, + ndt=self.mcc_subcycling_steps, + scattering_processes=electron_scattering_processes, + ) + electron_colls = [electron_colls_mcc, electron_colls_dsmc] + else: + electron_colls_mcc = picmi.MCCCollisions( + name="coll_elec", + species=self.electrons, + background_density=self.gas_density, + background_temperature=self.gas_temp, + background_mass=self.ions.mass, + ndt=self.mcc_subcycling_steps, + scattering_processes=electron_scattering_processes, + ) + electron_colls = [electron_colls_mcc] ion_scattering_processes = { "elastic": {"cross_section": cross_sec_direc + "ion_scattering.dat"}, @@ -316,6 +339,7 @@ def setup_run(self): ndt=self.mcc_subcycling_steps, scattering_processes=ion_scattering_processes, ) + ion_colls = [ion_colls] ####################################################################### # Initialize simulation # @@ -325,7 +349,7 @@ def setup_run(self): solver=self.solver, time_step_size=self.dt, max_steps=self.max_steps, - warpx_collisions=[electron_colls, ion_colls], + warpx_collisions=electron_colls + ion_colls, verbose=self.test, ) self.solver.sim = self.sim diff --git a/Examples/Tests/CMakeLists.txt b/Examples/Tests/CMakeLists.txt index 5ff1d4a9a70..8bb92fea2d3 100644 --- a/Examples/Tests/CMakeLists.txt +++ b/Examples/Tests/CMakeLists.txt @@ -26,6 +26,7 @@ add_subdirectory(gaussian_beam) add_subdirectory(implicit) add_subdirectory(initial_distribution) add_subdirectory(initial_plasma_profile) +add_subdirectory(ionization_dsmc) add_subdirectory(field_ionization) add_subdirectory(ion_stopping) add_subdirectory(langmuir) diff --git a/Examples/Tests/ionization_dsmc/CMakeLists.txt b/Examples/Tests/ionization_dsmc/CMakeLists.txt new file mode 100644 index 00000000000..db8d3cb7b2a --- /dev/null +++ b/Examples/Tests/ionization_dsmc/CMakeLists.txt @@ -0,0 +1,12 @@ +# Add tests (alphabetical order) ############################################## +# + +add_warpx_test( + test_3d_ionization_dsmc # name + 3 # dims + 2 # nprocs + inputs_test_3d_ionization_dsmc # inputs + "analysis_ionization_dsmc_3d.py" # analysis + "analysis_default_regression.py --path diags/diag1000250" # checksum + OFF # dependency +) diff --git a/Examples/Tests/ionization_dsmc/analysis_default_regression.py b/Examples/Tests/ionization_dsmc/analysis_default_regression.py new file mode 120000 index 00000000000..d8ce3fca419 --- /dev/null +++ b/Examples/Tests/ionization_dsmc/analysis_default_regression.py @@ -0,0 +1 @@ +../../analysis_default_regression.py \ No newline at end of file diff --git a/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py new file mode 100755 index 00000000000..294896a1761 --- /dev/null +++ b/Examples/Tests/ionization_dsmc/analysis_ionization_dsmc_3d.py @@ -0,0 +1,222 @@ +#!/usr/bin/env python3 +# DSMC ionization test script: +# - compares WarpX simulation results with theoretical model predictions. + +import matplotlib.pyplot as plt +import numpy as np +import tqdm +from openpmd_viewer import OpenPMDTimeSeries +from scipy.stats import qmc + +from pywarpx import picmi + +constants = picmi.constants + +ts = OpenPMDTimeSeries("diags/diag2/") + +q_e = constants.q_e +m_p = constants.m_p +m_e = constants.m_e +k_B = constants.kb +ep0 = constants.ep0 +clight = constants.c + +plasma_density = 1e14 +neutral_density = 1e20 +dt = 1e-9 +electron_temp = 10 +neutral_temp = 300 +max_steps = 250 +max_time = max_steps * dt + +L = [0.1] * 3 + +sigma_iz_file = "../../../../warpx-data/MCC_cross_sections/Xe/ionization.dat" +iz_data = np.loadtxt(sigma_iz_file) + +energy_eV = iz_data[:, 0] +sigma_m2 = iz_data[:, 1] +iz_energy = energy_eV[0] + + +def get_Te(ts): + T_e = [] + for iteration in tqdm.tqdm(ts.iterations): + ux, uy, uz = ts.get_particle( + ["ux", "uy", "uz"], species="electrons", iteration=iteration + ) + v_std_x = np.std(ux * clight) + v_std_y = np.std(uy * clight) + v_std_z = np.std(uz * clight) + v_std = (v_std_x + v_std_y + v_std_z) / 3 + T_e.append(m_e * v_std**2 / q_e) + return T_e + + +def get_density(ts): + number_data = np.genfromtxt("diags/counts.txt") + Te = get_Te(ts) + total_volume = L[0] * L[1] * L[2] + electron_weight = number_data[:, 8] + neutral_weight = number_data[:, 9] + ne = electron_weight / total_volume + nn = neutral_weight / total_volume + return [ne, nn, ne * Te] + + +def compute_rate_coefficients(temperatures_eV, energy_eV, sigma_m2, num_samples=1024): + """Integrate cross sections over maxwellian VDF to obtain reaction rate coefficients + Given electron energy in eV (`energy_eV`) and reaction cross sections at those energies (`sigma_m2`), + this function computes the reaction rate coefficient $k(T_e)$ for maxwellian electrons at + a provided list of electron temperatures `temperatures_eV`. + + The rate coefficient is given by + + $$ + k(T_e) = \int \sigma(E) E dE + $$ + + where the energy $E$ is drawn from a maxwellian distribution function with zero speed and temperature $T_e$. + We solve this using a quasi-monte carlo approach, by drawing a large number of low-discrepancy samples from + the appropriate distribution and obtaining the average of $\sigma(E) E$. + """ + thermal_speed_scale = np.sqrt(q_e / m_e) + k = np.zeros(temperatures_eV.size) + + # obtain low-discrepancy samples of normal dist + dist = qmc.MultivariateNormalQMC(np.zeros(3), np.eye(3)) + v = dist.random(num_samples) + + for i, T in enumerate(temperatures_eV): + # scale velocities to proper temperature + # compute energies corresponding to each sampled velocity vector + speed_squared = (v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2) * T + e = 0.5 * speed_squared + speed = np.sqrt(speed_squared) * thermal_speed_scale + # get cross section by interpolating on table + sigma = np.interp(e, energy_eV, sigma_m2, left=0) + k[i] = np.mean(sigma * speed) + return k + + +def rhs(state, params): + """Compute the right-hand side of ODE system that solves the global model described below. + The global model solves for the evolution of plasma density ($n_e$), neutral density ($n_n), + and electron temperature ($T_e$) in the presence of ionization. + The model equations consist of a continuity equation for electrons and neutrals, + combined with an energy equation for electrons. + + $$ + \frac{\partial n_e}{\partial t} = \dot{n} + \frac{\partial n_n}{\partial t} = -\dot{n} + \frac{3}{2}\frac{\partial n_e T_e}{\partial t} = -\dot{n} \epsilon_{iz}, + $$ + + where + + $$ + \dot{n} = n_n n_e k_{iz}(T_e), + $$ + + $k_iz$ is the ionization rate coefficient as a function of electron temperature in eV, and $\epsilon_{iz}$ is the ionization energy cost in eV. + """ + # unpack parameters + E_iz, Te_table, kiz_table = params + ne, nn, energy = state[0], state[1], state[2] + + # compute rhs + Te = energy / ne / 1.5 + rate_coeff = np.interp(Te, Te_table, kiz_table, left=0) + ndot = ne * nn * rate_coeff + f = np.empty(3) + + # fill in ionization rate + f[0] = ndot # d(ne)/dt + f[1] = -ndot # d(nn)/dt + f[2] = -ndot * iz_energy # -d(ne*eps) / dt + return f + + +def solve_theory_model(): + # integrate cross-section table to get rate coefficients + Te_table = np.linspace(0, 2 * electron_temp, 256) + kiz_table = compute_rate_coefficients( + Te_table, energy_eV, sigma_m2, num_samples=4096 + ) + + # set up system + num_steps = max_steps + 1 + state_vec = np.zeros((num_steps, 3)) + state_vec[0, :] = np.array( + [plasma_density, neutral_density, 1.5 * plasma_density * electron_temp] + ) + t = np.linspace(0, max_time, num_steps) + params = (iz_energy, Te_table, kiz_table) + + # solve the system (use RK2 integration) + for i in range(1, num_steps): + u = state_vec[i - 1, :] + k1 = rhs(u, params) + k2 = rhs(u + k1 * dt, params) + state_vec[i, :] = u + (k1 + k2) * dt / 2 + + # return result + ne = state_vec[:, 0] + nn = state_vec[:, 1] + Te = state_vec[:, 2] / (1.5 * ne) + return t, [ne, nn, ne * Te] + + +t_warpx = np.loadtxt("diags/counts.txt")[:, 1] +data_warpx = get_density(ts) + +t_theory, data_theory = solve_theory_model() + +fig, axs = plt.subplots(1, 3, figsize=(12, 4)) + +# Plot 1 +method = "dsmc" +labels = ["$n_e$ [m$^{-3}$]", "$n_n$ [m${-3}$]", "$n_e T_e$ [eVm$^{-3}$]"] +titles = ["Plasma density", "Neutral density", "Normalized electron temperature"] + +for i, (title, label, field_warpx, field_theory) in enumerate( + zip(titles, labels, data_warpx, data_theory) +): + axs[i].set_ylabel(label) + axs[i].set_title(title) + axs[i].plot(t_warpx, field_warpx, label="WarpX (" + method + ")") + axs[i].plot(t_theory, field_theory, label="theory", color="red", ls="--") + + axs[i].legend() +plt.tight_layout() +plt.savefig("ionization_dsmc_density_Te.png", dpi=150) + + +tolerances = [4e-2, 1e-6, 4e-2] + + +def check_tolerance(array, tolerance): + assert np.all(array <= tolerance), ( + f"Test did not pass: one or more elements exceed the tolerance of {tolerance}." + ) + print("All elements of are within the tolerance.") + + +plt.figure() +labels = [ + "Plasma density $(n_e$)", + "Neutral density $(n_n$)", + "Normalized electron temperature $(T_en_e$)", +] +for i, (label, field_warpx, field_theory, tolerance) in enumerate( + zip(labels, data_warpx, data_theory, tolerances) +): + relative_error = np.array( + abs((data_warpx[i] - data_theory[i][::5]) / data_theory[i][::5]) + ) + plt.plot(t_warpx, relative_error, label=label) + plt.ylabel("Relative error") + plt.xlabel("Time [s]") + plt.legend() + check_tolerance(relative_error, tolerance) +plt.savefig("./relative_error_density_Te.png", dpi=150) diff --git a/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc b/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc new file mode 100644 index 00000000000..aed4cf6c52d --- /dev/null +++ b/Examples/Tests/ionization_dsmc/inputs_test_3d_ionization_dsmc @@ -0,0 +1,116 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 250 +amr.n_cell = 8 8 8 +amr.max_level = 0 + +warpx.do_electrostatic = labframe +algo.particle_shape = 1 +amrex.verbose = 1 +geometry.coord_sys = 0 +################################# +############ CONSTANTS ############# +################################# +my_constants.Te = 10 +warpx.const_dt = 1e-09 + +################################# +####### GENERAL PARAMETERS ###### +################################# +geometry.dims = 3 +geometry.prob_hi = 0.1 0.1 0.1 +geometry.prob_lo = 0 0 0 +amr.max_grid_size = 8 + +################################# +###### BOUNDARY CONDITIONS ###### +################################# +geometry.is_periodic = 1 1 1 +boundary.field_hi = periodic periodic periodic +boundary.field_lo = periodic periodic periodic +boundary.particle_hi = periodic periodic periodic +boundary.particle_lo = periodic periodic periodic + +################################# +############ PLASMA ############# +################################# +particles.species_names = electrons ions neutrals + +electrons.charge = -q_e +electrons.density = 1e+14 +electrons.initialize_self_fields = 0 +electrons.injection_style = nuniformpercell +electrons.mass = m_e +electrons.momentum_distribution_type = gaussian +electrons.num_particles_per_cell_each_dim = 4 4 4 +electrons.profile = constant +electrons.ux_m = 0.0 +electrons.ux_th = sqrt(q_e * Te / m_e) / clight +electrons.uy_m = 0.0 +electrons.uy_th = sqrt(q_e * Te / m_e)/ clight +electrons.uz_m = 0.0 +electrons.uz_th = sqrt(q_e * Te / m_e)/ clight + +ions.charge = q_e +ions.density = 1e+14 +ions.initialize_self_fields = 0 +ions.injection_style = nuniformpercell +ions.mass = 2.196035502270312e-25 +ions.momentum_distribution_type = gaussian +ions.num_particles_per_cell_each_dim = 4 4 4 +ions.profile = constant +ions.ux_m = 0.0 +ions.ux_th = 4.5810168302300867e-07 +ions.uy_m = 0.0 +ions.uy_th = 4.5810168302300867e-07 +ions.uz_m = 0.0 +ions.uz_th = 4.5810168302300867e-07 + +neutrals.charge = 0 +neutrals.density = 1e+20 +neutrals.initialize_self_fields = 0 +neutrals.injection_style = nuniformpercell +neutrals.mass = 2.196035502270312e-25 +neutrals.momentum_distribution_type = gaussian +neutrals.num_particles_per_cell_each_dim = 4 4 4 +neutrals.profile = constant +neutrals.ux_m = 0.0 +neutrals.ux_th = 4.5810168302300867e-07 +neutrals.uy_m = 0.0 +neutrals.uy_th = 4.5810168302300867e-07 +neutrals.uz_m = 0.0 +neutrals.uz_th = 4.5810168302300867e-07 + +collisions.collision_names = coll_elec + +coll_elec.ionization_cross_section = ../../../../warpx-data/MCC_cross_sections/Xe/ionization.dat +coll_elec.ionization_energy = 12.1298431 +coll_elec.product_species = ions electrons +coll_elec.ionization_target_species = neutrals +coll_elec.ndt = 1 +coll_elec.scattering_processes = ionization +coll_elec.species = electrons neutrals +coll_elec.type = dsmc + +################################# +############ DIAGNOSTICS ############# +################################# +diagnostics.diags_names = diag1 diag2 +warpx.reduced_diags_names = counts +counts.intervals = 5 +counts.path = diags/ +counts.type = ParticleNumber + +# Diagnostics +diag1.intervals = 250 +diag1.diag_type = Full +diag1.electrons.variables = ux uy uz +diag1.neutrals.variables = ux uy uz +diag1.format = plotfile + +diag2.intervals = 5 +diag2.diag_type = Full +diag2.electrons.variables = ux uy uz +diag2.neutrals.variables = ux uy uz +diag2.format = openpmd diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index da673671953..90631cecf89 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -2481,14 +2481,21 @@ class DSMCCollisions(picmistandard.base._ClassWithInit): scattering_processes: dictionary The scattering process to use and any needed information + product_species: list + The species produced by collision processes (currently only ionization + products are supported). + ndt: integer, optional The collisions will be applied every "ndt" steps. Must be 1 or larger. """ - def __init__(self, name, species, scattering_processes, ndt=None, **kw): + def __init__( + self, name, species, scattering_processes, product_species=None, ndt=None, **kw + ): self.name = name self.species = species self.scattering_processes = scattering_processes + self.product_species = product_species self.ndt = ndt self.handle_init(kw) @@ -2497,12 +2504,16 @@ def collision_initialize_inputs(self): collision = pywarpx.Collisions.newcollision(self.name) collision.type = "dsmc" collision.species = [species.name for species in self.species] + if self.product_species is not None: + collision.product_species = [ + species.name for species in self.product_species + ] collision.ndt = self.ndt collision.scattering_processes = self.scattering_processes.keys() for process, kw in self.scattering_processes.items(): for key, val in kw.items(): - if key == "species": + if "species" in key: val = val.name collision.add_new_attr(process + "_" + key, val) diff --git a/Regression/Checksum/benchmarks_json/test_1d_dsmc_picmi.json b/Regression/Checksum/benchmarks_json/test_1d_dsmc_picmi.json index 62f915cf2e8..98872b583b8 100644 --- a/Regression/Checksum/benchmarks_json/test_1d_dsmc_picmi.json +++ b/Regression/Checksum/benchmarks_json/test_1d_dsmc_picmi.json @@ -1,27 +1,27 @@ { "lev=0": { - "rho_electrons": 0.004437338851654305, - "rho_he_ions": 0.005200276265886133 - }, - "he_ions": { - "particle_momentum_x": 2.768463746716725e-19, - "particle_momentum_y": 2.7585450668167785e-19, - "particle_momentum_z": 3.6189671443598533e-19, - "particle_position_x": 2201.408357891233, - "particle_weight": 17190472656250.002 + "rho_electrons": 0.004438215169514585, + "rho_he_ions": 0.005201960539530351 }, "electrons": { - "particle_momentum_x": 3.523554668287801e-20, - "particle_momentum_y": 3.515628626179393e-20, - "particle_momentum_z": 1.258711379033217e-19, - "particle_position_x": 2140.8168584833174, - "particle_weight": 14588988281250.002 + "particle_momentum_x": 3.533561243341172e-20, + "particle_momentum_y": 3.5386193532044873e-20, + "particle_momentum_z": 1.2559082152752288e-19, + "particle_position_x": 2141.145138810434, + "particle_weight": 14592390625000.002 + }, + "he_ions": { + "particle_momentum_x": 2.77081570738744e-19, + "particle_momentum_y": 2.757422176946563e-19, + "particle_momentum_z": 3.6216485183828697e-19, + "particle_position_x": 2201.938389152835, + "particle_weight": 17195707031250.002 }, "neutrals": { - "particle_momentum_x": 1.4054952479597137e-19, - "particle_momentum_y": 1.403311018061206e-19, - "particle_momentum_z": 1.411491089895956e-19, - "particle_position_x": 1119.82858839282, - "particle_weight": 6.4588e+19 + "particle_momentum_x": 1.4076569258763281e-19, + "particle_momentum_y": 1.4053910675178145e-19, + "particle_momentum_z": 1.4126450876863223e-19, + "particle_position_x": 1122.2997023272349, + "particle_weight": 6.4587999954199224e+19 } } diff --git a/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json new file mode 100644 index 00000000000..d9037876f89 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_3d_ionization_dsmc.json @@ -0,0 +1,38 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.0, + "Bz": 0.0, + "Ex": 22604.886154206244, + "Ey": 25871.77915984331, + "Ez": 24818.05831785292, + "jx": 2268.615138333634, + "jy": 326.15256830784614, + "jz": 288.7225344616736 + }, + "neutrals": { + "particle_momentum_x": 7.849329019484996e-19, + "particle_momentum_y": 7.871937065975134e-19, + "particle_momentum_z": 7.866670113934242e-19, + "particle_position_x": 1638.398241914082, + "particle_position_y": 1638.411011051806, + "particle_position_z": 1638.4009672441146 + }, + "ions": { + "particle_momentum_x": 1.1081540791730902e-18, + "particle_momentum_y": 1.1084000949226488e-18, + "particle_momentum_z": 1.1045654617165205e-18, + "particle_position_x": 2303.068003441504, + "particle_position_y": 2301.6312817994844, + "particle_position_z": 2301.7510230334133, + "particle_weight": 140463256835.93753 + }, + "electrons": { + "particle_momentum_x": 3.213371643156214e-20, + "particle_momentum_y": 3.2189558527727626e-20, + "particle_momentum_z": 2.7592636558840016e-20, + "particle_position_x": 2305.726186451743, + "particle_position_y": 2301.7018490320256, + "particle_position_z": 2303.7728642413595 + } +} diff --git a/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp b/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp index 8becd7d231a..5caa949055f 100644 --- a/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp +++ b/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp @@ -276,7 +276,7 @@ BackgroundMCCCollision::doCollisions (amrex::Real cur_time, amrex::Real dt, Mult } amrex::Print() << Utils::TextMsg::Info( - "Setting up collisions for " + m_species_names[0] + " with:\n" + "Setting up Monte-Carlo collisions for " + m_species_names[0] + " with:\n" + " total non-ionization collision probability: " + std::to_string(m_total_collision_prob) + "\n total ionization collision probability: " diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index b64c6d4b1fa..b73588a2b6e 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -102,10 +102,40 @@ public: const amrex::ParmParse pp_collision_name(collision_name); pp_collision_name.queryarr("product_species", m_product_species); - // if DSMC the colliding species are also product species - // Therefore, we insert the colliding species at the beginning of `m_product_species` + // If DSMC the colliding species are also product species. + // Therefore, we insert the colliding species at the beginning of `m_product_species`. if (collision_type == CollisionType::DSMC) { + // If the scattering process is ionization ensure that the + // explicitly specified "target" species, i.e., the species that + // undergoes ionization, is second in the species list for this + // collision set. The reason for this is that during the collision + // operation, an outgoing particle of the first species type will + // be created. + std::string target_species; + pp_collision_name.query("ionization_target_species", target_species); + if (!target_species.empty()) { + if (m_species_names[0] == target_species) { + std::swap(m_species_names[0], m_species_names[1]); + } else if (m_species_names[1] != target_species) { + WARPX_ABORT_WITH_MESSAGE("DSMC: Ionization target species, " + target_species + " must be one of the colliding species."); + } + } + m_product_species.insert( m_product_species.begin(), m_species_names.begin(), m_species_names.end() ); + + // During an ionization event, the non-ionizing target species could + // also be a product (for example, electron impact ionization producing another + // electron). We check that if this is the case, the product species is + // only listed once. The number of product particles is appropriately + // handled in `SplitAndScatterFunc` in that case. + if (m_product_species.size() > 2) { + if (m_product_species[2] == m_product_species[0]) { + m_product_species.erase(m_product_species.begin() + 2); + } + else if (m_product_species[3] == m_product_species[0]) { + m_product_species.erase(m_product_species.begin() + 3); + } + } } m_have_product_species = !m_product_species.empty(); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H index c5bd2e1cec6..fd4fd6fd66e 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H @@ -65,7 +65,6 @@ void CollisionPairFilter (const amrex::ParticleReal u1x, const amrex::ParticleRe // Evaluate the cross-section for each scattering process to determine // the total collision probability. - // The size of the arrays below is a compile-time constant (template parameter) // for performance reasons: it avoids dynamic memory allocation on the GPU. int coll_type[max_process_count] = {0}; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 6051aab1b59..6f5a6561b1e 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -211,6 +211,7 @@ public: private: amrex::Vector m_scattering_processes; amrex::Gpu::DeviceVector m_scattering_processes_exe; + bool m_isSameSpecies; Executor m_exe; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp index cf5f8de8d3c..b79741a536a 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -32,48 +32,54 @@ DSMCFunc::DSMCFunc ( // create a vector of ScatteringProcess objects from each scattering // process name + bool ionization_flag = false; for (const auto& scattering_process : scattering_process_names) { const std::string kw_cross_section = scattering_process + "_cross_section"; std::string cross_section_file; pp_collision_name.query(kw_cross_section.c_str(), cross_section_file); - // if the scattering process is excitation or ionization get the - // energy associated with that process + // if the scattering process is excitation, ionization or forward get + // the energy associated with that process + // (note that this allows forward scattering to be used both with and + // without a fixed energy loss) amrex::ParticleReal energy = 0._prt; if (scattering_process.find("excitation") != std::string::npos || - scattering_process.find("ionization") != std::string::npos) { + scattering_process.find("ionization") != std::string::npos || + scattering_process.find("forward") != std::string::npos ) { const std::string kw_energy = scattering_process + "_energy"; utils::parser::getWithParser( pp_collision_name, kw_energy.c_str(), energy); } - // if the scattering process is forward scattering get the energy - // associated with the process if it is given (this allows forward - // scattering to be used both with and without a fixed energy loss) - else if (scattering_process.find("forward") != std::string::npos) { - const std::string kw_energy = scattering_process + "_energy"; - utils::parser::queryWithParser( - pp_collision_name, kw_energy.c_str(), energy); - } ScatteringProcess process(scattering_process, cross_section_file, energy); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(process.type() != ScatteringProcessType::INVALID, "Cannot add an unknown scattering process type"); + // Only one ionization process is currently supported as part of a given + // collision set. + if (process.type() == ScatteringProcessType::IONIZATION) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !ionization_flag, + "DSMC only supports a single ionization process" + ); + ionization_flag = true; + + // And add a check that the ionization species has the same mass + // (and a positive charge), compared to the target species + } m_scattering_processes.push_back(std::move(process)); } - const int process_count = static_cast(m_scattering_processes.size()); - // Store ScatteringProcess::Executor(s). #ifdef AMREX_USE_GPU amrex::Gpu::HostVector h_scattering_processes_exe; for (auto const& p : m_scattering_processes) { h_scattering_processes_exe.push_back(p.executor()); } - m_scattering_processes_exe.resize(process_count); + m_scattering_processes_exe.resize(h_scattering_processes_exe.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_scattering_processes_exe.begin(), - h_scattering_processes_exe.end(), m_scattering_processes_exe.begin()); + h_scattering_processes_exe.end(), m_scattering_processes_exe.begin()); amrex::Gpu::streamSynchronize(); #else for (auto const& p : m_scattering_processes) { @@ -83,6 +89,6 @@ DSMCFunc::DSMCFunc ( // Link executor to appropriate ScatteringProcess executors m_exe.m_scattering_processes_data = m_scattering_processes_exe.data(); - m_exe.m_process_count = process_count; + m_exe.m_process_count = static_cast(m_scattering_processes_exe.size()); m_exe.m_isSameSpecies = m_isSameSpecies; } diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index e4b4d8a6a3a..59be9944077 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -67,28 +67,72 @@ public: { using namespace amrex::literals; + // Return a vector of zeros, indicating that for all the "product" species + // there were no new particles added. if (n_total_pairs == 0) { return amrex::Vector(m_num_product_species, 0); } - amrex::Gpu::DeviceVector offsets(n_total_pairs); - index_type* AMREX_RESTRICT offsets_data = offsets.data(); - const index_type* AMREX_RESTRICT p_offsets = offsets.dataPtr(); + // The following is used to calculate the appropriate offsets for + // non-product producing processes (i.e., non ionization processes). + // Note that a standard cummulative sum is not appropriate since the + // mask is also used to specify the type of collision and can therefore + // have values >1 + amrex::Gpu::DeviceVector no_product_offsets(n_total_pairs); + index_type* AMREX_RESTRICT no_product_offsets_data = no_product_offsets.data(); + const index_type* AMREX_RESTRICT no_product_p_offsets = no_product_offsets.dataPtr(); + auto const no_product_total = amrex::Scan::PrefixSum(n_total_pairs, + [=] AMREX_GPU_DEVICE (index_type i) -> index_type { + return ((mask[i] > 0) & (mask[i] != int(ScatteringProcessType::IONIZATION))) ? 1 : 0; + }, + [=] AMREX_GPU_DEVICE (index_type i, index_type s) { no_product_offsets_data[i] = s; }, + amrex::Scan::Type::exclusive, amrex::Scan::retSum + ); + + amrex::Vector num_added_vec(m_num_product_species, 0); + for (int i = 0; i < m_num_product_species; i++) + { + // Record the number of non product producing events lead to new + // particles for species1 and 2. Only 1 particle is created for + // each species (the piece that breaks off to have equal weight) + // particles. + num_added_vec[i] = static_cast(no_product_total); + } - // The following is used to calculate the appropriate offsets. Note that - // a standard cummulative sum is not appropriate since the mask is also - // used to specify the type of collision and can therefore have values >1 - auto const total = amrex::Scan::PrefixSum(n_total_pairs, - [=] AMREX_GPU_DEVICE (index_type i) -> index_type { return mask[i] ? 1 : 0; }, - [=] AMREX_GPU_DEVICE (index_type i, index_type s) { offsets_data[i] = s; }, + // The following is used to calculate the appropriate offsets for + // product producing processes (i.e., ionization). + // Note that a standard cummulative sum is not appropriate since the + // mask is also used to specify the type of collision and can therefore + // have values >1 + amrex::Gpu::DeviceVector with_product_offsets(n_total_pairs); + index_type* AMREX_RESTRICT with_product_offsets_data = with_product_offsets.data(); + const index_type* AMREX_RESTRICT with_product_p_offsets = with_product_offsets.dataPtr(); + auto const with_product_total = amrex::Scan::PrefixSum(n_total_pairs, + [=] AMREX_GPU_DEVICE (index_type i) -> index_type { + return (mask[i] == int(ScatteringProcessType::IONIZATION)) ? 1 : 0; + }, + [=] AMREX_GPU_DEVICE (index_type i, index_type s) { with_product_offsets_data[i] = s; }, amrex::Scan::Type::exclusive, amrex::Scan::retSum ); - amrex::Vector num_added_vec(m_num_product_species); for (int i = 0; i < m_num_product_species; i++) { - // How many particles of product species i are created. - const index_type num_added = total * m_num_products_host[i]; - num_added_vec[i] = static_cast(num_added); - tile_products[i]->resize(products_np[i] + num_added); + // Add the number of product producing events to the species involved + // in those processes. For the two colliding particles, if either is set to + // have just 1 copy in the products it indicates that that species is not a + // product of the product producing reaction (instead it is just tracked as + // an outgoing particle in non-product producing reactions), and therefore + // it does not count in the products. + int num_products = m_num_products_host[i]; + if ((i < 2) & (num_products == 1)) { + num_products = 0; + } + const index_type num_added = with_product_total * num_products; + num_added_vec[i] += static_cast(num_added); + } + + // resize the particle tiles to accomodate the new particles + for (int i = 0; i < m_num_product_species; i++) + { + tile_products[i]->resize(products_np[i] + num_added_vec[i]); } const auto soa_1 = ptile1.getParticleTileData(); @@ -119,27 +163,65 @@ public: const index_type* AMREX_RESTRICT products_np_data = products_np.data(); #endif - const int* AMREX_RESTRICT p_num_products_device = m_num_products_device.data(); + const int num_product_species = m_num_product_species; + const auto ionization_energy = m_ionization_energy; + // Store the list indices for ionization products, ensuring that + // the first product species is always an electron (which is assumed + // during the scattering operation). + // Also, get the starting index for the first ionization product (if ionization + // is present). If species1 is also a product, this would start the + // indexing for product particles after the particles created from + // fragmentation. + int ioniz_product1_list_index = 0, ioniz_product2_list_index = 0; + index_type ioniz_product1_offset = 0, ioniz_product2_offset = 0; + if (num_product_species == 3) { + if (pc_products[0]->getCharge() < 0.0) { + ioniz_product1_list_index = 0; + ioniz_product2_list_index = 2; + ioniz_product1_offset = products_np[0] + no_product_total + with_product_total; + ioniz_product2_offset = products_np[2]; + } else { + ioniz_product1_list_index = 2; + ioniz_product2_list_index = 0; + ioniz_product1_offset = products_np[2]; + ioniz_product2_offset = products_np[0] + no_product_total + with_product_total; + } + } else if (num_product_species == 4) { + if (pc_products[2]->getCharge() < 0.0) { + ioniz_product1_list_index = 2; + ioniz_product2_list_index = 3; + ioniz_product1_offset = products_np[2]; + ioniz_product2_offset = products_np[3]; + } else { + ioniz_product1_list_index = 3; + ioniz_product2_list_index = 2; + ioniz_product1_offset = products_np[3]; + ioniz_product2_offset = products_np[2]; + } + } + // Grab the masses of ionization products + amrex::ParticleReal m_ioniz_product1 = 0; + amrex::ParticleReal m_ioniz_product2 = 0; + if (num_product_species > 2) { + m_ioniz_product1 = pc_products[ioniz_product1_list_index]->getMass(); + m_ioniz_product2 = pc_products[ioniz_product2_list_index]->getMass(); + } + + // First perform all non-product producing collisions amrex::ParallelForRNG(n_total_pairs, [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept { - if (mask[i]) + if ((mask[i] > 0) & (mask[i] != int(ScatteringProcessType::IONIZATION))) { - // for now we ignore the possibility of having actual reaction - // products - only duplicating (splitting) of the colliding - // particles is supported. - - const auto product1_index = products_np_data[0] + - (p_offsets[i]*p_num_products_device[0] + 0); + const auto product1_index = products_np_data[0] + no_product_p_offsets[i]; // Make a copy of the particle from species 1 copy_species1[0](soa_products_data[0], soa_1, static_cast(p_pair_indices_1[i]), static_cast(product1_index), engine); // Set the weight of the new particles to p_pair_reaction_weight[i] soa_products_data[0].m_rdata[PIdx::w][product1_index] = p_pair_reaction_weight[i]; - const auto product2_index = products_np_data[1] + - (p_offsets[i]*p_num_products_device[1] + 0); + const auto product2_index = products_np_data[1] + no_product_p_offsets[i]; // Make a copy of the particle from species 2 copy_species2[1](soa_products_data[1], soa_2, static_cast(p_pair_indices_2[i]), static_cast(product2_index), engine); @@ -235,6 +317,184 @@ public: uy2 += uCOM_y; uz2 += uCOM_z; +#if (defined WARPX_DIM_RZ) + /* Undo the earlier velocity rotation. */ + amrex::ParticleReal const ux1buf_new = ux1; + ux1 = ux1buf_new*std::cos(-theta) - uy1*std::sin(-theta); + uy1 = ux1buf_new*std::sin(-theta) + uy1*std::cos(-theta); +#endif + } + + // Next perform all product producing collisions + else if (mask[i] == int(ScatteringProcessType::IONIZATION)) + { + const auto species1_index = products_np_data[0] + no_product_total + with_product_p_offsets[i]; + // Make a copy of the particle from species 1 + copy_species1[0](soa_products_data[0], soa_1, static_cast(p_pair_indices_1[i]), + static_cast(species1_index), engine); + // Set the weight of the new particles to p_pair_reaction_weight[i] + soa_products_data[0].m_rdata[PIdx::w][species1_index] = p_pair_reaction_weight[i]; + + // create a copy of the first product species at the location of species 2 + const auto product1_index = ioniz_product1_offset + with_product_p_offsets[i]; + copy_species1[ioniz_product1_list_index](soa_products_data[ioniz_product1_list_index], soa_2, static_cast(p_pair_indices_2[i]), + static_cast(product1_index), engine); + // Set the weight of the new particle to p_pair_reaction_weight[i] + soa_products_data[ioniz_product1_list_index].m_rdata[PIdx::w][product1_index] = p_pair_reaction_weight[i]; + + // create a copy of the other product species at the location of species 2 + const auto product2_index = ioniz_product2_offset + with_product_p_offsets[i]; + copy_species1[ioniz_product2_list_index](soa_products_data[ioniz_product2_list_index], soa_2, static_cast(p_pair_indices_2[i]), + static_cast(product2_index), engine); + // Set the weight of the new particle to p_pair_reaction_weight[i] + soa_products_data[ioniz_product2_list_index].m_rdata[PIdx::w][product2_index] = p_pair_reaction_weight[i]; + + // Grab the colliding particle velocities to calculate the COM + // Note that the two product particles currently have the same + // velocity as the "target" particle + auto& ux1 = soa_products_data[0].m_rdata[PIdx::ux][species1_index]; + auto& uy1 = soa_products_data[0].m_rdata[PIdx::uy][species1_index]; + auto& uz1 = soa_products_data[0].m_rdata[PIdx::uz][species1_index]; + auto& ux_p1 = soa_products_data[ioniz_product1_list_index].m_rdata[PIdx::ux][product1_index]; + auto& uy_p1 = soa_products_data[ioniz_product1_list_index].m_rdata[PIdx::uy][product1_index]; + auto& uz_p1 = soa_products_data[ioniz_product1_list_index].m_rdata[PIdx::uz][product1_index]; + auto& ux_p2 = soa_products_data[ioniz_product2_list_index].m_rdata[PIdx::ux][product2_index]; + auto& uy_p2 = soa_products_data[ioniz_product2_list_index].m_rdata[PIdx::uy][product2_index]; + auto& uz_p2 = soa_products_data[ioniz_product2_list_index].m_rdata[PIdx::uz][product2_index]; + +#if (defined WARPX_DIM_RZ) + /* In RZ geometry, macroparticles can collide with other macroparticles + * in the same *cylindrical* cell. For this reason, collisions between macroparticles + * are actually not local in space. In this case, the underlying assumption is that + * particles within the same cylindrical cell represent a cylindrically-symmetry + * momentum distribution function. Therefore, here, we temporarily rotate the + * momentum of one of the macroparticles in agreement with this cylindrical symmetry. + * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation; + * there is a corresponding assert statement at initialization.) + */ + amrex::ParticleReal const theta = ( + soa_products_data[ioniz_product1_list_index].m_rdata[PIdx::theta][product1_index] + - soa_products_data[0].m_rdata[PIdx::theta][species1_index] + ); + amrex::ParticleReal const ux1buf = ux1; + ux1 = ux1buf*std::cos(theta) - uy1*std::sin(theta); + uy1 = ux1buf*std::sin(theta) + uy1*std::cos(theta); +#endif + + // for simplicity (for now) we assume non-relativistic particles + // and simply calculate the center-of-momentum velocity from the + // rest masses + auto const uCOM_x = (m1 * ux1 + m2 * ux_p2) / (m1 + m2); + auto const uCOM_y = (m1 * uy1 + m2 * uy_p2) / (m1 + m2); + auto const uCOM_z = (m1 * uz1 + m2 * uz_p2) / (m1 + m2); + + // transform to COM frame + ux1 -= uCOM_x; + uy1 -= uCOM_y; + uz1 -= uCOM_z; + ux_p1 -= uCOM_x; + uy_p1 -= uCOM_y; + uz_p1 -= uCOM_z; + ux_p2 -= uCOM_x; + uy_p2 -= uCOM_y; + uz_p2 -= uCOM_z; + + if (mask[i] == int(ScatteringProcessType::IONIZATION)) { + // calculate kinetic energy of the collision (in eV) + const amrex::ParticleReal E1 = ( + 0.5_prt * m1 * (ux1*ux1 + uy1*uy1 + uz1*uz1) / PhysConst::q_e + ); + const amrex::ParticleReal E2 = ( + 0.5_prt * m2 * (ux_p2*ux_p2 + uy_p2*uy_p2 + uz_p2*uz_p2) / PhysConst::q_e + ); + const amrex::ParticleReal E_coll = E1 + E2; + + // subtract the energy cost for ionization + const amrex::ParticleReal E_out = (E_coll - ionization_energy) * PhysConst::q_e; + + // Energy division after the ionization event is done as follows: + // The ion product energy is obtained from the target energy as + // E2_prime = min(E2 / E_coll * E_out, 0.5 * E_out) + // The energy division for the remaining two particles + // must be done such that velocity vectors exist with net + // zero linear momentum in the current frame. A sufficient + // condition for this is that E1_prime, E2_prime and E3_prime + // are valid edge lengths for a triangle - effectively that + // a ellipse can be drawn from the energy components. + // That ellipse has semi-major and semi-minor axis: + // a = (E_out - E2_prime) / 2.0 + // b = 0.5 * sqrt(E_out^2 - 2 * E_out * E2_prime) + // The energy components are found by randomly sampling an + // x value between -a and a, and finding the corresponding + // y value that falls on the ellipse: y^2 = b^2 - b^2/a^2 * x^2. + // The secondary electron's energy is then: + // E0_prime = sqrt(y^2 + (x - E2_prime/2)^2) + // and the final particle's is: + // E1_prime = E_out - E0_prime - E0_prime + + // The product ordering ensures that product 2 is the + // ion product. + const amrex::ParticleReal E2_prime = std::min(E2 / E_coll * E_out, 0.5_prt * E_out); + + // find ellipse semi-major and minor axis + const amrex::ParticleReal a = 0.5_prt * (E_out - E2_prime); + const amrex::ParticleReal b = 0.5_prt * std::sqrt(E_out*E_out - 2.0_prt * E_out * E2_prime); + + // sample random x value and calculate y + const amrex::ParticleReal x = (2._prt * amrex::Random(engine) - 1.0_prt) * a; + const amrex::ParticleReal y2 = b*b - b*b/(a*a) * x*x; + const amrex::ParticleReal E0_prime = std::sqrt(y2 + x*x - x*E2_prime + 0.25_prt*E2_prime*E2_prime); + const amrex::ParticleReal E1_prime = E_out - E0_prime - E2_prime; + + // Now that appropriate energies are set for each outgoing species + // the directions for the velocity vectors must be chosen such + // that the net linear momentum in the current frame is 0. + // This is achieved by arranging the momentum vectors in + // a triangle and finding the required angles between the vectors. + const amrex::ParticleReal p0 = std::sqrt(2.0_prt * m1 * E0_prime); + const amrex::ParticleReal p1 = std::sqrt(2.0_prt * m_ioniz_product1 * E1_prime); + const amrex::ParticleReal p2 = std::sqrt(2.0_prt * m_ioniz_product2 * E2_prime); + + const amrex::ParticleReal cos_alpha = (p0*p0 + p1*p1 - p2*p2) / (2.0_prt * p0 * p1); + const amrex::ParticleReal sin_alpha = std::sqrt(1.0_prt - cos_alpha*cos_alpha); + const amrex::ParticleReal cos_gamma = (p0*p0 + p2*p2 - p1*p1) / (2.0_prt * p0 * p2); + const amrex::ParticleReal sin_gamma = std::sqrt(1.0_prt - cos_gamma*cos_gamma); + + // choose random theta and phi values (orientation of the triangle) + const amrex::ParticleReal Theta = amrex::Random(engine) * 2.0_prt * MathConst::pi; + const amrex::ParticleReal phi = amrex::Random(engine) * MathConst::pi; + + const amrex::ParticleReal cos_theta = std::cos(Theta); + const amrex::ParticleReal sin_theta = std::sin(Theta); + const amrex::ParticleReal cos_phi = std::cos(phi); + const amrex::ParticleReal sin_phi = std::sin(phi); + + ux1 = p0 / m1 * cos_theta * cos_phi; + uy1 = p0 / m1 * cos_theta * sin_phi; + uz1 = -p0 / m1 * sin_theta; + + ux_p1 = p1 / m_ioniz_product1 * (-cos_alpha * cos_theta * cos_phi - sin_alpha * sin_phi); + uy_p1 = p1 / m_ioniz_product1 * (-cos_alpha * cos_theta * sin_phi + sin_alpha * cos_phi); + uz_p1 = p1 / m_ioniz_product1 * (cos_alpha * sin_theta); + + ux_p2 = p2 / m_ioniz_product2 * (-cos_gamma * cos_theta * cos_phi + sin_gamma * sin_phi); + uy_p2 = p2 / m_ioniz_product2 * (-cos_gamma * cos_theta * sin_phi - sin_gamma * cos_phi); + uz_p2 = p2 / m_ioniz_product2 * (cos_gamma * sin_theta); + } + else { + amrex::Abort("Unknown scattering process."); + } + // transform back to labframe + ux1 += uCOM_x; + uy1 += uCOM_y; + uz1 += uCOM_z; + ux_p1 += uCOM_x; + uy_p1 += uCOM_y; + uz_p1 += uCOM_z; + ux_p2 += uCOM_x; + uy_p2 += uCOM_y; + uz_p2 += uCOM_z; + #if (defined WARPX_DIM_RZ) /* Undo the earlier velocity rotation. */ amrex::ParticleReal const ux1buf_new = ux1; @@ -272,10 +532,10 @@ public: private: // How many different type of species the collision produces int m_num_product_species; + // If ionization collisions are included, what is the energy cost + amrex::ParticleReal m_ionization_energy = 0.0; // Vectors of size m_num_product_species storing how many particles of a given species are - // produced by a collision event. These vectors are duplicated (one version for host and one - // for device) which is necessary with GPUs but redundant on CPU. - amrex::Gpu::DeviceVector m_num_products_device; + // produced by a collision event. amrex::Gpu::HostVector m_num_products_host; CollisionType m_collision_type; }; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp index de8de9b505d..8c93697a5a2 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp @@ -13,31 +13,62 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name, MultiParticleContainer const * const mypc): m_collision_type{BinaryCollisionUtils::get_collision_type(collision_name, mypc)} { - const amrex::ParmParse pp_collision_name(collision_name); - if (m_collision_type == CollisionType::DSMC) { - // here we can add logic to deal with cases where products are created, - // for example with impact ionization - m_num_product_species = 2; - m_num_products_host.push_back(1); - m_num_products_host.push_back(1); -#ifndef AMREX_USE_GPU - // On CPU, the device vector can be filled immediately - m_num_products_device.push_back(1); - m_num_products_device.push_back(1); -#endif + const amrex::ParmParse pp_collision_name(collision_name); + + // Check if ionization is one of the scattering processes by querying + // for any specified product species (ionization is the only current + // DSMC process with products) + amrex::Vector product_species; + pp_collision_name.queryarr("product_species", product_species); + + const bool ionization_flag = (!product_species.empty()); + + // if ionization is one of the processes, check if one of the colliding + // species is also used as a product species + if (ionization_flag) { + // grab the colliding species + amrex::Vector colliding_species; + pp_collision_name.getarr("species", colliding_species); + // grab the target species (i.e., the species that looses an + // electron during the collision) + std::string target_species; + pp_collision_name.query("ionization_target_species", target_species); + // find the index of the non-target species (the one that could + // also be used as a product species) + int non_target_idx = 0; + if (colliding_species[0] == target_species) { + non_target_idx = 1; + } + + // check if the non-target species is in ``product_species`` + auto it = std::find(product_species.begin(), product_species.end(), colliding_species[non_target_idx]); + + 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 + } + + // get the ionization energy + pp_collision_name.get("ionization_energy", m_ionization_energy); + + } else { + m_num_product_species = 2; + m_num_products_host.push_back(1); + m_num_products_host.push_back(1); + } } else { WARPX_ABORT_WITH_MESSAGE("Unknown collision type in SplitAndScatterFunc"); } - -#ifdef AMREX_USE_GPU - m_num_products_device.resize(m_num_product_species); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_num_products_host.begin(), - m_num_products_host.end(), - m_num_products_device.begin()); - amrex::Gpu::streamSynchronize(); -#endif }