Alternative injection condition

docs/source/run/parameters.rst

@@ -484,7 +484,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`)
@@ -639,6 +639,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.
+
 * ``<beam name> or beams.output_ratio`` (`int`) optional (default `1`)
     Set the fraction of beam particles that should be written to the openPMD output.
     For example, an output ratio of 100 will output every 100th beam particle.

src/Hipace.H

@@ -238,6 +238,7 @@ public:
     inline static bool m_do_beam_jz_minus_rho = false;
     /** Whether to deposit rho (plasma) for diagnostics */
     inline static bool m_deposit_rho = false;
+    inline static bool m_deposit_rho_beam = false;
     /** Whether to deposit rho for every individual plasma for diagnostics */
     inline static bool m_deposit_rho_individual = false;
     /** Whether to deposit temperature for every individual plasma for diagnostics */

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>
@@ -126,6 +127,10 @@ Hipace::ReadParameters ()
     queryWithParser(pph, "do_beam_jz_minus_rho", m_do_beam_jz_minus_rho);
     m_deposit_rho = m_diags.needsRho();
     queryWithParser(pph, "deposit_rho", m_deposit_rho);
+    queryWithParser(pph, "deposit_rho_beam", m_deposit_rho_beam);
+    if (m_deposit_rho_beam) {
+        m_deposit_rho = true;
+    }
     m_deposit_rho_individual = m_diags.needsRhoIndividual();
     queryWithParser(pph, "deposit_rho_individual", m_deposit_rho_individual);
     m_deposit_temp_individual = m_diags.needsTempIndividual();
@@ -545,7 +550,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);
 
@@ -805,6 +810,9 @@ Hipace::SolveOneSlice (int islice, int step)
         m_multi_plasma.AdvanceParticles(m_fields, m_3D_geom, false, lev, current_N_level);
     }
 
+    // injection
+    m_multi_plasma.DoLaserInjection(m_3D_geom, islice);
+
     if (m_depos_order_z == 2) {
         CalculateEzNext(current_N_level, step);
     }
@@ -813,7 +821,7 @@ Hipace::SolveOneSlice (int islice, int step)
     m_adaptive_time_step.GatherMinAccSlice(m_multi_beam, m_3D_geom[0], m_fields);
 
     // Push beam particles
-    m_multi_beam.AdvanceBeamParticlesSlice(m_fields, m_3D_geom, islice, current_N_level);
+    m_multi_beam.AdvanceBeamParticlesSlice(m_fields, m_3D_geom, islice, current_N_level, step);
 
     m_multi_beam.shiftSlippedParticles(islice, m_3D_geom[0]);
 

src/diagnostics/OpenPMDWriter.cpp

@@ -387,7 +387,7 @@ OpenPMDWriter::CopyBeams (MultiBeam& beams, const amrex::Vector< std::string > b
             }
 
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-                int(m_real_beam_data[ibeam].size()) == soa.NumRealComps(),
+                int(m_real_beam_data[ibeam].size()) == (soa.NumRealComps() - 1),
                 "List of real names in openPMD Writer class does not match the beam");
 
             for (std::size_t idx=0; idx<m_real_beam_data[ibeam].size(); idx++) {
@@ -397,9 +397,10 @@ OpenPMDWriter::CopyBeams (MultiBeam& beams, const amrex::Vector< std::string > b
                         std::max<uint64_t>(old_size+old_size/4, m_offset[ibeam] + np)
                     );
                 }
+                const int get_idx = idx >= BeamIdx::nsubcycles ? idx+1 : idx;
                 amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost,
-                    soa.GetRealData(idx).begin(),
-                    soa.GetRealData(idx).begin() + np,
+                    soa.GetRealData(get_idx).begin(),
+                    soa.GetRealData(get_idx).begin() + np,
                     m_real_beam_data[ibeam][idx].data() + m_offset[ibeam]);
             }
         }

src/particles/beam/BeamParticleContainer.H

@@ -26,17 +26,15 @@ struct BeamIdx
         x=0, y, z,  // position
         w,          // weight
         ux, uy, uz, // momentum
-        real_nattribs_in_buffer,
-        real_nattribs=real_nattribs_in_buffer
+        nsubcycles, // nsubcycles: by how many subcycles was this particle pushed already
+        real_nattribs,
+        real_nattribs_in_buffer = real_nattribs
     };
     enum {
-        // no extra components stored in MultiBuffer, besides 64bit idcpu
-        int_nattribs_in_buffer,
-        // nsubcycles: by how many subcycles was this particle pushed already
-        // nsubcycles is not stored or communicated in MultiBuffer
-        nsubcycles=int_nattribs_in_buffer,
         mr_level,
-        int_nattribs
+        int_nattribs,
+        // no extra components stored in MultiBuffer, besides 64bit idcpu
+        int_nattribs_in_buffer = mr_level
     };
 };
 
@@ -200,11 +198,11 @@ public:
 
     bool communicateIdCpuComponent () const { return true; }
     bool communicateRealComponent (int rcomp) const {
-        // communicate all compile-time and runtime real components
+        // otherwise communicate all compile-time and runtime real components
         return rcomp < numRealComponents();
     }
     bool communicateIntComponent (int icomp) const {
-        // don't communicate nsubcycles or mr_level
+        // don't communicate mr_level
         return icomp < BeamIdx::int_nattribs_in_buffer;
     }
 

src/particles/beam/BeamParticleContainer.cpp

@@ -399,7 +399,7 @@ BeamParticleContainer::initializeSlice (int slice, int which_slice) {
                     ptd.m_runtime_rdata[2][ip] = ptd_init.m_runtime_rdata[2][idx_src];
                 }
                 ptd.idcpu(ip) = ptd_init.idcpu(idx_src);
-                ptd.idata(BeamIdx::nsubcycles)[ip] = 0;
+                ptd.rdata(BeamIdx::nsubcycles)[ip] = 0;
                 ptd.idata(BeamIdx::mr_level)[ip] = 0;
             }
         );

src/particles/beam/BeamParticleContainerInit.cpp

@@ -107,7 +107,7 @@ namespace
         ptd.rdata(BeamIdx::uz)[ip] = uz;
         ptd.rdata(BeamIdx::w )[ip] = is_valid ? std::abs(weight) : amrex::Real{0};
 
-        ptd.idata(BeamIdx::nsubcycles)[ip] = 0;
+        ptd.rdata(BeamIdx::nsubcycles)[ip] = 0;
         ptd.idata(BeamIdx::mr_level)[ip] = 0;
 
         ptd.idcpu(ip) = pid + ip;

src/particles/beam/MultiBeam.H

@@ -44,10 +44,11 @@ public:
      * \param[in] gm Geometry object at level lev
      * \param[in] slice longitudinal slice
      * \param[in] current_N_level number of MR levels active on the current slice
+     * \param[in] step time step
      */
     void AdvanceBeamParticlesSlice (
         const Fields& fields, amrex::Vector<amrex::Geometry> const& gm, const int slice,
-        int const current_N_level);
+        int const current_N_level, int step);
 
     /** Compute reduced beam diagnostics of current slice, store in member variable.
      * \param[in] step time step of simulation

src/particles/beam/MultiBeam.cpp

@@ -71,11 +71,11 @@ MultiBeam::shiftSlippedParticles (const int slice, amrex::Geometry const& geom)
 void
 MultiBeam::AdvanceBeamParticlesSlice (
     const Fields& fields, amrex::Vector<amrex::Geometry> const& gm, const int slice,
-    int const current_N_level)
+    int const current_N_level, int step)
 {
     for (int i=0; i<m_nbeams; i++) {
         if (m_all_beams[i].m_do_push){
-            ::AdvanceBeamParticlesSlice(m_all_beams[i], fields, gm, slice, current_N_level);
+            ::AdvanceBeamParticlesSlice(m_all_beams[i], fields, gm, slice, current_N_level, step);
         }
     }
 

src/particles/deposition/BeamDepositCurrent.cpp

@@ -58,6 +58,7 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields,
     const int     jyb_cmp = do_beam_jx_jy_deposition  ? Comps[which_slice]["jy"    +beam_str] : -1;
     const int     jzb_cmp = do_beam_jz_deposition     ? Comps[which_slice]["jz"    +beam_str] : -1;
     const int rhomjzb_cmp = do_beam_rhomjz_deposition ? Comps[which_slice]["rhomjz"+beam_str] : -1;
+    const int     rho_cmp = (which_slice == WhichSlice::This) && Hipace::m_deposit_rho_beam ? Comps[which_slice]["rho"] : -1;
 
     // Offset for converting positions to indexes
     amrex::Real const x_pos_offset = GetPosOffset(0, gm[lev], isl_fab.box());
@@ -83,6 +84,7 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields,
 
     const amrex::Real clight = phys_const.c;
     const amrex::Real q = beam.m_charge;
+    const amrex::Real n_subcycles = static_cast<amrex::Real>(beam.m_n_subcycles);
 
     amrex::AnyCTO(
         // use compile-time options
@@ -100,7 +102,7 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields,
                 isl_fab.array(), isl_fab.box(),
                 beam.getBeamSlice(which_beam_slice).getParticleTileData(),
                 amrex::GpuArray<int, 0>{},
-                amrex::GpuArray<int, 4>{jxb_cmp, jyb_cmp, jzb_cmp, rhomjzb_cmp});
+                amrex::GpuArray<int, 5>{jxb_cmp, jyb_cmp, jzb_cmp, rhomjzb_cmp, rho_cmp});
         },
         // is_valid
         // return whether the particle is valid and should deposit
@@ -111,7 +113,8 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields,
             return ptd.id(ip).is_valid() &&
                 (only_highest ?
                     (ptd.idata(BeamIdx::mr_level)[ip] == lev) :
-                    (ptd.idata(BeamIdx::mr_level)[ip] >= lev));
+                    (ptd.idata(BeamIdx::mr_level)[ip] >= lev)) &&
+                ptd.rdata(BeamIdx::nsubcycles)[ip] <= n_subcycles;
         },
         // get_cell
         // return the lowest cell index that the particle deposits into
@@ -187,6 +190,11 @@ DepositCurrentSlice (BeamParticleContainer& beam, Fields& fields,
                             arr.ptr(i_cell+ix, j_cell+iy, depos_idx[3]),
                             sx_cell[ix]*sy_cell[iy]*wqrhomjz);
                     }
+                    if (depos_idx[4] != -1) { // deposit_rho
+                        amrex::Gpu::Atomic::Add(
+                            arr.ptr(i_cell+ix, j_cell+iy, depos_idx[4]),
+                            sx_cell[ix]*sy_cell[iy]*wq);
+                    }
                 }
             }
         });

src/particles/deposition/PlasmaDepositCurrent.cpp

@@ -36,7 +36,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
     " (WhichSlice::Next), for the ion charge deposition (WhichSLice::RhomJzIons)"
     " or for the Salame slice (WhichSlice::Salame)");
 
-    const amrex::Real max_qsa_weighting_factor = plasma.m_max_qsa_weighting_factor;
+    // const amrex::Real max_qsa_weighting_factor = plasma.m_max_qsa_weighting_factor;
     const amrex::Real charge = (which_slice == WhichSlice::RhomJzIons) ? -plasma.m_charge : plasma.m_charge;
     const amrex::Real mass = plasma.m_mass;
     // only deposit rho individual on WhichSlice::This
@@ -201,16 +201,17 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields,
                     + 1._rt
                 );
 
-                if (gamma_psi < 0.0_rt || gamma_psi > max_qsa_weighting_factor || psi_inv < 0.0_rt)
-                {
-                    // This particle violates the QSA, discard it and do not deposit its current
-                    if (p_n_qsa_violation) {
-                        amrex::Gpu::Atomic::Add(p_n_qsa_violation, 1);
-                    }
-                    ptd.rdata(PlasmaIdx::w)[ip] = 0.0_rt;
-                    ptd.id(ip).make_invalid();
-                    return;
-                }
+                // if (gamma_psi < 0.0_rt || (gamma_psi > max_qsa_weighting_factor && ptd.id(ip)!=3) || psi_inv < 0.0_rt)
+                // {
+                //     // AMREX_DEVICE_PRINTF("gamma_psi of plasma particle: %f psi_inv: %f\n", gamma_psi, psi_inv);
+                //     // This particle violates the QSA, discard it and do not deposit its current
+                //     if (p_n_qsa_violation) {
+                //         amrex::Gpu::Atomic::Add(p_n_qsa_violation, 1);
+                //     }
+                //     ptd.rdata(PlasmaIdx::w)[ip] = 0.0_rt;
+                //     ptd.id(ip).make_invalid();
+                //     return;
+                // }
 
                 for (int iy=0; iy <= depos_order; ++iy) {
                     for (int ix=0; ix <= depos_order; ++ix) {

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
 {
@@ -26,10 +27,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
@@ -115,6 +118,13 @@ 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] gm Geometry object for the whole domain
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void DoLaserInjection (amrex::Vector<amrex::Geometry> const& gm, const int islice);
+
     /** \brief whether any plasma species uses a neutralizing background, e.g. no ion motion */
     [[nodiscard]] 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/TemperatureDeposition.H"
 #include "particles/deposition/ExplicitDeposition.H"
@@ -41,7 +42,8 @@ MultiPlasma::ReadParameters ()
 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
@@ -60,7 +62,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
@@ -151,6 +164,15 @@ MultiPlasma::DoLaserIonization (
     }
 }
 
+void
+MultiPlasma::DoLaserInjection (
+    amrex::Vector<amrex::Geometry> const& gm, const int islice)
+{
+    for (auto& plasma : m_all_plasmas) {
+        plasma.PlasmaToBeam(gm, islice);
+    }
+}
+
 bool
 MultiPlasma::AnySpeciesNeutralizeBackground () const
 {

src/particles/plasma/PlasmaParticleContainer.H

@@ -16,6 +16,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) */
@@ -31,6 +33,7 @@ struct PlasmaIdx
         ux_half_step,       // momentum half a step behind the current slice for leapfrog pusher
         uy_half_step,       // at the same step for AB5 pusher
         psi_half_step,      // never effected by temp slice
+        time_integral,
 #ifdef HIPACE_USE_AB5_PUSH
         Fx1, Fx2, Fx3, Fx4, Fx5,             // AB5 force terms
         Fy1, Fy2, Fy3, Fy4, Fy5,             //
@@ -111,6 +114,13 @@ 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] gm Geometry object for the whole domain
+     * \param[in] islice Slice index corresponding to the current simulation slice
+     */
+    void PlasmaToBeam (amrex::Vector<amrex::Geometry> const& gm, const int islice);
+
     /** Reorder particles to speed-up current deposition
      * \param[in] islice zeta slice index
      */
@@ -219,6 +229,10 @@ public:
     amrex::Gpu::Buffer<amrex::Real> m_laser_dp_prefactor;
     /** to calculate laser ionization momentum width at the second order for linear polarization */
     amrex::Gpu::Buffer<amrex::Real> m_laser_dp_second_prefactor;
+    // laser injection:
+    bool m_can_laser_injection = false; /**< whether laser injection is enabled */
+    std::string m_product_beam_name = ""; /**< name of Injection product beam */
+    BeamParticleContainer* m_product_beam_pc = nullptr; /**< Injection product beam */
 
     // plasma sorting/reordering: