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

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"},

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<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();
 

Source/Particles/Collision/BinaryCollision/DSMC/CollisionFilterFunc.H

@@ -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);

Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H

@@ -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;

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<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) {
@@ -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;
 }

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

@@ -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);
@@ -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<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];
@@ -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);
                 }
                 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.

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<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:
+            // 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