[WIP] Laser Injection

docs/source/run/parameters.rst

@@ -595,6 +595,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 unabled. 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>
@@ -698,6 +699,9 @@ Hipace::SolveOneSlice (int islice, int step)
     // plasma laser ionization
     m_multi_plasma.DoLaserIonization(islice, m_multi_laser.GetLaserGeom(), m_multi_laser);
 
+    // injection
+    m_multi_plasma.DoLaserInjection(m_multi_beam, getDiagBeamNames(), m_multi_laser.GetLaserGeom(), m_multi_laser, 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

@@ -258,7 +258,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/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
 {
@@ -109,6 +110,17 @@ 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] beams MultiBeam object containing all beams
+     * \param[in] beamnames Vector of strings representing the names of the beams
+     * \param[in] laser_geom Geometry of the laser used for ionization
+     * \param[in] laser MultiLaser object containing the laser details
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void DoLaserInjection (MultiBeam& beams, const amrex::Vector< std::string > beamnames, const amrex::Geometry& laser_geom,
+                 const MultiLaser& laser, 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"
@@ -57,7 +58,19 @@ MultiPlasma::InitData (amrex::Vector<amrex::BoxArray> slice_ba,
             plasma.InitIonizationModule(gm[0],
                 Hipace::m_background_density_SI); // geometry only for dz
         }
+
+        MultiBeam multi_beam;
+        if(plasma.m_can_laser_injection) {
+            for (int i=0; i<m_names.size(); ++i) {
+                if(m_names[i] == plasma.m_product_beam_name) {
+                    plasma.m_product_beam_pc = &multi_beam.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 (
+    MultiBeam& beams, const amrex::Vector< std::string > beamnames, const amrex::Geometry& laser_geom,
+                 const MultiLaser& laser, const int islice)
+{
+    for (auto& plasma : m_all_plasmas) {
+        plasma.PlasmaToBeam(beams, beamnames, laser_geom,laser, 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,18 @@ public:
                           const MultiLaser& laser,
                           const amrex::Real background_density_SI);
 
+    /** Do the transfer of particles from the plasma container to the beam container
+     *
+     * \param[in] beams MultiBeam object containing all beams
+     * \param[in] beamnames Vector of strings representing the names of the beams
+     * \param[in] laser_geom Geometry of the laser used for ionization
+     * \param[in] laser MultiLaser object containing the laser details
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void PlasmaToBeam (MultiBeam& beams, const amrex::Vector< std::string > beamnames, const amrex::Geometry& laser_geom,
+                 const MultiLaser& laser, const int islice
+                 );
+
     /** Reorder particles to speed-up current deposition
      * \param[in] islice zeta slice index
      */
@@ -215,8 +229,11 @@ 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 unabled */
+    std::string m_product_beam_name = ""; /**< name of Injection product beam */
+    BeamParticleContainer* m_product_beam_pc = nullptr; /**< Injection product beam */
 
     // plasma sorting/reordering:
 

src/particles/plasma/PlasmaParticleContainer.cpp

@@ -20,6 +20,9 @@
 #include "particles/pusher/BeamParticleAdvance.H"
 #include "particles/particles_utils/FieldGather.H"
 #include "particles/pusher/GetAndSetPosition.H"
+#include "particles/beam/BeamParticleContainer.H"
+#include "particles/beam/MultiBeam.H"
+
 #include <cmath>
 #include <fstream>
 #include <sstream>
@@ -78,6 +81,8 @@ PlasmaParticleContainer::ReadParameters ()
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_init_ion_lev >= 0,
             "The initial ion level must be specified");
     }
+    m_can_laser_injection = pp.contains("injection_product");
+    queryWithParser(pp, "can_laser_injection", m_can_laser_injection);
     queryWithParserAlt(pp, "neutralize_background", m_neutralize_background, pp_alt);
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_can_ionize || !m_neutralize_background,
         "Cannot use neutralize_background when ionization is turned on");
@@ -91,6 +96,7 @@ PlasmaParticleContainer::ReadParameters ()
         m_charge *= m_init_ion_lev;
     }
     queryWithParser(pp, "ionization_product", m_product_name);
+    queryWithParser(pp, "injection_product", m_product_beam_name);
 
     std::string density_func_str = "0.";
     DeprecatedInput(m_name, "density", "density(x,y,z)");
@@ -426,7 +432,6 @@ IonizationModule (const int lev,
                 const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); // ensures thread-safe access when incrementing `p_ip_elec`
                 const long pidx = pid + old_size;
 
-                // Copy ion data to new electron
                 // Set the ionized electron ID to 2 (valid/invalid) for the ionized electrons
                 amrex::ParticleIDWrapper{idcpu_elec[pidx]} = 2;
                 amrex::ParticleCPUWrapper{idcpu_elec[pidx]} = lev; // current level
@@ -730,6 +735,81 @@ LaserIonization (const int islice,
     }
 }
 
+void
+PlasmaParticleContainer::
+PlasmaToBeam (MultiBeam& beams, const amrex::Vector< std::string > beamnames, const amrex::Geometry& laser_geom,
+                 const MultiLaser& laser, const int islice)
+{
+    if (!m_can_laser_ionize || !laser.UseLaser(islice) || !m_can_laser_injection) return;
+    HIPACE_PROFILE("PlasmaParticleContainer::PlasmaToBeam()");
+
+    amrex::Gpu::DeviceScalar<uint32_t> num_new_beam_part(0);
+    uint32_t* AMREX_RESTRICT p_num_new_beam_part = num_new_beam_part.dataPtr();
+
+    // Loop over plasma particle boxes
+    for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti)
+    {
+        // Loading the data from the plasma containers
+        const auto ptd_plasma = pti.GetParticleTile().getParticleTileData();
+
+        amrex::Long const num_particles = pti.numParticles();
+
+        amrex::ReduceOps<amrex::ReduceOpSum> reduce_op;
+        amrex::ReduceData<uint64_t> reduce_data(reduce_op);
+        using ReduceTuple = typename decltype(reduce_data)::Type;
+
+        // This kernel calculates the number of ionized electrons in the plasma container and make them invalid
+        reduce_op.eval(
+            num_particles, reduce_data,
+            [=] AMREX_GPU_DEVICE (int ip) -> ReduceTuple
+            {
+                if (ptd_plasma.id(ip) != 2) // whether the plasma particle is from ionization
+                {
+                    ptd_plasma.id(ip).make_invalid(); // make the particle invalid in the plasma container
+                    return {1};
+                } else {
+                    return {0};
+                }
+        });
+
+        auto [sum_new_beam_part] = reduce_data.value();
+
+        amrex::Gpu::Atomic::Add(p_num_new_beam_part, static_cast<uint32_t>(sum_new_beam_part));
+
+        auto& beam_elec = m_product_beam_pc;
+
+        // Resize the beam container
+        auto& beam_soa = beam_elec->getBeamSlice(WhichBeamSlice::This).GetStructOfArrays();
+        auto old_size = beam_soa.size();
+        auto new_size = old_size + num_new_beam_part.dataValue();
+        beam_soa.resize(new_size);
+
+        auto ptd_beam = beam_elec->getBeamSlice(WhichBeamSlice::This).getParticleTileData();
+
+        // This kernel does the transfer of the ionized electrons from the plasma container
+        // to the beam container
+        amrex::Gpu::DeviceScalar<uint32_t> ip_elec(0);
+        uint32_t * AMREX_RESTRICT p_ip_elec = ip_elec.dataPtr();
+        amrex::ParallelFor(num_particles,
+            [=] AMREX_GPU_DEVICE (int ip) {
+                if (ptd_plasma.id(ip) != 2){
+                    const long pid_beam = amrex::Gpu::Atomic::Add(p_ip_elec, 1u);
+                    const long pidx_beam = pid_beam + old_size;
+                    ptd_beam.rdata(BeamIdx::x)[pidx_beam] = ptd_plasma.pos(0, ip);
+                    ptd_beam.rdata(BeamIdx::y)[pidx_beam] = ptd_plasma.pos(1, ip);
+                    ptd_beam.rdata(BeamIdx::z)[pidx_beam] = ptd_plasma.pos(2, ip);
+                    ptd_beam.rdata(BeamIdx::ux)[pidx_beam] = ptd_plasma.rdata(PlasmaIdx::ux)[ip];
+                    ptd_beam.rdata(BeamIdx::uy)[pidx_beam] = ptd_plasma.rdata(PlasmaIdx::uy)[ip];
+                    amrex::Real ux = ptd_plasma.rdata(PlasmaIdx::ux)[ip];
+                    amrex::Real uy = ptd_plasma.rdata(PlasmaIdx::uy)[ip];
+                    amrex::Real psi = ptd_plasma.rdata(PlasmaIdx::psi)[ip];
+                    ptd_beam.rdata(BeamIdx::uz)[pidx_beam] = (1+ux*ux+uy*uy-psi*psi)/(2.*psi);
+                    ptd_beam.rdata(BeamIdx::w)[pidx_beam] = ptd_plasma.rdata(PlasmaIdx::w)[ip];
+                }
+            });
+    }
+}
+
 void
 PlasmaParticleContainer::InSituComputeDiags (int islice)
 {