Skip to content

Commit

Permalink
WIP update collision documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
roelof-groenewald committed Sep 3, 2023
1 parent 20e256e commit d90d7c2
Show file tree
Hide file tree
Showing 7 changed files with 128 additions and 53 deletions.
28 changes: 28 additions & 0 deletions Docs/source/refs.bib
Original file line number Diff line number Diff line change
Expand Up @@ -483,3 +483,31 @@ @article{Stanier2020
author = {A. Stanier and L. Chacón and A. Le},
keywords = {Hybrid, Particle-in-cell, Plasma, Asymptotic-preserving, Cancellation problem, Space weather},
}

@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}
}
73 changes: 61 additions & 12 deletions Docs/source/theory/collisions.rst
Original file line number Diff line number Diff line change
Expand Up @@ -3,18 +3,36 @@
Collisions
==========

Monte Carlo 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).

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:`ScatteringProcess` 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.
.. _theory-collisions-mcc:

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.
Monte Carlo Collisions
----------------------

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
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::
Expand All @@ -23,19 +41,44 @@ 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
-----------------------------

(Under construction)

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 hard-sphere scattering, 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}
Expand Down Expand Up @@ -74,6 +117,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 bencmark of the MCC
implementation against literature results.

Particle cooling due to elastic collisions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Expand Down
2 changes: 2 additions & 0 deletions Docs/source/usage/examples.rst
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,8 @@ Uniform plasma
:download:`2D case <../../../Examples/Physics_applications/uniform_plasma/inputs_2d>`
:download:`3D case <../../../Examples/Physics_applications/uniform_plasma/inputs_3d>`

.. _examples-mcc-turner:

Capacitive discharge
--------------------

Expand Down
34 changes: 15 additions & 19 deletions Docs/source/usage/parameters.rst
Original file line number Diff line number Diff line change
Expand Up @@ -1565,6 +1565,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.
Expand All @@ -1574,29 +1575,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.
Expand All @@ -1619,7 +1620,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``.
Expand All @@ -1630,8 +1631,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_multipler``.
(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``.
Expand Down Expand Up @@ -1701,22 +1702,17 @@ 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`.
Expand Down
38 changes: 19 additions & 19 deletions Regression/Checksum/benchmarks_json/Python_dsmc_1d.json
Original file line number Diff line number Diff line change
@@ -1,27 +1,27 @@
{
"lev=0": {
"rho_electrons": 0.0044367956712501894,
"rho_he_ions": 0.005200487707892467
"rho_electrons": 0.004436602398896733,
"rho_he_ions": 0.0052003262664415285
},
"he_ions": {
"particle_momentum_x": 2.7735512966774165e-19,
"particle_momentum_y": 2.7574111491186894e-19,
"particle_momentum_z": 3.620520352986572e-19,
"particle_position_x": 2201.236370518716,
"particle_weight": 17190734375000.002
},
"electrons": {
"particle_momentum_x": 3.5197641187987634e-20,
"particle_momentum_y": 3.5395348457892096e-20,
"particle_momentum_z": 1.2572070275928055e-19,
"particle_position_x": 2139.282134158907,
"particle_weight": 14584539062500.002
"particle_momentum_x": 3.50212700099208e-20,
"particle_momentum_y": 3.5368926859820716e-20,
"particle_momentum_z": 1.2588108956625115e-19,
"particle_position_x": 2139.6498177543617,
"particle_weight": 14582968750000.002
},
"neutrals": {
"particle_momentum_x": 1.4062123597534044e-19,
"particle_momentum_y": 1.4053058130860795e-19,
"particle_momentum_z": 1.4094729915175569e-19,
"particle_position_x": 1122.8200012433333,
"particle_weight": 6.458799999999999e+19
},
"he_ions": {
"particle_momentum_x": 2.772918878614073e-19,
"particle_momentum_y": 2.7576133543773225e-19,
"particle_momentum_z": 3.619033046968002e-19,
"particle_position_x": 2201.5707517530504,
"particle_weight": 17191519531250.002
"particle_momentum_x": 1.405588503355727e-19,
"particle_momentum_y": 1.408077689882847e-19,
"particle_momentum_z": 1.4024616940779626e-19,
"particle_position_x": 1121.2330379095083,
"particle_weight": 6.4588e+19
}
}
4 changes: 2 additions & 2 deletions Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ public:
*
* @param[in] collision_name the name of the collision
*/
DSMC (const std::string collision_name);
DSMC (std::string collision_name);

/** Perform the collisions
*
Expand All @@ -75,7 +75,7 @@ public:
*
*/
void doCollisionsWithinTile (
amrex::Real dt, int const lev, amrex::MFIter const& mfi,
amrex::Real dt, int lev, amrex::MFIter const& mfi,
WarpXParticleContainer& species_1,
WarpXParticleContainer& species_2,
SmartCopy& copy_species1,
Expand Down
2 changes: 1 addition & 1 deletion Source/Particles/Collision/BinaryCollision/DSMC/DSMC.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ DSMC::doCollisions (amrex::Real /*cur_time*/, amrex::Real dt, MultiParticleConta

// Enable tiling
amrex::MFItInfo info;
if (amrex::Gpu::notInLaunchRegion()) info.EnableTiling(species1.tile_size);
if (amrex::Gpu::notInLaunchRegion()) info.EnableTiling(WarpXParticleContainer::tile_size);

// Loop over refinement levels
for (int lev = 0; lev <= species1.finestLevel(); ++lev){
Expand Down

0 comments on commit d90d7c2

Please sign in to comment.