WIP: electron impact ionization

Author: emarchet

docs/source/run/parameters.rst

@@ -500,20 +500,23 @@ When both are specified, the per-species value is used.
     the specific ionization energy of each state.
     Options are: ``electron``, ``positron``, ``H``, ``D``, ``T``, ``He``, ``Li``, ``Be``, ``B``, ….
 
-* ``<plasma name>.can_ionize`` (`bool`) optional (default `0`)
-    Whether this plasma can ionize. Can also be set to 1 by specifying ``<plasma name>.ionization_product``.
+* ``<plasma name>.can_field_ionize`` (`bool`) optional (default `0`)
+    Whether this plasma can be ionized by fields, as for example from a driving particle beam. 
 
-* ``<plasma name>.can_laser_ionize`` (`bool`) optional (default `<plasma name>.can_ionize`)
+* ``<plasma name>.can_laser_ionize`` (`bool`) optional (default `0`)
     Whether this plasma can be ionized by a laser.
 
+* ``<plasma name>.can_impact_ionize`` (`bool`) optional (default `0`)
+    Whether this plasma can be ionized by impact ionization between plasma particles.
+
 * ``<plasma name>.initial_ion_level`` (`int`) optional (default `-1`)
     The initial ionization state of the plasma. `0` for neutral gasses.
     If set, the plasma charge gets multiplied by this number. If the plasma species is not ionizable,
     the initial ionization level is set to 1.
 
 * ``<plasma name>.ionization_product`` (`string`) optional (default "")
     Name of the plasma species that contains the new electrons that are produced
-    when this plasma gets ionized. Only needed if this plasma is ionizable.
+    when this plasma gets field or laser ionized. Only needed if this plasma is ionizable.
 
 * ``<plasma name> or plasmas.neutralize_background`` (`bool`) optional (default `1`)
     Whether to add a neutralizing background of immobile particles of opposite charge.
@@ -1307,9 +1310,10 @@ WARNING: this module is in development.
 HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX,
 for collisions between plasma-plasma and beam-plasma.
 As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified.
+Please note that the impact ionization module is currently under development, and SI units are always required.
 
 * ``hipace.collisions`` (list of `strings`) optional
-    List of names of binary Coulomb collisions.
+    List of names of binary elastic Coulomb collisions.
     Each will represent collisions between 2 species.
 
 * ``<collision name>.species`` (two `strings`) optional
@@ -1321,6 +1325,32 @@ As collisions depend on the physical density, in normalized units `hipace.backgr
     Coulomb logarithm used for this collision.
     If not specified, the Coulomb logarithm is determined from the temperature in each cell.
 
+* ``hipace.impact_ionization`` (list of `strings`) optional
+    List of names of binary inelastic Coulomb collisions.
+    Each will represent collisions between 2 plasma species.
+    Three species must be specified: the projectile, the target (ion or neutral atom), and the product (electron).
+    Currently only impact ionization between electrons and another species, which might be ions or neutral atoms, is implemented.
+    Please note that the impact ionization module is currently under development, and SI units are always required.
+
+* ``<impact_name>.type`` (string) optional
+    Type of impact ionization. Currently, only `electron_impact` is implemented.
+
+* ``<impact name>.projectile`` (string) optional
+    The name of the projectile species for impact ionization.
+    The species must be defined according to `Plasma parameters`.
+
+* ``<impact name>.target`` (string) optional
+    The name of the target species for impact ionization.
+    The species must be defined according to `Plasma parameters`.
+
+* ``<impact name>.new_electron`` (string) optional
+    The name of the product species for impact ionization.
+    The species must be defined according to `Plasma parameters` and should be an electron species.
+    In electron impacts, initial and newly created electrons can be treated as on single species 
+    simply by using the same name for the projectile and the new electrons.
+    However, to treat the new electrons as a separate species, a different name can be used.
+    In that case, the new electrons won't be contributing to further ionizations.
+
 Radiation reaction
 ^^^^^^^^^^^^^^^^^^
 

src/Hipace.H

@@ -364,8 +364,9 @@ private:
     std::vector<std::string> m_collision_names;
     /** Vector of binary collisions */
     amrex::Vector< CoulombCollision > m_all_collisions;
-    amrex::Vector<std::string> m_collision_names2;
-    amrex::Vector<Collision> m_all_collisions2;
+    /** Impact ionization */
+    amrex::Vector<std::string> m_impact_names;
+    amrex::Vector<Collision> m_all_impacts;
 
     void InitDiagnostics (const int step, const amrex::Real time, const bool is_last_step);
     void FillFieldDiagnostics (const int current_N_level, int islice);

src/Hipace.cpp

@@ -261,10 +261,10 @@ Hipace::ReadParameters ()
             "be specified via 'hipace.background_density_SI'");
     }
 
-    queryWithParser(pph, "collisions2", m_collision_names2);
-    for (int i = 0; i != m_collision_names2.size(); ++i) {
-        m_all_collisions2.emplace_back(Collision());
-        m_all_collisions2.back().ReadParameters(m_multi_plasma.m_names, m_collision_names2[i]);
+    queryWithParser(pph, "impact_ionization", m_impact_names);
+    for (int i = 0; i != m_impact_names.size(); ++i) {
+        m_all_impacts.emplace_back(Collision());
+        m_all_impacts.back().ReadParameters(m_multi_plasma.m_names, m_impact_names[i]);
     }
 
     // external fields applied to the grid
@@ -850,8 +850,8 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st
     // collisions for plasmas and beams
     doCoulombCollision();
 
-    for (auto& collision : m_all_collisions2) {
-        collision.doCollision(0, m_slice_geom[0], m_multi_plasma);
+    for (auto& impact : m_all_impacts) {
+        impact.doCollision(0, m_slice_geom[0], m_multi_plasma);
     }
 
     // get minimum beam uz after push

src/particles/collisions/Collision.H

@@ -13,9 +13,9 @@
  */
 class Collision
 {
-    std::string m_inout_species1_name = "";
-    std::string m_inout_species2_name = "";
-    std::string m_out_species3_name = "";
+    std::string m_projectile_name = "";
+    std::string m_target_name = "";
+    std::string m_new_electron_name = "";
     std::string m_collision_type = "";
     std::string m_physical_element ="";
     bool m_has_collision_product = false;

src/particles/collisions/Collision.cpp

@@ -68,11 +68,11 @@ Collision::ReadParameters(
     // The user is asked to always input the target type as the second specie
     m_has_collision_product = m_collision_type == "electron_impact" || m_collision_type == "ion_impact";
 
-    getWithParser(pp, "species1", m_inout_species1_name);
-    getWithParser(pp, "species2", m_inout_species2_name);
+    getWithParser(pp, "projectile", m_projectile_name);
+    getWithParser(pp, "target", m_target_name);
 
     if (m_has_collision_product) {
-        getWithParser(pp, "species3", m_out_species3_name);
+        getWithParser(pp, "new_electron", m_new_electron_name);
     }
 
     if (m_collision_type == "electron_impact") {
@@ -90,7 +90,7 @@ Collision::ReadParameters(
 
     }
     // Plasma physical element name
-    amrex::ParmParse pp_s2(m_inout_species2_name);
+    amrex::ParmParse pp_s2(m_target_name);
     getWithParser(pp_s2, "element", m_physical_element);
 }
 
@@ -120,9 +120,9 @@ Collision::doCollisionA (
     const amrex::Real* p_binding_energies = m_binding_energies.data();
     const amrex::Real* p_ionization_energies = m_ionization_energies.data();
 
-    auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name);
-    auto& species2 = multi_plasma.GetPlasma(m_inout_species2_name);
-    auto& species3 = multi_plasma.GetPlasma(m_out_species3_name);
+    auto& species1 = multi_plasma.GetPlasma(m_projectile_name);
+    auto& species2 = multi_plasma.GetPlasma(m_target_name);
+    auto& species3 = multi_plasma.GetPlasma(m_new_electron_name);
 
     auto m1 = species1.GetMass();
     auto m2 = species2.GetMass();
@@ -134,7 +134,7 @@ Collision::doCollisionA (
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
         ion_atomic_number == 1 || 
         ion_atomic_number == 18, 
-        "The current implementation of electron-impact ionization only supports Hydrogen and Argon. Please check the input file and the physical element specified for species2."
+        "The current implementation of electron-impact ionization only supports Hydrogen and Argon. Please check the input file and the physical element specified for target."
     );
 
     const amrex::Real inv_dV = geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2);
@@ -374,8 +374,8 @@ Collision::doCollisionImp (
         G const& ionization_function)
 {
     // assume the two species are different for now
-    auto& species1 = multi_plasma.GetPlasma(m_inout_species1_name);
-    auto& species2 = multi_plasma.GetPlasma(m_inout_species2_name);
+    auto& species1 = multi_plasma.GetPlasma(m_projectile_name);
+    auto& species2 = multi_plasma.GetPlasma(m_target_name);
 
     for (PlasmaParticleIterator pti(species1); pti.isValid(); ++pti) {
         amrex::removeInvalidParticles(species1.ParticlesAt(0, pti));
@@ -442,34 +442,6 @@ Collision::doCollisionImp (
         auto p_flag1 = flag1.dataPtr();
         auto p_flag2 = flag2.dataPtr();
 
-        amrex::Gpu::DeviceVector<amrex::Real> cell_weight1(num_cells);
-        amrex::Gpu::DeviceVector<amrex::Real> cell_weight2(num_cells);
-        auto p_cell_weight1 = cell_weight1.dataPtr();
-        auto p_cell_weight2 = cell_weight2.dataPtr();
-
-        amrex::ParallelFor(2*num_cells,
-            [=] AMREX_GPU_DEVICE (int icell) {
-                if (icell < num_cells) {
-                    auto start = offset1[icell];
-                    auto stop = offset1[icell+1];
-                    amrex::Real loc_cell_weight1 = 0;
-                    for (int idx1 = start; idx1 < stop; ++idx1) {
-                        loc_cell_weight1 += ptd1.rdata(PlasmaIdx::w)[perm1[idx1]];
-                    }
-                    p_cell_weight1[icell] = loc_cell_weight1;
-                } else {
-                    icell -= num_cells;
-                    auto start = offset2[icell];
-                    auto stop = offset2[icell+1];
-                    amrex::Real loc_cell_weight2 = 0;
-                    for (int idx2 = start; idx2 < stop; ++idx2) {
-                        loc_cell_weight2 += ptd2.rdata(PlasmaIdx::w)[perm2[idx2]];
-                    }
-                    p_cell_weight2[icell] = loc_cell_weight2;
-                }
-            }
-        );
-
         // loop over independent pairs
         amrex::ParallelForRNG(total_ind_pairs,
             [=] AMREX_GPU_DEVICE (int ipair, amrex::RandomEngine const& engine){
@@ -485,9 +457,6 @@ Collision::doCollisionImp (
                 const int N1 = offset1_stop - offset1_start;
                 const int N2 = offset2_stop - offset2_start;
 
-                //const amrex::Real loc_cell_weight1 = p_cell_weight1[icell];
-                //const amrex::Real loc_cell_weight2 = p_cell_weight2[icell];
-
                 int idx1 = icoll;
                 int idx2 = icoll;
                 while (idx1 < N1 && idx2 < N2) {
@@ -554,12 +523,12 @@ Collision::doCollisionImp (
             amrex::Scan::Type::exclusive, amrex::Scan::retSum
         );
 
-        auto& species3 = multi_plasma.GetPlasma(m_out_species3_name);
+        auto& species3 = multi_plasma.GetPlasma(m_new_electron_name);
         auto& ptile3 = species3.ParticlesAt(0, pti);
         int old_size3 = ptile3.size();
         ptile3.resize(old_size3 + num_new_particles);
 
-        // get new ptd after resize in case species3 is the same as species1 or 2
+        // get new ptd after resize in case species3 is the same as species1 or species2
         ptd1 = ptile1.getParticleTileData();
         ptd2 = ptile2.getParticleTileData();
         auto ptd3 = ptile3.getParticleTileData();

src/particles/plasma/MultiPlasma.cpp

@@ -54,7 +54,7 @@ void
 MultiPlasma::InitIonization (amrex::Vector<amrex::Geometry> gm)
 {
     for (auto& plasma : m_all_plasmas) {
-        if(plasma.m_can_ionize) {
+        if(plasma.m_can_field_ionize || plasma.m_can_laser_ionize) {
             for (int i=0; i<m_names.size(); ++i) {
                 if(m_names[i] == plasma.m_product_name) {
                     plasma.m_product_pc = &m_all_plasmas[i];

src/particles/plasma/PlasmaParticleContainer.H

@@ -279,6 +279,7 @@ public:
     bool m_can_ionize = false; /**< whether this plasma can ionize from anything */
     bool m_can_field_ionize = false; /**< whether this plasma can ionize from the field */
     bool m_can_laser_ionize = false; /**< whether this plasma can ionize from a laser */
+    bool m_can_impact_ionize = false; /**< whether this plasma can ionize from particle impact */
     int m_init_ion_lev = -1; /**< initial Ion level of each particle */
     int m_max_ion_lev = 0; /**< maximum ionization level */
     std::string m_product_name = ""; /**< name of Ionization product plasma */

src/particles/plasma/PlasmaParticleContainer.cpp

@@ -65,14 +65,17 @@ PlasmaParticleContainer::ReadParameters ()
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_mass != 0, "The plasma particle mass must be specified");
 
     bool ion_lev_specified = queryWithParser(pp, "initial_ion_level", m_init_ion_lev);
-    // m_can_field_ionize = pp.contains("ionization_product");
 
-    // queryWithParser(pp, "can_ionize", m_can_field_ionize);
-    queryWithParser(pp, "can_ionize", m_can_ionize);
-    m_can_laser_ionize = false;
+    queryWithParser(pp, "can_field_ionize", m_can_field_ionize);
     queryWithParser(pp, "can_laser_ionize", m_can_laser_ionize);
+    queryWithParser(pp, "can_impact_ionize", m_can_impact_ionize);
 
-    // m_can_ionize = m_can_field_ionize || m_can_laser_ionize;    // !!!
+    m_can_ionize = m_can_field_ionize || m_can_laser_ionize || m_can_impact_ionize;
+
+    DeprecatedInput(m_name, "can_ionize", "can_field_ionize");
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_can_field_ionize == pp.contains("ionization_product") || m_can_laser_ionize == pp.contains("ionization_product"),
+        "When specifying ionization_product, can_*_ionize must be set to 1 via field or laser");
 
     if(m_can_ionize) {
         m_neutralize_background = false; // change default