apply code review changes

Co-authored-by: RemiLehe <[email protected]>

Docs/source/theory/collisions.rst

@@ -5,13 +5,13 @@ Collisions
 
 WarpX includes several different models to capture collisional processes
 including collisions between kinetic particles (Coulomb collisions, DSMC,
-nuclear fusion) as well as collisions between kinetic particles and a fluid
-background species (MCC, background stopping).
+nuclear fusion) as well as collisions between kinetic particles and a fixed
+(i.e. non-evolving) background species (MCC, background stopping).
 
 .. _theory-collisions-mcc:
 
-Monte Carlo Collisions
-----------------------
+Background Monte Carlo Collisions (MCC)
+---------------------------------------
 
 Several types of collisions between simulation particles and a neutral
 background gas are supported including elastic scattering, back scattering,
@@ -54,10 +54,21 @@ Once a particle is selected for a specific collision process, that process deter
 
 .. _theory-collisions-dsmc:
 
-Direct Simulation Monte Carlo
------------------------------
+Direct Simulation Monte Carlo (DSMC)
+------------------------------------
 
-(Under construction)
+The algorithm by which binary collisions are treated is outlined below. The
+description assumes collisions between different species.
+
+1. Particles from both species are sorted by grid-cells.
+2. The order of the particles in each cell is shuffled.
+3. Within each cell, particles are paired to form collision partners. Particles
+   of the species with fewer members in a given cell is split in half so that
+   each particle has exactly one partner of the other species.
+4. Each collision pair is considered for a collision using the same logic as in
+   the MCC description above.
+5. Particles that are chosen for collision are scattered according to the
+   selected collision process.
 
 Scattering processes
 --------------------
@@ -72,7 +83,7 @@ with the neutral velocity and vice-versa.
 Elastic scattering
 ^^^^^^^^^^^^^^^^^^
 
-The ``elastic`` option uses hard-sphere scattering, with a differential
+The ``elastic`` option uses isotropic scattering, i.e., with a differential
 cross section that is independent of angle.
 This scattering process as well as the ones below that relate to it, are all
 performed in the center-of-momentum (COM) frame. Designating the COM velocity of
@@ -120,7 +131,7 @@ The process is also the same as for elastic scattering except the excitation ene
 Benchmarks
 ----------
 
-See the :ref:`MCC example <examples-mcc-turner>` for a bencmark of the MCC
+See the :ref:`MCC example <examples-mcc-turner>` for a benchmark of the MCC
 implementation against literature results.
 
 Particle cooling due to elastic collisions

Docs/source/usage/python.rst

@@ -121,6 +121,8 @@ Other operations related to particles:
 
 .. autoclass:: pywarpx.picmi.CoulombCollisions
 
+.. autoclass:: pywarpx.picmi.DSMCCollisions
+
 .. autoclass:: pywarpx.picmi.MCCCollisions
 
 .. autoclass:: pywarpx.picmi.FieldIonization

Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H

@@ -42,17 +42,16 @@ namespace BinaryCollisionUtils{
     /**
      * \brief Return (relativistic) collision energy, collision speed and
      * Lorentz factor for transforming between the lab and center-of-momentum
-     * frames. Note the use of `double` for energy since this calculation is
-     * prone to error with single precision.
+     * frames.
      */
     AMREX_GPU_HOST_DEVICE AMREX_INLINE
     void get_collision_parameters (
         const amrex::ParticleReal& u1x, const amrex::ParticleReal& u1y,
         const amrex::ParticleReal& u1z, const amrex::ParticleReal& u2x,
         const amrex::ParticleReal& u2y, const amrex::ParticleReal& u2z,
         const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
-        amrex::ParticleReal& E_coll, amrex::ParticleReal& v_coll,
-        amrex::ParticleReal& lab_to_COM_factor )
+        amrex::ParticleReal& E_kin_COM, amrex::ParticleReal& v_rel_COM,
+        amrex::ParticleReal& lab_to_COM_lorentz_factor )
     {
         // General notations in this function:
         //     x_sq denotes the square of x
@@ -86,6 +85,8 @@ namespace BinaryCollisionUtils{
         );
 
         // Total energy in the lab frame
+        // Note the use of `double` for energy since this calculation is
+        // prone to error with single precision.
         const auto m1_dbl = static_cast<double>(m1);
         const auto m2_dbl = static_cast<double>(m2);
         const double E_lab = (m1_dbl * g1 + m2_dbl * g2) * c_sq;
@@ -95,7 +96,7 @@ namespace BinaryCollisionUtils{
 
         // Kinetic energy in the center of mass frame
         const double E_star = std::sqrt(E_star_sq);
-        E_coll = static_cast<amrex::ParticleReal>(E_star - (m1_dbl + m2_dbl)*c_sq);
+        E_kin_COM = static_cast<amrex::ParticleReal>(E_star - (m1_dbl + m2_dbl)*c_sq);
 
         // Square of the norm of the momentum of one of the particles in the center of mass frame
         // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle
@@ -112,14 +113,16 @@ namespace BinaryCollisionUtils{
         const auto g2_star = std::sqrt(one_pr + p_star_sq / static_cast<amrex::ParticleReal>(m2_sq*c_sq));
 
         // relative velocity in the center of mass frame
-        v_coll = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + one_pr/(m2*g2_star));
+        v_rel_COM = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + one_pr/(m2*g2_star));
 
         // Cross sections and relative velocity are computed in the center of mass frame.
-        // On the other hand, the particle densities (weight over volume) in the lab frame are used. To
-        // take into account this discrepancy, we need to multiply the fusion probability by the ratio
-        // between the Lorentz factors in the COM frame and the Lorentz factors in the lab frame
-        // (see Perez et al., Phys.Plasmas.19.083104 (2012))
-        lab_to_COM_factor = g1_star*g2_star/static_cast<amrex::ParticleReal>(g1*g2);
+        // On the other hand, the particle densities (weight over volume) in the lab frame are used.
+        // To take this disrepancy into account, it is needed to multiply the
+        // collision probability by the ratio between the Lorentz factors in the
+        // COM frame and the Lorentz factors in the lab frame (see
+        // Perez et al., Phys.Plasmas.19.083104 (2012)). The correction factor
+        // is calculated here.
+        lab_to_COM_lorentz_factor = g1_star*g2_star/static_cast<amrex::ParticleReal>(g1*g2);
     }
 }
 

Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H

@@ -27,7 +27,7 @@
 
 
 /**
- * \brief This class perfoms DSMC (direct simulation Monte Carlo) collisions
+ * \brief This class performs DSMC (direct simulation Monte Carlo) collisions
  * within a cell. Particles are paired up and for each pair a stochastic process
  * determines whether a collision occurs. The algorithm is similar to the one
  * used for binary Coulomb collisions and the nuclear fusion module.

Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp

@@ -19,45 +19,45 @@
 
 ParticleCreationFunc::ParticleCreationFunc (const std::string collision_name,
                                             MultiParticleContainer const * const mypc)
-    {
-        const amrex::ParmParse pp_collision_name(collision_name);
+{
+    const amrex::ParmParse pp_collision_name(collision_name);
 
-        m_collision_type = BinaryCollisionUtils::get_collision_type(collision_name, mypc);
+    m_collision_type = BinaryCollisionUtils::get_collision_type(collision_name, mypc);
 
-        if (m_collision_type == CollisionType::ProtonBoronToAlphasFusion)
-        {
-            // Proton-Boron fusion only produces alpha particles
-            m_num_product_species = 1;
-            // Proton-Boron fusion produces 3 alpha particles per fusion reaction
-            m_num_products_host.push_back(3);
+    if (m_collision_type == CollisionType::ProtonBoronToAlphasFusion)
+    {
+        // Proton-Boron fusion only produces alpha particles
+        m_num_product_species = 1;
+        // Proton-Boron fusion produces 3 alpha particles per fusion reaction
+        m_num_products_host.push_back(3);
 #ifndef AMREX_USE_GPU
-            // On CPU, the device vector can be filled immediatly
-            m_num_products_device.push_back(3);
+        // On CPU, the device vector can be filled immediately
+        m_num_products_device.push_back(3);
 #endif
-        }
-        else if ((m_collision_type == CollisionType::DeuteriumTritiumToNeutronHeliumFusion)
-              || (m_collision_type == CollisionType::DeuteriumDeuteriumToProtonTritiumFusion)
-              || (m_collision_type == CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion))
-        {
-            m_num_product_species = 2;
-            m_num_products_host.push_back(1);
-            m_num_products_host.push_back(1);
+    }
+    else if ((m_collision_type == CollisionType::DeuteriumTritiumToNeutronHeliumFusion)
+             || (m_collision_type == CollisionType::DeuteriumDeuteriumToProtonTritiumFusion)
+             || (m_collision_type == CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion))
+    {
+        m_num_product_species = 2;
+        m_num_products_host.push_back(1);
+        m_num_products_host.push_back(1);
 #ifndef AMREX_USE_GPU
-            // On CPU, the device vector can be filled immediately
-            m_num_products_device.push_back(1);
-            m_num_products_device.push_back(1);
+        // On CPU, the device vector can be filled immediately
+        m_num_products_device.push_back(1);
+        m_num_products_device.push_back(1);
 #endif
-        }
-        else
-        {
-          WARPX_ABORT_WITH_MESSAGE("Unknown collision type in ParticleCreationFunc");
-        }
+    }
+    else
+    {
+        WARPX_ABORT_WITH_MESSAGE("Unknown collision type in ParticleCreationFunc");
+    }
 
 #ifdef AMREX_USE_GPU
-        m_num_products_device.resize(m_num_product_species);
-        amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_num_products_host.begin(),
-                              m_num_products_host.end(),
-                              m_num_products_device.begin());
-        amrex::Gpu::streamSynchronize();
+     m_num_products_device.resize(m_num_product_species);
+     amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_num_products_host.begin(),
+                           m_num_products_host.end(),
+                           m_num_products_device.begin());
+     amrex::Gpu::streamSynchronize();
 #endif
-    }
+}