Add beam-plasma collisions (#983)

* prepare functions

* add basic framework

* add transverse beam particle binning per slice

* add beam collisions

* fix doc

* fix CI

* fix automatic temperature calculation for beams

* add benchmark CI test for collisions with beam

* reset benchmark for new CI test

* cleaning and doc

* fix doc

* Apply suggestions from code review

* Update docs/source/run/parameters.rst

* Update docs/source/run/parameters.rst

---------

Co-authored-by: Tools <hipace.tools@desy.de>

CMakeLists.txt

@@ -345,6 +345,12 @@ if(BUILD_TESTING)
                     WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
             )
 
+            add_test(NAME collisions_beam.SI.1Rank
+                    COMMAND bash ${HiPACE_SOURCE_DIR}/tests/collisions_beam.SI.1Rank.sh
+                            $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
+                    WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+            )
+
             add_test(NAME ion_motion.SI.1Rank
                     COMMAND bash ${HiPACE_SOURCE_DIR}/tests/ion_motion.SI.1Rank.sh
                             $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}

docs/source/run/parameters.rst

@@ -810,15 +810,18 @@ Binary collisions
 
 WARNING: this module is in development.
 
-HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX, between plasma species.
+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.
 
-* ``plasmas.collisions`` (list of `strings`) optional
-    List of names of types binary Coulomb collisions.
-    Each will represent collisions between 2 plasma species (potentially the same).
+* ``hipace.collisions`` (list of `strings`) optional
+    List of names of binary Coulomb collisions.
+    Each will represent collisions between 2 species.
 
 * ``<collision name>.species`` (two `strings`) optional
-    The name of the two plasma species for which collisions should be included.
+    The name of the two species for which collisions should be included.
+    This can either be plasma-plasma or beam-plasma collisions. For plasma-plasma collisions, the species can be the same to model collisions within a species.
+    The names must be in `plasmas.names` or `beams.names` (for beam-plasma collisions).
 
 * ``<collision name>.CoulombLog`` (`float`) optional (default `-1.`)
     Coulomb logarithm used for this collision.

src/Hipace.H

@@ -162,6 +162,12 @@ public:
     /** \brief reset all laser slices to 0 */
     void ResetLaser ();
 
+    /**
+     * \brief does Coulomb collisions between plasmas and beams
+     * \param[in] islice_local slice on which collisions are calculated, needed for beam binning
+     */
+    void doCoulombCollision (const int islice_local);
+
     /**
        \brief Returns true on the head rank, otherwise false.
      */
@@ -310,6 +316,8 @@ public:
     inline static amrex::Real m_background_density_SI = 0.;
     /** Time step for the beam evolution */
     amrex::Real m_dt = 0.0;
+    /** Number of binary collisions */
+    inline static int m_ncollisions = 0;
     /** Adaptive time step instance */
     AdaptiveTimeStep m_adaptive_time_step;
     /** Laser instance (soon to be multi laser container) */
@@ -368,6 +376,11 @@ private:
     /** Diagnostics */
     Diagnostic m_diags;
 
+    /** User-input names of the binary collisions to be used */
+    std::vector<std::string> m_collision_names;
+    /** Vector of binary collisions */
+    amrex::Vector< CoulombCollision > m_all_collisions;
+
     /** \brief resizes the diagnostic fab to the correct box in a loop over boxes
      *
      * \param[in] it index of box to be resized to

src/Hipace.cpp

@@ -170,6 +170,18 @@ Hipace::Hipace () :
     MakeGeometry();
 
     m_use_laser = m_multi_laser.m_use_laser;
+
+    queryWithParser(pph, "collisions", m_collision_names);
+    /** Initialize the collision objects */
+    m_ncollisions = m_collision_names.size();
+     for (int i = 0; i < m_ncollisions; ++i) {
+         m_all_collisions.emplace_back(CoulombCollision(m_multi_plasma.m_names, m_multi_beam.m_names, m_collision_names[i]));
+     }
+     if (m_normalized_units && m_ncollisions > 0) {
+         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_background_density_SI!=0,
+             "For collisions with normalized units, a background plasma density must "
+             "be specified via 'hipace.background_density_SI'");
+     }
 }
 
 Hipace::~Hipace ()
@@ -618,8 +630,8 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step)
                                            m_external_E_uniform, m_external_B_uniform,
                                            m_external_E_slope, m_external_B_slope);
 
-    // collisions for all particles calculated on level 0
-    m_multi_plasma.doCoulombCollision(0, m_slice_geom[0].Domain(), m_slice_geom[0]);
+    // collisions for plasmas and beams
+    doCoulombCollision(islice_local);
 
     // Advance laser slice by 1 step and store result to 3D array
     // no MR for laser
@@ -1350,6 +1362,39 @@ Hipace::NotifyFinish (const int it, bool only_ghost, bool only_time)
 #endif
 }
 
+void
+Hipace::doCoulombCollision (const int islice_local)
+{
+
+    // collisions for all particles calculated on level 0
+    const int lev = 0;
+
+    for (int i = 0; i < m_ncollisions; ++i)
+    {
+        if (m_all_collisions[i].m_nbeams == 1) {
+            // do beam-plasma collisions
+            auto& species1 = m_multi_beam.m_all_beams[ m_all_collisions[i].m_species1_index ];
+            auto& species2 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species2_index ];
+
+            // TODO: enable tiling
+
+            CoulombCollision::doBeamPlasmaCoulombCollision(
+                lev, m_slice_geom[0].Domain(), m_slice_geom[0], islice_local, species1, species2,
+                m_all_collisions[i].m_CoulombLog, m_background_density_SI);
+        } else {
+            // do plasma-plasma collisions
+            auto& species1 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species1_index ];
+            auto& species2 = m_multi_plasma.m_all_plasmas[ m_all_collisions[i].m_species2_index ];
+
+            // TODO: enable tiling
+
+            CoulombCollision::doPlasmaPlasmaCoulombCollision(
+                lev, m_slice_geom[0].Domain(), m_slice_geom[0], species1, species2, m_all_collisions[i].m_isSameSpecies,
+                m_all_collisions[i].m_CoulombLog, m_background_density_SI);
+        }
+    }
+}
+
 void
 Hipace::ResizeFDiagFAB (const int it, const int step)
 {

src/particles/beam/BeamParticleContainer.H

@@ -157,6 +157,10 @@ public:
     std::string get_name () const {return m_name;}
     amrex::Real m_charge; /**< charge of each particle of this species */
     amrex::Real m_mass; /**< mass of each particle of this species */
+    /** Returns elementary charge q_e (or -q_e for electrons). */
+    amrex::Real GetCharge () const {return m_charge;}
+    /** Returns mass of physical species */
+    amrex::Real GetMass () const {return m_mass;}
     bool m_do_z_push {true}; /**< Pushing beam particles in z direction */
     int m_n_subcycles {10}; /**< Number of sub-cycles in the beam pusher */
     bool m_do_radiation_reaction {false}; /**< whether to calculate radiation losses */

src/particles/beam/MultiBeam.H

@@ -193,11 +193,11 @@ public:
     /** \brief returns if any beam uses the SALAME algorithm
      */
     bool AnySpeciesSalame ();
+    amrex::Vector<std::string> m_names {"no_beam"}; /**< names of all beam containers */
+    amrex::Vector<BeamParticleContainer> m_all_beams; /**< contains all beam containers */
 
 private:
 
-    amrex::Vector<BeamParticleContainer> m_all_beams; /**< contains all beam containers */
-    amrex::Vector<std::string> m_names {"no_beam"}; /**< names of all beam containers */
     int m_nbeams {0}; /**< number of beam containers */
     /** number of real particles per beam, as opposed to ghost particles */
     amrex::Vector<amrex::Long> m_n_real_particles;

src/particles/collisions/ComputeTemperature.H

@@ -13,7 +13,7 @@ AMREX_GPU_HOST_DEVICE
 T_R ComputeTemperature (
     T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I,
     T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT psi,
-    T_R const m, T_R const clight, T_R const inv_c2 )
+    T_R const m, T_R const clight, T_R const inv_c2, bool is_beam_coll )
 {
     using namespace amrex::literals;
 
@@ -27,9 +27,11 @@ T_R ComputeTemperature (
     for (int i = Is; i < (int) Ie; ++i)
     {
         // particle's Lorentz factor
-        gm = (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2
-              + psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] );
-        const amrex::Real uz = clight * (gm - psi[I[i]]);
+        gm = is_beam_coll ? std::sqrt( 1._rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]]
+                                                + psi[I[i]]*psi[I[i]])*inv_c2 )
+                          : (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2
+                                       + psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] );
+        const amrex::Real uz = is_beam_coll ? psi[I[i]] : clight * (gm - psi[I[i]]);
         us = ( ux[ I[i] ] * ux[ I[i] ] +
                uy[ I[i] ] * uy[ I[i] ] +
                uz         * uz);

src/particles/collisions/CoulombCollision.H

@@ -2,25 +2,29 @@
 #define HIPACE_COULOMB_COLLISION_H_
 
 #include "particles/plasma/PlasmaParticleContainer.H"
+#include "particles/beam/BeamParticleContainer.H"
 
 #include <AMReX_DenseBins.H>
 #include <AMReX_REAL.H>
 #include <AMReX_ParmParse.H>
 
 /**
- * \brief This class handles Coulomb collisions between 2 plasma species (can be the same).
+ * \brief This class handles Coulomb collisions between 2 particle species
+ * (can be plasma-plasma or beam-plasma, the species can be the same)
  */
 class CoulombCollision
 {
 public:
     int  m_species1_index {-1};
     int  m_species2_index {-1};
+    int  m_nbeams {0};
     bool m_isSameSpecies {false};
     amrex::Real m_CoulombLog {-1.};
 
     /** Constructor */
     CoulombCollision(
-        const std::vector<std::string>& species_names,
+        const std::vector<std::string>& plasma_species_names,
+        const std::vector<std::string>& beam_species_names,
         std::string const collision_name);
 
     /**
@@ -37,11 +41,30 @@ public:
      *            Coulomb logarithm is deduced from the plasma temperature, measured in each cell.
      * \param[in] background_density_SI background plasma density (only needed for normalized units)
      **/
-    static void doCoulombCollision (
+    static void doPlasmaPlasmaCoulombCollision (
         int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1,
         PlasmaParticleContainer& species2, bool is_same_species, amrex::Real CoulombLog,
         amrex::Real background_density_SI);
 
+    /**
+     * \brief Perform Coulomb collisions of a beam with a plasma species over a push by one beam time step
+     *        Particles of both species are sorted per cell, paired, and collided pairwise.
+     *
+     * \param[in] lev MR level
+     * \param[in] bx transverse box (plasma particles will be sorted per-cell on this box)
+     * \param[in] geom corresponding geometry object
+     * \param[in] islice_local local slice on which collisions are calculated
+     * \param[in,out] species1 beam species
+     * \param[in,out] species2 plasma species
+     * \param[in] CoulombLog Value of the Coulomb logarithm used for the collisions. If <0, the
+     *            Coulomb logarithm is deduced from the plasma temperature, measured in each cell.
+     * \param[in] background_density_SI background plasma density (only needed for normalized units)
+     **/
+    static void doBeamPlasmaCoulombCollision (
+        int lev, const amrex::Box& bx, const amrex::Geometry& geom, int islice_local,
+        BeamParticleContainer& species1, PlasmaParticleContainer& species2, amrex::Real CoulombLog,
+        amrex::Real background_density_SI);
+
 };
 
 #endif // HIPACE_COULOMB_COLLISION_H_