Add DSMC module

Docs/source/refs.bib

@@ -244,3 +244,31 @@ @inproceedings{SandbergIPAC23
     url = {https://indico.jacow.org/event/41/contributions/2276},
     language = {English}
 }
+
+@article{Perez2012,
+    author = {Pérez, F. and Gremillet, L. and Decoster, A. and Drouin, M. and Lefebvre, E.},
+    title = "{Improved modeling of relativistic collisions and collisional ionization in particle-in-cell codes}",
+    journal = {Physics of Plasmas},
+    volume = {19},
+    number = {8},
+    pages = {083104},
+    year = {2012},
+    month = {08},
+    issn = {1070-664X},
+    doi = {10.1063/1.4742167},
+    url = {https://doi.org/10.1063/1.4742167},
+    eprint = {https://pubs.aip.org/aip/pop/article-pdf/doi/10.1063/1.4742167/13891570/083104\_1\_online.pdf},
+}
+
+@article{Higginson2019,
+  doi = {10.1016/j.jcp.2019.03.020},
+  url = {https://doi.org/10.1016/j.jcp.2019.03.020},
+  year = {2019},
+  month = jul,
+  publisher = {Elsevier {BV}},
+  volume = {388},
+  pages = {439--453},
+  author = {Drew Pitney Higginson and Anthony Link and Andrea Schmidt},
+  title = {A pairwise nuclear fusion algorithm for weighted particle-in-cell plasma simulations},
+  journal = {Journal of Computational Physics}
+}

Docs/source/theory/collisions.rst

@@ -3,18 +3,36 @@
 Collisions
 ==========
 
-Monte Carlo Collisions
-----------------------
-
-The Monte Carlo collisions (MCC) module can be found in *Source/Particles/Collisions/BackgroundMCC/*.
-Several types of collisions between simulation particles and a neutral background gas are supported including elastic scattering, back scattering, charge exchange, excitation collisions and impact ionization.
-An instance of the class :cpp:class:`MCCProcess` is created for each type of collision included in a simulation. This class saves information about the type of collision and the collision cross-section as a function of energy.
-
-The so-called null collision strategy is used in order to minimize the computational burden of the MCC module.
-This strategy is standard in PIC-MCC and a detailed description can be found elsewhere, for example in :cite:t:`b-Birdsall1991`.
-In short the maximum collision probability is found over a sensible range of energies and is used to pre-select the appropriate number of macroparticles for collision consideration. Only these pre-selected particles are then individually considered for a collision based on their energy and the cross-sections of all the different collisional processes included.
-
-The MCC implementation assumes that the background neutral particles are **thermal**, and are moving at non-relativistic velocities in the lab frame. For each simulation particle considered for a collision, a velocity vector for a neutral particle is randomly chosen. The particle velocity is then boosted to the stationary frame of the neutral through a Galilean transformation. The energy of the collision is calculated using the particle utility function, ``ParticleUtils::getCollisionEnergy()``, as
+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 fixed
+(i.e. non-evolving) background species (MCC, background stopping).
+
+.. _theory-collisions-mcc:
+
+Background Monte Carlo Collisions (MCC)
+---------------------------------------
+
+Several types of collisions between simulation particles and a neutral
+background gas are supported including elastic scattering, back scattering,
+charge exchange, excitation collisions and impact ionization.
+
+The so-called null collision strategy is used in order to minimize the
+computational burden of the MCC module. This strategy is standard in PIC-MCC and
+a detailed description can be found elsewhere, for example in :cite:t:`b-Birdsall1991`.
+In short the maximum collision probability is found over a sensible range of
+energies and is used to pre-select the appropriate number of macroparticles for
+collision consideration. Only these pre-selected particles are then individually
+considered for a collision based on their energy and the cross-sections of all
+the different collisional processes included.
+
+The MCC implementation assumes that the background neutral particles are **thermal**,
+and are moving at non-relativistic velocities in the lab frame. For each
+simulation particle considered for a collision, a velocity vector for a neutral
+particle is randomly chosen given the user specified neutral temperature. The
+particle velocity is then boosted to the stationary frame of the neutral through
+a Galilean transformation. The energy of the collision is calculated using the
+particle utility function, ``ParticleUtils::getCollisionEnergy()``, as
 
     .. math::
 
@@ -23,19 +41,55 @@ The MCC implementation assumes that the background neutral particles are **therm
                  &= \frac{2Mmu^2}{M + m + \sqrt{M^2+m^2+2\gamma mM}}\frac{1}{\gamma + 1}
        \end{aligned}
 
-where :math:`u` is the speed of the particle as tracked in WarpX (i.e. :math:`u = \gamma v` with :math:`v` the particle speed), while :math:`m` and :math:`M` are the rest masses of the simulation and background species, respectively. The Lorentz factor is defined in the usual way, :math:`\gamma = \sqrt{1 + u^2/c^2}`. Note that if :math:`\gamma\to1` the above expression clearly reduces to the classical equation :math:`E_{coll} = \frac{1}{2}\frac{Mm}{M+m} u^2`. The collision cross-sections for all scattering processes are evaluated at the energy as calculated above.
+where :math:`u` is the speed of the particle as tracked in WarpX (i.e.
+:math:`u = \gamma v` with :math:`v` the particle speed), while :math:`m` and
+:math:`M` are the rest masses of the simulation and background species,
+respectively. The Lorentz factor is defined in the usual way,
+:math:`\gamma \def \sqrt{1 + u^2/c^2}`. Note that if :math:`\gamma\to1` the above
+expression reduces to the classical equation
+:math:`E_{coll} = \frac{1}{2}\frac{Mm}{M+m} u^2`. The collision cross-sections
+for all scattering processes are evaluated at the energy as calculated above.
 
 Once a particle is selected for a specific collision process, that process determines how the particle is scattered as outlined below.
 
+.. _theory-collisions-dsmc:
+
+Direct Simulation Monte Carlo (DSMC)
+------------------------------------
+
+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
+--------------------
+
 Charge exchange
 ^^^^^^^^^^^^^^^
 
-This is the simplest scattering process. Under charge exchange the simulation particle's velocity is simply replaced with the sampled velocity of the neutral particle. Charge exchange usually has a cooling effect on the ions in the simulation by exchanging energetic ions for lower energy neutrals.
+This process can occur when an ion and neutral (of the same species) collide
+and results in the exchange of an electron. The ion velocity is simply replaced
+with the neutral velocity and vice-versa.
 
 Elastic scattering
 ^^^^^^^^^^^^^^^^^^
 
-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 the particle as :math:`\vec{u}_c` and its labframe velocity as :math:`\vec{u}_l`, the transformation from lab frame to COM frame is done with a general Lorentz boost (see function ``ParticleUtils::doLorentzTransform()``):
+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
+the particle as :math:`\vec{u}_c` and its labframe velocity as :math:`\vec{u}_l`,
+the transformation from lab frame to COM frame is done with a general Lorentz
+boost (see function ``ParticleUtils::doLorentzTransform()``):
 
     .. math::
             \begin{bmatrix}
@@ -74,6 +128,12 @@ Excitation
 
 The process is also the same as for elastic scattering except the excitation energy cost is subtracted from the particle energy. This is done by reducing the velocity before a scattering angle is chosen.
 
+Benchmarks
+----------
+
+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/parameters.rst

@@ -1663,6 +1663,7 @@ Collision models
 ----------------
 
 WarpX provides several particle collision models, using varying degrees of approximation.
+Details about the collision models can be found in the :ref:`theory section <theory-collisions>`.
 
 * ``collisions.collision_names`` (`strings`, separated by spaces)
     The name of each collision type.
@@ -1672,29 +1673,29 @@ WarpX provides several particle collision models, using varying degrees of appro
 * ``<collision_name>.type`` (`string`) optional
     The type of collision. The types implemented are:
 
-    - ``pairwisecoulomb`` for pairwise Coulomb collisions, the default if unspecified.
+    - ``pairwisecoulomb`` for pair-wise Coulomb collisions, the default if unspecified.
       This provides a pair-wise relativistic elastic Monte Carlo binary Coulomb collision model,
-      following the algorithm given by `Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`__.
+      following the algorithm given by :cite:t:`param-Perez2012`.
       When the RZ mode is used, `warpx.n_rz_azimuthal_modes` must be set to 1 at the moment,
       since the current implementation of the collision module assumes axisymmetry.
     - ``nuclearfusion`` for fusion reactions.
-      This implements the pair-wise fusion model by `Higginson et al. (JCP 388, 439-453, 2019) <https://doi.org/10.1016/j.jcp.2019.03.020>`__.
+      This implements the pair-wise fusion model by :cite:t:`param-Higginson2019`.
       Currently, WarpX supports deuterium-deuterium, deuterium-tritium, deuterium-helium and proton-boron fusion.
       When initializing the reactant and product species, you need to use ``species_type`` (see the documentation
       for this parameter), so that WarpX can identify the type of reaction to use.
       (e.g. ``<species_name>.species_type = 'deuterium'``)
+    - ``dsmc`` for pair-wise, non-Coulomb collisions between kinetic species.
+      This is a "direct simulation Monte Carlo" treatment of collisions between
+      kinetic species. See :ref:`DSMC section <theory-collisions-dsmc>`.
     - ``background_mcc`` for collisions between particles and a neutral background.
       This is a relativistic Monte Carlo treatment for particles colliding
-      with a neutral background gas. The implementation follows the so-called
-      null collision strategy discussed for example in `Birdsall (IEEE Transactions on
-      Plasma Science, vol. 19, no. 2, pp. 65-85, 1991) <https://ieeexplore.ieee.org/document/106800>`_.
-      See also :ref:`collisions section <theory-collisions>`.
+      with a neutral background gas. See :ref:`MCC section <theory-collisions-mcc>`.
     - ``background_stopping`` for slowing of ions due to collisions with electrons or ions.
       This implements the approximate formulae as derived in Introduction to Plasma Physics,
       from Goldston and Rutherford, section 14.2.
 
 * ``<collision_name>.species`` (`strings`)
-    If using ``pairwisecoulomb`` or ``nuclearfusion``, this should be the name(s) of the species,
+    If using ``dsmc``, ``pairwisecoulomb`` or ``nuclearfusion``, this should be the name(s) of the species,
     between which the collision will be considered. (Provide only one name for intra-species collisions.)
     If using ``background_mcc`` or ``background_stopping`` type this should be the name of the
     species for which collisions with a background will be included.
@@ -1717,7 +1718,7 @@ WarpX provides several particle collision models, using varying degrees of appro
     :math:`A` is the mass number.
     If this is not provided, or if a non-positive value is provided,
     a Coulomb logarithm will be computed automatically according to the algorithm in
-    `Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`__.
+    :cite:t:`param-Perez2012`.
 
 * ``<collision_name>.fusion_multiplier`` (`float`) optional.
     Only for ``nuclearfusion``.
@@ -1728,8 +1729,8 @@ WarpX provides several particle collision models, using varying degrees of appro
     More specifically, in a fusion reaction between two macroparticles with weight ``w_1`` and ``w_2``,
     the weight of the product macroparticles will be ``min(w_1,w_2)/fusion_multiplier``.
     (And the weights of the reactant macroparticles are reduced correspondingly after the reaction.)
-    See `Higginson et al. (JCP 388, 439-453, 2019) <https://doi.org/10.1016/j.jcp.2019.03.020>`__
-    for more details. The default value of ``fusion_multiplier`` is 1.
+    See :cite:t:`param-Higginson2019` for more details.
+    The default value of ``fusion_multiplier`` is 1.
 
 * ``<collision_name>.fusion_probability_threshold`` (`float`) optional.
     Only for ``nuclearfusion``.
@@ -1799,25 +1800,21 @@ WarpX provides several particle collision models, using varying degrees of appro
     where :math:`\beta` is the term on the r.h.s except :math:`W_b`.
 
 * ``<collision_name>.scattering_processes`` (`strings` separated by spaces)
-    Only for ``background_mcc``. The MCC scattering processes that should be
+    Only for ``dsmc`` and ``background_mcc``. The scattering processes that should be
     included. Available options are ``elastic``, ``back`` & ``charge_exchange``
     for ions and ``elastic``, ``excitationX`` & ``ionization`` for electrons.
-    The ``elastic`` option uses hard-sphere scattering, with a differential
-    cross section that is independent of angle.
-    With ``charge_exchange``, the ion velocity is replaced with the neutral
-    velocity, chosen from a Maxwellian based on the value of
-    ``<collision_name>.background_temperature``.
     Multiple excitation events can be included for electrons corresponding to
     excitation to different levels, the ``X`` above can be changed to a unique
     identifier for each excitation process. For each scattering process specified
-    a path to a cross-section data file must  also be given. We use
+    a path to a cross-section data file must also be given. We use
     ``<scattering_process>`` as a placeholder going forward.
 
 * ``<collision_name>.<scattering_process>_cross_section`` (`string`)
-    Only for ``background_mcc``. Path to the file containing cross-section data
+    Only for ``dsmc`` and ``background_mcc``. Path to the file containing cross-section data
     for the given scattering processes. The cross-section file must have exactly
     2 columns of data, the first containing equally spaced energies in eV and the
-    second the corresponding cross-section in :math:`m^2`.
+    second the corresponding cross-section in :math:`m^2`. The energy column should
+    represent the kinetic energy of the colliding particles in the center-of-mass frame.
 
 * ``<collision_name>.<scattering_process>_energy`` (`float`)
     Only for ``background_mcc``. If the scattering process is either

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