From 7890b50e6e1c4f52c3969b739c7260f68c9c877d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 16:50:46 -0700 Subject: [PATCH 01/14] Activate DMSC for ionization --- .../inputs_base_1d_picmi.py | 58 +++++++++++-------- 1 file changed, 33 insertions(+), 25 deletions(-) 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..5fd19bc3b4b 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.ions, 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"}, From 9b11492514db6934ea6430e1ee7cadbca95fbf48 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 17:02:33 -0700 Subject: [PATCH 02/14] Read name of the ionization species --- .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 3 +++ .../BinaryCollision/DSMC/DSMCFunc.cpp | 19 ++++++++++++++++++- 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index a692d2cbb9e..4c813ff1b5d 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -205,7 +205,10 @@ public: private: amrex::Vector m_scattering_processes; + amrex::Vector m_ionization_processes; amrex::Gpu::DeviceVector m_scattering_processes_exe; + amrex::Gpu::DeviceVector m_ionization_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..299b356404a 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -52,7 +52,24 @@ DSMCFunc::DSMCFunc ( WARPX_ALWAYS_ASSERT_WITH_MESSAGE(process.type() != ScatteringProcessType::INVALID, "Cannot add an unknown scattering process type"); - m_scattering_processes.push_back(std::move(process)); + // if the scattering process is ionization get the secondary species + // only one ionization process is supported, the vector + // m_ionization_processes is only used to make it simple to calculate + // the maximum collision frequency with the same function used for + // particle conserving processes + if (process.type() == ScatteringProcessType::IONIZATION) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!ionization_flag, + "Background MCC only supports a single ionization process"); + ionization_flag = true; + + std::string secondary_species; + pp_collision_name.get("ionization_species", secondary_species); + m_species_names.push_back(secondary_species); + + m_ionization_processes.push_back(std::move(process)); + } else { + m_scattering_processes.push_back(std::move(process)); + } } const int process_count = static_cast(m_scattering_processes.size()); From 1ac2c44c9a16e5c1ad5c4697fe43679ae10912e4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 17:23:36 -0700 Subject: [PATCH 03/14] Create ionization executors --- .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 3 +++ .../BinaryCollision/DSMC/DSMCFunc.cpp | 18 +++++++++++++----- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 4c813ff1b5d..1e296414b95 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -209,6 +209,9 @@ private: amrex::Gpu::DeviceVector m_scattering_processes_exe; amrex::Gpu::DeviceVector m_ionization_processes_exe; + bool ionization_flag = false; + + amrex::Vector m_species_names; // TODO: is this the right place to store this? 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 299b356404a..cdcbf61173b 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -72,26 +72,34 @@ DSMCFunc::DSMCFunc ( } } - 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; + amrex::Gpu::HostVector h_ionization_processes_exe; for (auto const& p : m_scattering_processes) { h_scattering_processes_exe.push_back(p.executor()); } - m_scattering_processes_exe.resize(process_count); + for (auto const& p : m_ionization_processes) { + h_ionization_processes_exe.push_back(p.executor()); + } + m_scattering_processes_exe.resize(h_scattering_processes_exe.size()); + m_ionization_processes_exe.resize(h_ionization_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::copyAsync(amrex::Gpu::hostToDevice, h_ionization_processes_exe.begin(), + h_ionization_processes_exe.end(), m_ionization_processes_exe.begin()); amrex::Gpu::streamSynchronize(); #else for (auto const& p : m_scattering_processes) { m_scattering_processes_exe.push_back(p.executor()); } + for (auto const& p : m_ionization_processes) { + m_ionization_processes_exe.push_back(p.executor()); + } #endif // 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; } From 014fdedbb3ca2fb4795176dc101329224094b0af Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 18:28:13 -0700 Subject: [PATCH 04/14] Move ionization into the same vector as the other processes --- .../DSMC/CollisionFilterFunc.H | 1 + .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 2 -- .../BinaryCollision/DSMC/DSMCFunc.cpp | 20 ++----------------- 3 files changed, 3 insertions(+), 20 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H index 46b228b049e..1283a086012 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H @@ -65,6 +65,7 @@ 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? AMREX_ASSERT_WITH_MESSAGE( (process_count < 4), "Too many scattering processes in DSMC routine." ); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 1e296414b95..191cf698df1 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -205,9 +205,7 @@ public: private: amrex::Vector m_scattering_processes; - amrex::Vector m_ionization_processes; amrex::Gpu::DeviceVector m_scattering_processes_exe; - amrex::Gpu::DeviceVector m_ionization_processes_exe; bool ionization_flag = false; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp index cdcbf61173b..e4d17a97e19 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -53,10 +53,7 @@ DSMCFunc::DSMCFunc ( "Cannot add an unknown scattering process type"); // if the scattering process is ionization get the secondary species - // only one ionization process is supported, the vector - // m_ionization_processes is only used to make it simple to calculate - // the maximum collision frequency with the same function used for - // particle conserving processes + // 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"); @@ -65,37 +62,24 @@ DSMCFunc::DSMCFunc ( std::string secondary_species; pp_collision_name.get("ionization_species", secondary_species); m_species_names.push_back(secondary_species); - - m_ionization_processes.push_back(std::move(process)); - } else { - m_scattering_processes.push_back(std::move(process)); } + m_scattering_processes.push_back(std::move(process)); } // Store ScatteringProcess::Executor(s). #ifdef AMREX_USE_GPU amrex::Gpu::HostVector h_scattering_processes_exe; - amrex::Gpu::HostVector h_ionization_processes_exe; for (auto const& p : m_scattering_processes) { h_scattering_processes_exe.push_back(p.executor()); } - for (auto const& p : m_ionization_processes) { - h_ionization_processes_exe.push_back(p.executor()); - } m_scattering_processes_exe.resize(h_scattering_processes_exe.size()); - m_ionization_processes_exe.resize(h_ionization_processes_exe.size()); amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_scattering_processes_exe.begin(), h_scattering_processes_exe.end(), m_scattering_processes_exe.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_ionization_processes_exe.begin(), - h_ionization_processes_exe.end(), m_ionization_processes_exe.begin()); amrex::Gpu::streamSynchronize(); #else for (auto const& p : m_scattering_processes) { m_scattering_processes_exe.push_back(p.executor()); } - for (auto const& p : m_ionization_processes) { - m_ionization_processes_exe.push_back(p.executor()); - } #endif // Link executor to appropriate ScatteringProcess executors From 202fab0c2b58ea7e22482917ce6f777a75022ca7 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 10 Oct 2024 18:44:27 -0700 Subject: [PATCH 05/14] Update number of collisions --- .../BinaryCollision/DSMC/CollisionFilterFunc.H | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H index 1283a086012..d05cc4f8c4f 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H @@ -65,12 +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? + // 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); From f6a3d809d8dcc66469db214d9d4a03e2f3a735b7 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 10:15:54 -0700 Subject: [PATCH 06/14] Fix test: use electrons --- .../capacitive_discharge/inputs_base_1d_picmi.py | 2 +- Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) 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 5fd19bc3b4b..d5e0db7213c 100644 --- a/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py +++ b/Examples/Physics_applications/capacitive_discharge/inputs_base_1d_picmi.py @@ -288,7 +288,7 @@ def setup_run(self): if self.dsmc: electron_colls = picmi.DSMCCollisions( name="coll_elec", - species=[self.ions, self.neutrals], + species=[self.electrons, self.neutrals], ndt=5, scattering_processes=electron_scattering_processes, ) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp index e4d17a97e19..3fa1c9a74b1 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.cpp @@ -59,6 +59,9 @@ DSMCFunc::DSMCFunc ( "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) + std::string secondary_species; pp_collision_name.get("ionization_species", secondary_species); m_species_names.push_back(secondary_species); From b291097fdc2ee982e5b9349c2fead5f99cb61b75 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 10:49:27 -0700 Subject: [PATCH 07/14] Attempt to generalize DSMC --- .../DSMC/SplitAndScatterFunc.cpp | 34 +++++++++++++++---- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp index de8de9b505d..d145475119e 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp @@ -15,17 +15,37 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name, { const amrex::ParmParse pp_collision_name(collision_name); + // Check if ionization is one of the scattering processes + bool ionization_flag = false; + 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("excitation") != std::string::npos) { + 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 (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 From 12a5de47c454deacf56b6a4ffcde59e5821623fc Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 11:37:09 -0700 Subject: [PATCH 08/14] Outline implementation --- .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 2 +- .../DSMC/SplitAndScatterFunc.H | 37 +++++++++++++++++-- .../DSMC/SplitAndScatterFunc.cpp | 3 +- 3 files changed, 35 insertions(+), 7 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 191cf698df1..d3ef57e4c92 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -207,7 +207,7 @@ private: amrex::Vector m_scattering_processes; amrex::Gpu::DeviceVector m_scattering_processes_exe; - bool ionization_flag = false; + bool ionization_flag = false; // Not sure we actually need to store this amrex::Vector m_species_names; // TODO: is this the right place to store this? bool m_isSameSpecies; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index 473199a6b21..d7167a96c87 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,15 +130,14 @@ 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. + // TODO: by default: set all IDs to be invalid const auto product1_index = products_np_data[0] + (p_offsets[i]*p_num_products_device[0] + 0); @@ -146,6 +155,23 @@ 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) { + // The electron created from ionization is stored as the second particle in product 1 + // Copy its properties from species 2 (momentum, position, etc.) + + + // The ion created from ionization is stored as the first particle in product 3 + // Copy its properties from species 2 (momentum, position, etc.) + const auto product3_index = products_np_data[2] + + (p_offsets[i]*p_num_products_device[2] + 0); + + // Make a copy of the particle from species 3 + copy_species2[2](soa_products_data[2], soa_3, static_cast(p_pair_indices_3[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]; + } + // 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 +228,8 @@ public: else { amrex::Abort("Uneven mass charge-exchange not implemented yet."); } + } else if (mask[i] == int(ScatteringProcessType::IONIZATION)) { + // Same as in doBackgroundIonization } else { amrex::Abort("Unknown scattering process."); @@ -244,6 +272,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 d145475119e..b28871bd375 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp @@ -16,12 +16,11 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name, const amrex::ParmParse pp_collision_name(collision_name); // Check if ionization is one of the scattering processes - bool ionization_flag = false; 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("excitation") != std::string::npos) { - ionization_flag = true; + m_ionization_flag = true; } } From b5ebe9248281a8e3584e29c8ee3f19c4b56ca6fa Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 11:43:36 -0700 Subject: [PATCH 09/14] Continue outline --- .../Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index d7167a96c87..f580df48bda 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -158,17 +158,18 @@ public: if (ionization_flag) { // The electron created from ionization is stored as the second particle in product 1 // Copy its properties from species 2 (momentum, position, etc.) + // TODO + // Set the weight of the new particles to p_pair_reaction_weight[i] // The ion created from ionization is stored as the first particle in product 3 // Copy its properties from species 2 (momentum, position, etc.) const auto product3_index = products_np_data[2] + (p_offsets[i]*p_num_products_device[2] + 0); + // TODO - // Make a copy of the particle from species 3 - copy_species2[2](soa_products_data[2], soa_3, static_cast(p_pair_indices_3[i]), - static_cast(product3_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+1] = p_pair_reaction_weight[i]; soa_products_data[2].m_rdata[PIdx::w][product3_index] = p_pair_reaction_weight[i]; } From c02497773ff9dab52e5b7c1fd71431a1616c16b5 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 12:29:30 -0700 Subject: [PATCH 10/14] Implement creation of products of ionization --- .../BinaryCollision/BinaryCollision.H | 15 ++++++++ .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 3 -- .../BinaryCollision/DSMC/DSMCFunc.cpp | 7 ++-- .../DSMC/SplitAndScatterFunc.H | 35 +++++++++++-------- .../DSMC/SplitAndScatterFunc.cpp | 2 +- 5 files changed, 40 insertions(+), 22 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index b64c6d4b1fa..a8cb3e4970a 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("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(); diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index d3ef57e4c92..a26bb599447 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -207,9 +207,6 @@ private: amrex::Vector m_scattering_processes; amrex::Gpu::DeviceVector m_scattering_processes_exe; - bool ionization_flag = false; // Not sure we actually need to store this - - amrex::Vector m_species_names; // TODO: is this the right place to store this? 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 3fa1c9a74b1..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; @@ -61,10 +62,8 @@ DSMCFunc::DSMCFunc ( // TODO: Add a check that the first species is the electron species // (This should be done for impact ionization with MCC too) - - std::string secondary_species; - pp_collision_name.get("ionization_species", secondary_species); - m_species_names.push_back(secondary_species); + // 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)); } diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index f580df48bda..1171da2f1c7 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -137,8 +137,6 @@ public: { if (mask[i]) { - // TODO: by default: set all IDs to be invalid - 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 @@ -147,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 @@ -156,21 +155,29 @@ public: soa_products_data[1].m_rdata[PIdx::w][product2_index] = p_pair_reaction_weight[i]; if (ionization_flag) { - // The electron created from ionization is stored as the second particle in product 1 - // Copy its properties from species 2 (momentum, position, etc.) - // TODO - // Set the weight of the new particles to p_pair_reaction_weight[i] - - - // The ion created from ionization is stored as the first particle in product 3 - // Copy its properties from species 2 (momentum, position, etc.) const auto product3_index = products_np_data[2] + (p_offsets[i]*p_num_products_device[2] + 0); - // TODO + 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]; - // 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]; - soa_products_data[2].m_rdata[PIdx::w][product3_index] = 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 diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp index b28871bd375..070aa9caf58 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp @@ -26,7 +26,7 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name, if (m_collision_type == CollisionType::DSMC) { - if (ionization_flag) { + 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; From cfc9a9b6173d2a6c727e83b35985b30c43ce0018 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 12:48:53 -0700 Subject: [PATCH 11/14] Modify velocity of the electron --- .../BinaryCollision/DSMC/SplitAndScatterFunc.H | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index 1171da2f1c7..0e26c6da59c 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -237,7 +237,20 @@ public: amrex::Abort("Uneven mass charge-exchange not implemented yet."); } } else if (mask[i] == int(ScatteringProcessType::IONIZATION)) { - // Same as in doBackgroundIonization + // calculate kinetic energy of the incident electron + double E_coll; + double energy_cost = 0; // TODO: update this: get from scattering process + const 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 + ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine); + ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up, engine); } else { amrex::Abort("Unknown scattering process."); From d29b23f5d4d456342c1ad74e076fe0f11739817f Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 11 Oct 2024 12:51:12 -0700 Subject: [PATCH 12/14] Fix compilation errors --- .../BinaryCollision/DSMC/SplitAndScatterFunc.H | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index 0e26c6da59c..d46b7b56bac 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -240,7 +240,7 @@ public: // calculate kinetic energy of the incident electron double E_coll; double energy_cost = 0; // TODO: update this: get from scattering process - const ParticleReal u_coll2 = ux1*ux1 + uy1*uy1 + uz1*uz1; + 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); @@ -249,8 +249,13 @@ public: const amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m1; // isotropically scatter electrons - ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine); - ParticleUtils::RandomizeVelocity(e_ux, e_uy, e_uz, up, engine); + // - 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."); From 3b2001a7ec5909e8604d5fa53e355b6c8dd51bdd Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 12 Dec 2024 10:33:56 -0800 Subject: [PATCH 13/14] Update Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp --- .../Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp index 070aa9caf58..43e3351451e 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.cpp @@ -19,7 +19,7 @@ SplitAndScatterFunc::SplitAndScatterFunc (const std::string& collision_name, 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("excitation") != std::string::npos) { + if (scattering_process.find("ionization") != std::string::npos) { m_ionization_flag = true; } } From aa9ab75982715d47c2eeb50f7e0d2293254060f4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Mon, 16 Dec 2024 13:51:33 -0800 Subject: [PATCH 14/14] Update Source/Particles/Collision/BinaryCollision/BinaryCollision.H --- Source/Particles/Collision/BinaryCollision/BinaryCollision.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index a8cb3e4970a..67fe9443f27 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -112,7 +112,7 @@ public: 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) { + if (scattering_process.find("ionization") != std::string::npos) { ionization_flag = true; } }