[WIP] Complete version for ionization injection

docs/source/run/parameters.rst

@@ -453,7 +453,7 @@ When both are specified, the per-species value is used.
 * ``<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_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>.initial_ion_level`` (`int`) optional (default `-1`)
@@ -605,6 +605,13 @@ which are valid only for certain beam types, are introduced further below under
     The ideal index type is different for beam push and beam deposition so some experimentation
     may be required to find the overall fastest setting for a specific simulation.
 
+* ``<plasma name>.can_laser_injection`` (`bool`) optional (default `0`)
+    Whether the laser injection is enabled. Can also be set to 1 by specifying ``<plasma name>.injection_product``.
+
+* ``<plasma name>.injection_product`` (`string`) optional (default "")
+    Name of the beam species that contains the new electrons that are produced
+    when this plasma gets ionized. Only needed if this plasma is ionizable and the laser injection is unabled.
+
 Option: ``fixed_weight_pdf``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

src/Hipace.cpp

@@ -18,6 +18,7 @@
 #include "particles/pusher/GetAndSetPosition.H"
 #include "mg_solver/HpMultiGrid.H"
 #include "fields/fft_poisson_solver/fft/AnyFFT.H"
+#include "particles/beam/MultiBeam.H"
 
 #include <AMReX_ParmParse.H>
 #include <AMReX_IntVect.H>
@@ -458,7 +459,7 @@ Hipace::Evolve ()
         }
 
         // Only reset plasma after receiving time step, to use proper density
-        m_multi_plasma.InitData(m_slice_ba, m_slice_dm, m_slice_geom, m_3D_geom);
+        m_multi_plasma.InitData(m_slice_ba, m_slice_dm, m_slice_geom, m_3D_geom, m_multi_beam);
 
         m_multi_laser.SetInitialChi(m_multi_plasma);
 
@@ -714,6 +715,11 @@ Hipace::SolveOneSlice (int islice, int step)
     // plasma laser ionization
     m_multi_plasma.DoLaserIonization(islice, m_multi_laser.GetLaserGeom(), m_multi_laser);
 
+    // injection
+    for (int lev=0; lev<current_N_level; ++lev) {
+        m_multi_plasma.DoLaserInjection(lev, m_fields, m_multi_laser, m_3D_geom, islice);
+    }
+
     // Push plasma particles
     for (int lev=0; lev<current_N_level; ++lev) {
         m_multi_plasma.AdvanceParticles(m_fields, m_3D_geom, false, lev, current_N_level);

src/particles/beam/BeamParticleContainer.H

@@ -192,6 +192,10 @@ public:
         return m_total_num_particles;
     }
 
+    void updateSize (int n) {
+        m_total_num_particles = n;
+    }
+
     int numRealComponents () const {
         return BeamIdx::real_nattribs + (m_do_spin_tracking ? 3 : 0);
     }
@@ -259,7 +263,7 @@ private:
     amrex::IntVect m_reorder_idx_type = {0, 0, 0};
 
     // fixed_ppc:
-
+private:
     amrex::IntVect m_ppc {1, 1, 1}; /**< Number of particles per cell in each direction */
     amrex::RealVect m_position_mean {0., 0., 0.}; /**< mean position of the beam */
     amrex::Real m_min_density {0.}; /**< minimum density at which beam particles are generated */

src/particles/deposition/PlasmaDepositCurrent.cpp

@@ -196,7 +196,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
                     + 1._rt
                 );
 
-                if (gamma_psi < 0.0_rt || gamma_psi > max_qsa_weighting_factor || psi_inv < 0.0_rt)
+                if (gamma_psi < 0.0_rt || (gamma_psi > max_qsa_weighting_factor && ptd.id(ip)!=3) || psi_inv < 0.0_rt)
                 {
                     // This particle violates the QSA, discard it and do not deposit its current
                     amrex::Gpu::Atomic::Add(p_n_qsa_violation, 1);

src/particles/plasma/MultiPlasma.H

@@ -11,6 +11,7 @@
 #include "PlasmaParticleContainer.H"
 #include "fields/Fields.H"
 #include "particles/collisions/CoulombCollision.H"
+#include "particles/beam/MultiBeam.H"
 
 class MultiPlasma
 {
@@ -30,10 +31,12 @@ public:
      * \param[in] slice_dm DistributionMapping of the transverse slice domain
      * \param[in] slice_gm slice geometry
      * \param[in] gm Geometry of the simulation, to get the cell size
+     * \param[in] beams MultiBeam object containing all beams
      */
     void InitData (amrex::Vector<amrex::BoxArray> slice_ba,
                    amrex::Vector<amrex::DistributionMapping> slice_dm,
-                   amrex::Vector<amrex::Geometry> slice_gm, amrex::Vector<amrex::Geometry> gm);
+                   amrex::Vector<amrex::Geometry> slice_gm, amrex::Vector<amrex::Geometry> gm,
+                   MultiBeam& beams);
 
 
     /** Loop over plasma species and depose their currents into the current 2D slice in fields
@@ -109,6 +112,16 @@ public:
     void DoLaserIonization (const int islice, const amrex::Geometry& laser_geom,
                             const MultiLaser& laser);
 
+    /** Do the transfer of ionized plasma particles to beam particles
+     *
+     * \param[in] lev MR level
+     * \param[in] fields the general field class
+     * \param[in] laser MultiLaser object containing the laser details
+     * \param[in] gm Geometry object for the whole domain
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void DoLaserInjection (const int lev, const Fields& fields, const MultiLaser& laser, amrex::Vector<amrex::Geometry> const& gm, const int islice);
+
     /** \brief whether any plasma species uses a neutralizing background, e.g. no ion motion */
     bool AnySpeciesNeutralizeBackground () const;
 

src/particles/plasma/MultiPlasma.cpp

@@ -6,6 +6,7 @@
  * License: BSD-3-Clause-LBNL
  */
 #include "MultiPlasma.H"
+#include "particles/beam/MultiBeam.H"
 #include "particles/deposition/PlasmaDepositCurrent.H"
 #include "particles/deposition/ExplicitDeposition.H"
 #include "particles/pusher/PlasmaParticleAdvance.H"
@@ -38,7 +39,8 @@ MultiPlasma::MultiPlasma ()
 void
 MultiPlasma::InitData (amrex::Vector<amrex::BoxArray> slice_ba,
                        amrex::Vector<amrex::DistributionMapping> slice_dm,
-                       amrex::Vector<amrex::Geometry> slice_gm, amrex::Vector<amrex::Geometry> gm)
+                       amrex::Vector<amrex::Geometry> slice_gm, amrex::Vector<amrex::Geometry> gm,
+                       MultiBeam& beams)
 {
     for (auto& plasma : m_all_plasmas) {
         // make it think there is only level 0
@@ -57,7 +59,18 @@ MultiPlasma::InitData (amrex::Vector<amrex::BoxArray> slice_ba,
             plasma.InitIonizationModule(gm[0],
                 Hipace::m_background_density_SI); // geometry only for dz
         }
+
+        if(plasma.m_can_laser_injection) {
+            for (int i=0; i<beams.m_names.size(); ++i) {
+                if(beams.m_names[i] == plasma.m_product_beam_name) {
+                    plasma.m_product_beam_pc = &beams.m_all_beams[i];
+                }
+            }
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plasma.m_product_beam_pc != nullptr,
+                "Must specify a valid product beam for laser injection using ionization_product");
+        }
     }
+
 }
 
 amrex::Real
@@ -138,6 +151,16 @@ MultiPlasma::DoLaserIonization (
     }
 }
 
+void
+MultiPlasma::DoLaserInjection (
+    const int lev, const Fields& fields, const MultiLaser& laser, amrex::Vector<amrex::Geometry> const& gm, const int islice)
+{
+    for (auto& plasma : m_all_plasmas) {
+        plasma.InjectionCondition(lev, fields, laser, islice);
+        plasma.PlasmaToBeam(laser, gm, islice);
+    }
+}
+
 bool
 MultiPlasma::AnySpeciesNeutralizeBackground () const
 {

src/particles/plasma/PlasmaParticleContainer.H

@@ -15,6 +15,8 @@
 #include <AMReX_AmrParticles.H>
 #include <AMReX_Particles.H>
 #include <AMReX_AmrCore.H>
+#include "particles/beam/MultiBeam.H"
+#include "particles/beam/BeamParticleContainer.H"
 #include <map>
 
 /** \brief Map names and indices for plasma particles attributes (SoA data) */
@@ -113,6 +115,23 @@ public:
                           const MultiLaser& laser,
                           const amrex::Real background_density_SI);
 
+    /** Mark the plasma particles that has been injected
+     *
+     * \param[in] lev MR level
+     * \param[in] fields the general field class
+     * \param[in] laser MultiLaser object containing the laser details
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void InjectionCondition (const int lev, const Fields& fields, const MultiLaser& laser, const int islice);
+
+    /** Do the transfer of particles from the plasma container to the beam container
+     *
+     * \param[in] laser MultiLaser object containing the laser details
+     * \param[in] gm Geometry object for the whole domain
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void PlasmaToBeam (const MultiLaser& laser, amrex::Vector<amrex::Geometry> const& gm, const int islice);
+
     /** Reorder particles to speed-up current deposition
      * \param[in] islice zeta slice index
      */
@@ -215,8 +234,13 @@ public:
     amrex::Gpu::DeviceVector<amrex::Real> m_laser_adk_prefactor;
     /** to calculate laser ionization momentum width for linear polarization */
     amrex::Gpu::DeviceVector<amrex::Real> m_laser_dp_prefactor;
-    /** to calculate laser ionization momentum width at the second order for linear polarization */
     amrex::Gpu::DeviceVector<amrex::Real> m_laser_dp_second_prefactor;
+    // laser injection:
+    bool m_can_laser_injection = false; /**< whether laser injection is enabled */
+    amrex::Real m_uz_threshold = 1.;
+    amrex::Real m_injection_weight_factor = 1.;
+    std::string m_product_beam_name = ""; /**< name of Injection product beam */
+    BeamParticleContainer* m_product_beam_pc = nullptr; /**< Injection product beam */
 
     // plasma sorting/reordering: