Add beam-plasma collisions

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

@@ -207,6 +207,132 @@ Time step
     Only used when using adaptive time step (see ``hipace.dt`` above).
     Threshold beam momentum, below which the time step is not decreased (to avoid arbitrarily small time steps).
 
+* ``hipace.normalized_units`` (`bool`) optional (default `0`)
+    Using normalized units in the simulation.
+
+* ``hipace.verbose`` (`int`) optional (default `0`)
+    Level of verbosity.
+
+      * ``hipace.verbose = 1``, prints only the time steps, which are computed.
+
+      * ``hipace.verbose = 2`` additionally prints the number of iterations in the
+        predictor-corrector loop, as well as the B-Field error at each slice.
+
+      * ``hipace.verbose = 3`` also prints the number of particles, which violate the quasi-static
+        approximation and were neglected at each slice. It prints the number of ionized particles,
+        if ionization occurred. It also adds additional information if beams
+        are read in from file.
+
+* ``hipace.do_device_synchronize`` (`int`) optional (default `0`)
+    Level of synchronization on GPU.
+
+      * ``hipace.do_device_synchronize = 0``, synchronization happens only when necessary.
+
+      * ``hipace.do_device_synchronize = 1``, synchronizes most functions (all that are profiled
+        via ``HIPACE_PROFILE``)
+
+      * ``hipace.do_device_synchronize = 2`` additionally synchronizes low-level functions (all that
+        are profiled via ``HIPACE_DETAIL_PROFILE``)
+
+* ``hipace.depos_order_xy`` (`int`) optional (default `2`)
+    Transverse particle shape order. Currently, `0,1,2,3` are implemented.
+
+* ``hipace.depos_order_z`` (`int`) optional (default `0`)
+    Longitudinal particle shape order. Currently, only `0` is implemented.
+
+* ``hipace.depos_derivative_type`` (`int`) optional (default `2`)
+    Type of derivative used in explicit deposition. `0`: analytic, `1`: nodal, `2`: centered
+
+* ``hipace.outer_depos_loop`` (`bool`) optional (default `0`)
+    If the loop over depos_order is included in the loop over particles.
+
+* ``hipace.beam_injection_cr`` (`integer`) optional (default `1`)
+    Using a temporary coarsed grid for beam particle injection for a fixed particle-per-cell beam.
+    For very high-resolution simulations, where the number of grid points (`nx*ny*nz`)
+    exceeds the maximum `int (~2e9)`, it enables beam particle injection, which would
+    fail otherwise. As an example, a simulation with `2048 x 2048 x 2048` grid points
+    requires ``hipace.beam_injection_cr = 8``.
+
+* ``hipace.do_beam_jx_jy_deposition`` (`bool`) optional (default `1`)
+    Using the default, the beam deposits all currents ``Jx``, ``Jy``, ``Jz``. Using
+    ``hipace.do_beam_jx_jy_deposition = 0`` disables the transverse current deposition of the beams.
+
+* ``hipace.boxes_in_z`` (`int`) optional (default `1`)
+    Number of boxes along the z-axis. In serial runs, the arrays for 3D IO can easily exceed the
+    memory of a GPU. Using multiple boxes reduces the memory requirements by the same factor.
+    This option is only available in serial runs, in parallel runs, please use more GPU to achieve
+    the same effect.
+
+* ``hipace.openpmd_backend`` (`string`) optional (default `h5`)
+    OpenPMD backend. This can either be ``h5``, ``bp``, or ``json``. The default is chosen by what is
+    available. If both Adios2 and HDF5 are available, ``h5`` is used. Note that ``json`` is extremely
+    slow and is not recommended for production runs.
+
+* ``hipace.file_prefix`` (`string`) optional (default `diags/hdf5/`)
+    Path of the output.
+
+* ``hipace.do_tiling`` (`bool`) optional (default `true`)
+    Whether to use tiling, when running on CPU.
+    Currently, this option only affects plasma operations (gather, push and deposition).
+    The tile size can be set with ``plasmas.sort_bin_size``.
+
+* ``hipace.comms_buffer_on_gpu`` (`bool`) optional (default `0`)
+    Whether the buffers used for MPI communication should be allocated on the GPU (device memory).
+    By default they will be allocated on the CPU (pinned memory).
+    Setting this option to `1` is necessary to take advantage of GPU-Enabled MPI, however for this
+    additional enviroment variables need to be set depending on the system.
+
+* ``hipace.do_beam_jz_minus_rho`` (`bool`) optional (default `0`)
+    Whether the beam contribution to :math:`j_z-c\rho` is calculated and used when solving for Psi (used to caculate the transverse fields Ex-By and Ey+Bx).
+    if 0, this term is assumed to be 0 (a good approximation for an ultra-relativistic beam in the z direction with small transverse momentum).
+
+* ``hipace.deposit_rho`` (`bool`) optional (default `0`)
+    If the charge density ``rho`` of the plasma should be deposited so that it is available as a diagnostic.
+    Otherwise only ``rhomjz`` equal to :math:`\rho-j_z/c` will be available.
+    If ``rho`` is explicitly mentioned in ``diagnostic.field_data``, then the default will become `1`.
+
+* ``hipace.salame_n_iter`` (`int`) optional (default `3`)
+    Number of iterations the SALAME algorithm should do when it is used.
+
+* ``hipace.salame_do_advance`` (`bool`) optional (default `1`)
+    Whether the SALAME algorithm should calculate the SALAME-beam-only Ez field
+    by advancing plasma (if `1`) particles or by approximating it using the chi field (if `0`).
+
+* ``hipace.salame_Ez_target(zeta,zeta_initial,Ez_initial)`` (`string`) optional (default `Ez_initial`)
+    Parser function to specify the target Ez field at the witness beam for SALAME.
+    ``zeta``: position of the Ez field to set.
+    ``zeta_initial``: position where the SALAME algorithm first started.
+    ``Ez_initial``: field value at `zeta_initial`.
+    For `zeta` equal to `zeta_initial`, the function should return `Ez_initial`.
+    The default value of this function corresponds to a flat Ez field at the position of the SALAME beam.
+    Note: `zeta` is always less than or equal to `zeta_initial` and `Ez_initial` is typically below zero for electron beams.
+
+* ``hipace.background_density_SI`` (`float`) optional
+    Background plasma density in SI units. Certain physical modules (collisions, ionization, radiation reactions) depend on the actual background density.
+    Hence, in normalized units, they can only be included, if a background plasma density in SI units is provided using this input parameter.
+
+Binary collisions
+-----------------
+
+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.
+
+* ``plasmas.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 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.
+    If not specified, the Coulomb logarithm is determined from the temperature in each cell.
+
 Field solver parameters
 -----------------------
 
@@ -810,15 +936,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).
+    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_