Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[WIP] Add support for electron impact ionization in binary collisions (DSMC module) #5388

Draft
wants to merge 12 commits into
base: development
Choose a base branch
from
Draft

[WIP] Add support for electron impact ionization in binary collisions (DSMC module) #5388

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
7890b50
Activate DMSC for ionization
RemiLehe Oct 10, 2024
9b11492
Read name of the ionization species
RemiLehe Oct 11, 2024
1ac2c44
Create ionization executors
RemiLehe Oct 11, 2024
014fded
Move ionization into the same vector as the other processes
RemiLehe Oct 11, 2024
202fab0
Update number of collisions
RemiLehe Oct 11, 2024
f6a3d80
Fix test: use electrons
RemiLehe Oct 11, 2024
b291097
Attempt to generalize DSMC
RemiLehe Oct 11, 2024
12a5de4
Outline implementation
RemiLehe Oct 11, 2024
b5ebe92
Continue outline
RemiLehe Oct 11, 2024
c024977
Implement creation of products of ionization
RemiLehe Oct 11, 2024
cfc9a9b
Modify velocity of the electron
RemiLehe Oct 11, 2024
d29b23f
Fix compilation errors
RemiLehe Oct 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 33 additions & 25 deletions Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -268,32 +268,40 @@ 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:
electron_colls = picmi.DSMCCollisions(
name="coll_elec",
species=[self.electrons, self.neutrals],
ndt=5,
scattering_processes=electron_scattering_processes,
)
else:
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=electron_scattering_processes,
)

ion_scattering_processes = {
"elastic": {"cross_section": cross_sec_direc + "ion_scattering.dat"},
Expand Down
15 changes: 15 additions & 0 deletions Source/Particles/Collision/BinaryCollision/BinaryCollision.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,21 @@ public:
// Therefore, we insert the colliding species at the beginning of `m_product_species`
if (collision_type == CollisionType::DSMC) {
m_product_species.insert( m_product_species.begin(), m_species_names.begin(), m_species_names.end() );
// Check if ionization is one of the scattering processes
// TODO: This is reused in several places ; we should put it in a separate function
amrex::Vector<std::string> scattering_process_names;
bool ionization_flag = false;
pp_collision_name.queryarr("scattering_processes", scattering_process_names);
for (const auto& scattering_process : scattering_process_names) {
if (scattering_process.find("excitation") != std::string::npos) {
ionization_flag = true;
}
}
if (ionization_flag) {
std::string ionization_species;
pp_collision_name.get("ionization_species", ionization_species);
m_product_species.push_back(ionization_species);
}
}
m_have_product_species = !m_product_species.empty();

Expand Down
10 changes: 7 additions & 3 deletions Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,11 +65,15 @@ void CollisionPairFilter (const amrex::ParticleReal u1x, const amrex::ParticleRe

// Evaluate the cross-section for each scattering process to determine
// the total collision probability.
// TODO: Why only 4 previously?
// Note: this is brittle because it is not AMREX_ALWAYS_ASSERT_WITH_MESSAGE
// In practice, a user can specify e.g. 10 scattering processes and not get an error.
// We should check this in the constructor of DSMCFunc ahead of time.
AMREX_ASSERT_WITH_MESSAGE(
(process_count < 4), "Too many scattering processes in DSMC routine."
(process_count < 5), "Too many scattering processes in DSMC routine."
);
int coll_type[4] = {0, 0, 0, 0};
amrex::ParticleReal sigma_sums[4] = {0._prt, 0._prt, 0._prt, 0._prt};
int coll_type[5] = {0, 0, 0, 0, 0};
amrex::ParticleReal sigma_sums[5] = {0._prt, 0._prt, 0._prt, 0._prt, 0._prt};
for (int ii = 0; ii < process_count; ii++) {
auto const& scattering_process = scattering_processes[ii];
coll_type[ii] = int(scattering_process.m_type);
Expand Down
1 change: 1 addition & 0 deletions Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,6 +206,7 @@ public:
private:
amrex::Vector<ScatteringProcess> m_scattering_processes;
amrex::Gpu::DeviceVector<ScatteringProcess::Executor> m_scattering_processes_exe;

bool m_isSameSpecies;

Executor m_exe;
Expand Down
21 changes: 16 additions & 5 deletions Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ 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;
Expand All @@ -52,20 +53,30 @@ DSMCFunc::DSMCFunc (
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(process.type() != ScatteringProcessType::INVALID,
"Cannot add an unknown scattering process type");

// if the scattering process is ionization get the secondary species
// only one ionization process is supported
if (process.type() == ScatteringProcessType::IONIZATION) {
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!ionization_flag,
"Background MCC only supports a single ionization process");
ionization_flag = true;

// TODO: Add a check that the first species is the electron species
// (This should be done for impact ionization with MCC too)
// 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<int>(m_scattering_processes.size());

// Store ScatteringProcess::Executor(s).
#ifdef AMREX_USE_GPU
amrex::Gpu::HostVector<ScatteringProcess::Executor> 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) {
Expand All @@ -75,6 +86,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<int>(m_scattering_processes_exe.size());
m_exe.m_isSameSpecies = m_isSameSpecies;
}
65 changes: 60 additions & 5 deletions Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,17 @@ public:
amrex::Vector<int> 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.
// How many particles of product species i are potentially created.
// ("potentially" because the actual number of particles created depends on the
// type of reaction in the DSMC collision, which is stored in the mask. For
// instance, an IONIZATION reaction will create 2 electrons (the scattered
// electron, which will be split off from the incident electron, and the
// ionized electron), while an ELASTIC reaction will only 1 electron (the
// scattered electron). Therefore, here we allocate enough memory to store
// all potential particles that can be created (see also the constructor of
// SplitAndScatterFunc, where `m_num_products_host` is set), but set an
// invalid ID for the particles that should not actually be created, so that
// they are removed in the subsequent call to ReDistribute.
const index_type num_added = total * m_num_products_host[i];
num_added_vec[i] = static_cast<int>(num_added);
tile_products[i]->resize(products_np[i] + num_added);
Expand Down Expand Up @@ -120,16 +130,13 @@ public:
#endif

const int* AMREX_RESTRICT p_num_products_device = m_num_products_device.data();
const bool ionization_flag = m_ionization_flag;

amrex::ParallelForRNG(n_total_pairs,
[=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept
{
if (mask[i])
{
// 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);
// Make a copy of the particle from species 1
Expand All @@ -138,6 +145,7 @@ public:
// 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];

// TODO: when the reaction is ionization, this still creates a second particle of the target...
const auto product2_index = products_np_data[1] +
(p_offsets[i]*p_num_products_device[1] + 0);
// Make a copy of the particle from species 2
Expand All @@ -146,6 +154,32 @@ public:
// Set the weight of the new particles to p_pair_reaction_weight[i]
soa_products_data[1].m_rdata[PIdx::w][product2_index] = p_pair_reaction_weight[i];

if (ionization_flag) {
const auto product3_index = products_np_data[2] +
(p_offsets[i]*p_num_products_device[2] + 0);
if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
// The electron created from ionization is stored as the second particle in product 1
// Copy its properties from the target atoms, i.e. species 2 (momentum, position, etc.)
copy_species2[0](soa_products_data[0], soa_2, static_cast<int>(p_pair_indices_2[i]),
static_cast<int>(product1_index+1), engine);
// Set the weight of the new particles to p_pair_reaction_weight[i]
soa_products_data[0].m_rdata[PIdx::w][product1_index+1] = p_pair_reaction_weight[i];

// The ion created from ionization is stored as the first particle in product 3
// Copy its properties from the target atoms, i.e. species 2 (momentum, position, etc.)
copy_species2[2](soa_products_data[2], soa_2, static_cast<int>(p_pair_indices_2[i]),
static_cast<int>(product3_index), engine);
// Set the weight of the new particles to p_pair_reaction_weight[i]
soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i];
} else {
// Set IDs as invalid for slots that have been allocated for potential
// ionization events, so that the corresponding particles get removed in the next
// call to `ReDistribute` or `RemoveInvalidParticles`
soa_products_data[0].idcpu(product1_index+1) = amrex::ParticleIdCpus::Invalid;
soa_products_data[2].idcpu(product3_index) = amrex::ParticleIdCpus::Invalid;
}
}

// Set the child particle properties appropriately
auto& ux1 = soa_products_data[0].m_rdata[PIdx::ux][product1_index];
auto& uy1 = soa_products_data[0].m_rdata[PIdx::uy][product1_index];
Expand Down Expand Up @@ -202,6 +236,26 @@ public:
else {
amrex::Abort("Uneven mass charge-exchange not implemented yet.");
}
} else if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
// calculate kinetic energy of the incident electron
double E_coll;
double energy_cost = 0; // TODO: update this: get from scattering process
const amrex::ParticleReal u_coll2 = ux1*ux1 + uy1*uy1 + uz1*uz1;
ParticleUtils::getEnergy(u_coll2, m1, E_coll);
// each of the two electron (the scattered e) gets half the energy (could change this later)
const auto E_out = static_cast<amrex::ParticleReal>((E_coll - energy_cost) / 2.0_prt * PhysConst::q_e);
constexpr auto c2 = PhysConst::c * PhysConst::c;
const auto mc2 = m1*c2;
const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m1;

// isotropically scatter electrons
// - The incident electron
ParticleUtils::RandomizeVelocity(ux1, uy1, uz1, up, engine);
// - The newly created electron
auto& e_ux1 = soa_products_data[0].m_rdata[PIdx::ux][product1_index+1];
auto& e_uy1 = soa_products_data[0].m_rdata[PIdx::uy][product1_index+1];
auto& e_uz1 = soa_products_data[0].m_rdata[PIdx::uz][product1_index+1];
ParticleUtils::RandomizeVelocity(e_ux1, e_uy1, e_uz1, up, engine);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to do the transform back to labframe for e_u

}
else {
amrex::Abort("Unknown scattering process.");
Expand Down Expand Up @@ -244,6 +298,7 @@ public:
private:
// How many different type of species the collision produces
int m_num_product_species;
bool m_ionization_flag = false;
// 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.
Expand Down
33 changes: 26 additions & 7 deletions Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,17 +15,36 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name,
{
const amrex::ParmParse pp_collision_name(collision_name);

// Check if ionization is one of the scattering processes
amrex::Vector<std::string> scattering_process_names;
pp_collision_name.queryarr("scattering_processes", scattering_process_names);
for (const auto& scattering_process : scattering_process_names) {
if (scattering_process.find("excitation") != std::string::npos) {
m_ionization_flag = true;
}
}

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);
if (m_ionization_flag) {
// Product species include the ion
// TODO: as an alternative, we could use the runtime attribute `ionization_level` for this species
m_num_product_species = 3;
m_num_products_host.push_back(2); // electron species:
Copy link
Member Author

@RemiLehe RemiLehe Oct 11, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I attempted to fit impact ionization within the DSMC framework, but I now realize that this is making the book keeping complicated:

The amount of memory to be allocated (in SplitAndScatter) after the first pass (Executor, that checks which type of reaction is occurring for each pair) will depend on the type of reaction for each pair. (Because, if a pair does impact ionization, we need to allocate space for the scattered electron, the produced electron and the ionized species, whereas if a pair does one of the other reactions, we allocate space only for the scattered electron and the scattered target.) As a consequence, the offset (in SplitAndScatter will depend on the type of reaction.)

To simplify things, here I will always allocate space for the ionized species and the new electron, but if the reaction is not impact ionization, I will set the ID of the ionized species and the new electron to invalid.

// potentially 2 products per reaction: the scattered incoming electron, and the new electron from ionization
m_num_products_host.push_back(1);
m_num_products_host.push_back(1); // corresponds to the ionized species
} else {
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);
for (int i = 0; i < m_num_product_species; i++) {
m_num_products_device.push_back(m_num_products_host[i]);
}
#endif
}
else
Expand Down