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 3de88f3b3cb..d5e0db7213c 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,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"}, diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index b64c6d4b1fa..67fe9443f27 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -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 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("ionization") != 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(); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 5a3c925e9bd..bcb9e24240d 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -210,6 +210,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 e40a4e9822c..63373cce2b9 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -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; @@ -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(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) { @@ -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(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 473199a6b21..d46b7b56bac 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -85,7 +85,17 @@ public: 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. + // 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(num_added); tile_products[i]->resize(products_np[i] + num_added); @@ -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 @@ -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 @@ -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(p_pair_indices_2[i]), + static_cast(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(p_pair_indices_2[i]), + static_cast(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]; @@ -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((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); } else { amrex::Abort("Unknown scattering process."); @@ -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. diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp index de8de9b505d..43e3351451e 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp @@ -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 scattering_process_names; + pp_collision_name.queryarr("scattering_processes", scattering_process_names); + for (const auto& scattering_process : scattering_process_names) { + if (scattering_process.find("ionization") != 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: + // 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