Implement binary Coulomb collisions between plasma species

CMakeLists.txt

@@ -225,6 +225,12 @@ if(BUILD_TESTING)
                  WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
         )
 
+        add_test(NAME collisions.SI.1Rank
+                 COMMAND ${HiPACE_SOURCE_DIR}/tests/collisions.SI.1Rank.sh
+                         $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
+                 WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+        )
+
         #add_test(NAME hosing.2Rank
         #         COMMAND ${HiPACE_SOURCE_DIR}/tests/hosing.2Rank.sh
         #                 $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}

docs/source/run/parameters.rst

@@ -276,6 +276,24 @@ plasma parameters for each plasma are specified via `<plasma name>.plasma_proper
     arrays of size ``sort_bin_size`` (+ guard cells) that are atomic-added to the main current
     arrays.
 
+Binary collisions for plasma species
+------------------------------------
+
+WARNING: this module is in development. Currently only support electron-electron collisions in SI units.
+
+HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX, between plasma species.
+
+* ``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).
+
+* ``<collision name>.species`` (two `strings`) optional
+    The name of the two plasma species for which collisions should be included.
+
+* ``<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.
+
 Beam parameters
 ---------------
 

src/Hipace.cpp

@@ -597,7 +597,10 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox,
 
             FillDiagnostics(lev, islice);
 
+            m_multi_plasma.doCoulombCollision(lev, bx, geom[lev]);
+
             m_multi_plasma.DoFieldIonization(lev, geom[lev], m_fields);
+
             if (m_multi_plasma.IonizationOn() && m_do_tiling) m_multi_plasma.TileSort(bx, geom[lev]);
 
         } // end for (int isubslice = nsubslice-1; isubslice >= 0; --isubslice)

src/particles/CMakeLists.txt

@@ -17,6 +17,7 @@ target_sources(HiPACE
     BoxSort.cpp
     MultiBeam.cpp
     MultiPlasma.cpp
+    CoulombCollision.cpp
 )
 
 add_subdirectory(deposition)

src/particles/ComputeTemperature.H

@@ -0,0 +1,52 @@
+/* Copyright 2019-2020 Andrew Myers, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef HIPACE_COMPUTE_TEMPERATURE_H_
+#define HIPACE_COMPUTE_TEMPERATURE_H_
+
+template <typename T_index, typename T_R>
+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 )
+{
+    using namespace amrex::literals;
+    T_R constexpr c2i = 1._rt/(PhysConstSI::c * PhysConstSI::c);
+
+    int N = Ie - Is;
+    if ( N == 0 ) { return T_R(0.0); }
+
+    T_R vx = T_R(0.0);    T_R vy = T_R(0.0);
+    T_R vz = T_R(0.0);    T_R vs = T_R(0.0);
+    T_R gm = T_R(0.0);    T_R us = T_R(0.0);
+
+    for (int i = Is; i < Ie; ++i)
+    {
+        // particle's normalized vector potential
+        const amrex::Real psii = psi[I[i]] *
+            PhysConstSI::q_e / (PhysConstSI::m_e * PhysConstSI::c * PhysConstSI::c);
+        // particle's Lorentz factor
+        gm = (1.0_rt + ux[I[i]]*ux[I[i]]*c2i + uy[I[i]]*uy[I[i]]*c2i
+              + (psii+1._rt)*(psii+1._rt)) / (2.0_rt * (psii+1._rt) );
+        const amrex::Real uz = PhysConstSI::c * (gm - psii - 1._rt);
+        us = ( ux[ I[i] ] * ux[ I[i] ] +
+               uy[ I[i] ] * uy[ I[i] ] +
+               uz         * uz);
+        vx += ux[ I[i] ] / gm;
+        vy += uy[ I[i] ] / gm;
+        vz += uz         / gm;
+        vs += us / gm / gm;
+    }
+
+    vx = vx / N;    vy = vy / N;
+    vz = vz / N;    vs = vs / N;
+
+    return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz));
+}
+
+#endif // HIPACE_COMPUTE_TEMPERATURE_H_

src/particles/CoulombCollision.H

@@ -0,0 +1,45 @@
+#ifndef HIPACE_COULOMB_COLLISION_H_
+#define HIPACE_COULOMB_COLLISION_H_
+
+#include "PlasmaParticleContainer.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).
+ */
+class CoulombCollision
+{
+public:
+    int  m_species1_index;
+    int  m_species2_index;
+    bool m_isSameSpecies;
+    amrex::Real m_CoulombLog;
+
+    /** Constructor */
+    CoulombCollision(
+        const std::vector<std::string>& species_names,
+        std::string const collision_name);
+
+    /**
+     * \brief Perform Coulomb collisions of plasma species over longitudinal push by 1 cell.
+     *        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,out] species1 first plasma species
+     * \param[in,out] species2 second plasma species
+     * \param[in] is_same_species whether both species are the same (intra-species collisions)
+     * \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.
+     **/
+    static void doCoulombCollision (
+        int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1,
+        PlasmaParticleContainer& species2, const bool is_same_species, const amrex::Real CoulombLog);
+
+};
+
+#endif // HIPACE_COULOMB_COLLISION_H_

src/particles/CoulombCollision.cpp

@@ -0,0 +1,188 @@
+#include "CoulombCollision.H"
+#include "ShuffleFisherYates.H"
+#include "TileSort.H"
+#include "ElasticCollisionPerez.H"
+#include "utils/HipaceProfilerWrapper.H"
+
+CoulombCollision::CoulombCollision(
+    const std::vector<std::string>& species_names,
+    std::string const collision_name)
+{
+    using namespace amrex::literals;
+
+    // TODO: ionization level
+    // DOTO: Fix dt
+
+    // read collision species
+    std::vector<std::string> collision_species;
+    amrex::ParmParse pp(collision_name);
+    pp.getarr("species", collision_species);
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(collision_species.size() == 2,
+    "Collision species must name exactly two species.");
+
+    // default Coulomb log, if < 0, will be computed automatically
+    m_CoulombLog = -1.0;
+    pp.query("CoulombLog", m_CoulombLog);
+
+    for (int i=0; i<species_names.size(); i++)
+    {
+        if (species_names[i] == collision_species[0]) m_species1_index = i;
+        if (species_names[i] == collision_species[1]) m_species2_index = i;
+    }
+
+    m_isSameSpecies = collision_species[0] == collision_species[1] ? true : false;
+}
+
+void
+CoulombCollision::doCoulombCollision (
+    int lev, const amrex::Box& bx, const amrex::Geometry& geom, PlasmaParticleContainer& species1,
+    PlasmaParticleContainer& species2, const bool is_same_species, const amrex::Real CoulombLog)
+{
+    HIPACE_PROFILE("CoulombCollision::doCoulombCollision()");
+    AMREX_ALWAYS_ASSERT(lev == 0);
+    using namespace amrex::literals;
+
+    if ( is_same_species ) // species_1 == species_2
+    {
+        // Logically particles per-cell, and return indices of particles in each cell
+        PlasmaBins bins1 = findParticlesInEachTile(lev, bx, 1, species1, geom);
+        int const n_cells = bins1.numBins();
+
+        // Counter to check there is only 1 box
+        int count = 0;
+        for (PlasmaParticleIterator pti(species1, lev); pti.isValid(); ++pti) {
+
+            // Get particles SoA and AoS data
+            auto& aos1 = pti.GetArrayOfStructs();
+            const auto& pos1 = aos1.begin();
+            auto& soa1 = pti.GetStructOfArrays();
+            amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux).data();
+            amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy).data();
+            amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi).data();
+            const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data();
+            PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
+            PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
+            amrex::Real q1 = species1.GetCharge();
+            amrex::Real m1 = species1.GetMass();
+
+            const amrex::Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2);
+            const amrex::Real dt = geom.CellSize(2)/PhysConstSI::c;
+
+            amrex::ParallelForRNG(
+                n_cells,
+                [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
+                {
+                    // The particles from species1 that are in the cell `i_cell` are
+                    // given by the `indices_1[cell_start_1:cell_stop_1]`
+                    PlasmaBins::index_type const cell_start1 = offsets1[i_cell];
+                    PlasmaBins::index_type const cell_stop1  = offsets1[i_cell+1];
+                    PlasmaBins::index_type const cell_half1 = (cell_start1+cell_stop1)/2;
+
+                    if ( cell_stop1 - cell_start1 <= 1 ) return;
+                    // Do not collide if there is only one particle in the cell
+                    // shuffle
+                    ShuffleFisherYates(
+                        indices1, cell_start1, cell_half1, engine );
+
+                    // TODO: FIX DT
+                    // Call the function in order to perform collisions
+                    ElasticCollisionPerez(
+                        cell_start1, cell_half1,
+                        cell_half1, cell_stop1,
+                        indices1, indices1,
+                        ux1, uy1, psi1, ux1, uy1, psi1, w1, w1,
+                        q1, q1, m1, m1, -1.0_rt, -1.0_rt,
+                        dt, CoulombLog, dV, engine );
+                }
+                );
+            count++;
+        }
+        AMREX_ALWAYS_ASSERT(count == 1);
+
+    } else {
+
+        // Logically particles per-cell, and return indices of particles in each cell
+        PlasmaBins bins1 = findParticlesInEachTile(lev, bx, 1, species1, geom);
+        PlasmaBins bins2 = findParticlesInEachTile(lev, bx, 1, species2, geom);
+
+        int const n_cells = bins1.numBins();
+
+        // Counter to check there is only 1 box
+        int count = 0;
+        for (PlasmaParticleIterator pti(species1, lev); pti.isValid(); ++pti) {
+
+            // Get particles SoA and AoS data for species 1
+            auto& aos1 = pti.GetArrayOfStructs();
+            const auto& pos1 = aos1.begin();
+            auto& soa1 = pti.GetStructOfArrays();
+            amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux).data();
+            amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy).data();
+            amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi).data();
+            const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data();
+            PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
+            PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
+            amrex::Real q1 = species1.GetCharge();
+            amrex::Real m1 = species1.GetMass();
+
+            // Get particles SoA and AoS data for species 1
+            auto index = std::make_pair(pti.index(), pti.LocalTileIndex());
+            auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex());
+            auto& aos2 = ptile2.GetArrayOfStructs();
+            const auto& pos2 = aos2.begin();
+            auto& soa2 = ptile2.GetStructOfArrays();
+            amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux).data();
+            amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy).data();
+            amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi).data();
+            const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data();
+            PlasmaBins::index_type * const indices2 = bins2.permutationPtr();
+            PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr();
+            amrex::Real q2 = species2.GetCharge();
+            amrex::Real m2 = species2.GetMass();
+
+            const amrex::Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2);
+            const amrex::Real dt = geom.CellSize(2)/PhysConstSI::c;
+
+            // Extract particles in the tile that `mfi` points to
+            // ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi);
+            // ParticleTileType& ptile_2 = species_2->ParticlesAt(lev, mfi);
+            // Loop over cells, and collide the particles in each cell
+
+            // Loop over cells
+            amrex::ParallelForRNG(
+                n_cells,
+                [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
+                {
+                    // The particles from species1 that are in the cell `i_cell` are
+                    // given by the `indices_1[cell_start_1:cell_stop_1]`
+                    PlasmaBins::index_type const cell_start1 = offsets1[i_cell];
+                    PlasmaBins::index_type const cell_stop1  = offsets1[i_cell+1];
+                    // Same for species 2
+                    PlasmaBins::index_type const cell_start2 = offsets2[i_cell];
+                    PlasmaBins::index_type const cell_stop2  = offsets2[i_cell+1];
+
+                    // ux from species1 can be accessed like this:
+                    // ux_1[ indices_1[i] ], where i is between
+                    // cell_start_1 (inclusive) and cell_start_2 (exclusive)
+
+                    // Do not collide if one species is missing in the cell
+                    if ( cell_stop1 - cell_start1 < 1 ||
+                         cell_stop2 - cell_start2 < 1 ) return;
+                    // shuffle
+                    ShuffleFisherYates(indices1, cell_start1, cell_stop1, engine);
+                    ShuffleFisherYates(indices2, cell_start2, cell_stop2, engine);
+
+                    // TODO: FIX DT.
+                    // Call the function in order to perform collisions
+                    ElasticCollisionPerez(
+                        cell_start1, cell_stop1, cell_start2, cell_stop2,
+                        indices1, indices2,
+                        ux1, uy1, psi1, ux2, uy2, psi2, w1, w2,
+                        q1, q2, m1, m2, -1.0_rt, -1.0_rt,
+                        dt, CoulombLog, dV, engine );
+                }
+                );
+            count++;
+        }
+        AMREX_ALWAYS_ASSERT(count == 1);
+    }
+};