Fusion Extension implementation - inter-species

docs/source/index.rst

@@ -86,6 +86,7 @@ In case you are already fluent in compiling C++ projects and HPC, running PIC si
    models/collisional_ionization
    models/photons
    models/binary_collisions
+   models/fusion
    models/atomic_physics
 
 .. toctree::

docs/source/models/fusion.rst

@@ -0,0 +1,208 @@
+.. _model-fusion:
+
+Fusion
+======
+
+Overview
+--------
+
+The PIConGPU Fusion Extension enables Monte Carlo simulation of nuclear fusion reactions between macro-particles based on [Wu2021]_ and [Higginson2019]_. The extension implements fully relativistic binary collision algorithms with local charge and (if possible) mass conservation, producing physically accurate fusion products with proper energy-momentum distributions.
+
+See also: the workflow guide in :doc:`../usage/workflows/fusionReactions` for setup and parameter examples.
+
+Use Cases
+---------
+- Inertial confinement fusion (ICF) plasma simulations
+- High-energy density physics with fusion heating
+- Thermonuclear burn studies in laser-plasma interactions
+
+Key Features
+------------
+- Relativistic collision algorithm and cross-section evaluation
+- Stoichiometric multiplicity with compile-time weight validation
+- Local charge conservation always; local mass conservation when non-degenerate
+- Support for multiple simultaneous fusion channels
+
+Species Definition Requirements
+-------------------------------
+For species setup details and examples, see the workflow page:
+:doc:`../usage/workflows/fusionReactions`.
+
+Fusion Reaction Configuration (Summary)
+--------------------------------------
+For a minimal configuration pattern and parameter examples, see the workflow page:
+:doc:`../usage/workflows/fusionReactions`.
+
+Pipeline configuration and parameters are detailed in the workflow page.
+
+Fusion Particle Creation Algorithm
+----------------------------------
+
+Purpose
+-------
+The fusion particle creation algorithm enforces local charge conservation by calculating compile-time weights for product particles. Local charge conservation within each computational cell is essential for numerical stability in PIC simulations, preventing spurious electric fields from charge imbalances.
+
+Problem Definition
+------------------
+Input: Two reactant particles of some weight undergoing fusion. The Fusion algorithm determines the amount of fuel (reactants) to consume in the fusion process. This fuel has an associated weight W_p that represents the total weight of products (of each species) to be created. The fractional weights W₁, W₂, W₃, W₄ are then multiplied by W_p to determine the actual weights of the outgoing particles.
+
+Stoichiometric Multiplicity Limits and Invariants
+-------------------------------------------------
+We derive compile-time multiplicity limits (c₃, c₄) from species charges (Z) and mass numbers (A). These limits define how much of each product species must be created in total (across both reactant sites).
+
+- Non-degenerate (det ≠ 0), where det := q₃m₄ − q₄m₃:
+
+  .. math::
+
+     \begin{bmatrix} q_3 & q_4 \\ m_3 & m_4 \end{bmatrix}
+     \begin{bmatrix} c_3 \\ c_4 \end{bmatrix}
+     =
+     \begin{bmatrix} q_1{+}q_2 \\ m_1{+}m_2 \end{bmatrix}
+
+- Degenerate (q₃/m₃ = q₄/m₄): symmetric limits
+
+  .. math::
+
+     c_3 = c_4 = \frac{q_1{+}q_2}{q_3{+}q_4}\quad (q_3{+}q_4 \ne 0),\qquad
+     c_3 = c_4 = \frac{m_1{+}m_2}{m_3{+}m_4}\quad (\text{if } q_3{+}q_4 \approx 0)
+
+With multiplicity limits, the invariants are:
+  - Weight bounds: 0 ≤ Wᵢ ≤ cᵢ
+  - Weight conservation: W₁ + W₃ = c₃, W₂ + W₄ = c₄
+  - Local charge: q₁ = W₁q₃ + W₂q₄ and q₂ = W₃q₃ + W₄q₄
+  - Local mass (when Algorithm 1 is valid): m₁ = W₁m₃ + W₂m₄ and m₂ = W₃m₃ + W₄m₄
+
+Algorithm 1: Mass-Charge Conservation under Multiplicity Limits
+----------------------------------------------------------------
+Solve for site-1 weights (W₁, W₂) using:
+
+.. math::
+
+   \begin{bmatrix} q_3 & q_4 \\ m_3 & m_4 \end{bmatrix}
+   \begin{bmatrix} W_1 \\ W_2 \end{bmatrix} = \begin{bmatrix} q_1 \\ m_1 \end{bmatrix}
+
+Then set W₃ = c₃ − W₁ and W₄ = c₄ − W₂. Accept only if det ≠ 0, all 0 ≤ Wᵢ ≤ cᵢ, and site-2 mass/charge match.
+
+Algorithm 2: Charge-Only Conservation under Multiplicity Limits
+---------------------------------------------------------------
+If Algorithm 1 is invalid (det ≈ 0 or out-of-limit weights), enforce only local charge with multiplicity limits. Case handling is neutral-friendly and macroparticle-minimizing.
+
+Implementation Details
+---------------------
+In  ``include/picongpu/particles/fusion/detail/Creation.hpp`` both algorithms are evaluated at compile time using constexpr functions:
+
+- ``computeStoichiometryCaps()``: derives (c₃, c₄) from species
+- ``calculateMassChargeConservingWeightsWithCaps()``: local mass+charge under multiplicity limits
+- ``calculateChargeOnlyWithCaps()``: charge-only robust fallback under multiplicity limits
+
+Relativistic Kinematics and Fusion Sampling (Cannoni, 2016)
+-----------------------------------------------------------
+For each candidate pair we compute kinematics using Lorentz-invariant quantities to avoid loss of precision at high γ, following Cannoni’s formulation of invariant relative velocity [Cannoni2016]_. We adopt the metric signature (+, −, −, −) and four-momentum convention :math:`p_i^\mu = (E_i/c, \vec{p}_i)`.
+
+Key kinematic invariants (computed from lab-frame inputs)
+
+- Total four-momentum: :math:`P^\mu = p_1^\mu + p_2^\mu = (E_{\mathrm{tot}}/c, \vec{p}_{\mathrm{tot}})`.
+- Mandelstam s (energy-squared invariant):
+
+  .. math::
+
+     s \;:=\; c^2 P^\mu P_\mu 
+       \,=\, E_{\mathrm{tot}}^2 - |\vec{p}_{\mathrm{tot}}|^2 c^2 \,=\, (p_1^\mu + p_2^\mu)^2 c^2.
+
+- Center-of-mass (CM) energy and velocity:
+
+  .. math::
+
+     E_{\mathrm{cm}} = \sqrt{s},\quad
+     \vec{V}_{\mathrm{cm}} = \frac{c^2\,\vec{p}_{\mathrm{tot}}}{E_{\mathrm{tot}}},\quad
+     \gamma_{\mathrm{cm}} = \left(1- \frac{|\vec{V}_{\mathrm{cm}}|^2}{c^2}\right)^{-1/2}.
+
+- Relative Lorentz factor (Cannoni):
+
+  .. math::
+
+     \gamma_r \,=\, \frac{s - m_1^2 c^4 - m_2^2 c^4}{2 m_1 m_2 c^4}, \qquad
+     |v_{\mathrm{rel}}| \,=\, c\,\sqrt{1 - \gamma_r^{-2}}.
+
+Associated Møller (invariant) relative speed (used for rates):
+
+.. math::
+
+   v_{M} \,=\, |v_{\mathrm{rel}}|\,\gamma_{\mathrm{cm}}.
+
+Fusion probability and sampling
+-------------------------------
+- The microscopic cross section functor :math:`\sigma(E_{\mathrm{cm}})` (e.g., Bosch–Hale) is evaluated at the pair’s CM energy :math:`\sqrt{s}`.
+- Over a timestep :math:`\Delta t`, the acceptance probability scales as
+
+  .. math::
+
+     P \;\propto\; \sigma(E_{\mathrm{cm}})\, v_M\, \Delta t
+
+- On acceptance, sample an isotropic direction in the CM frame; compute product total energies
+
+  .. math::
+
+     E_{3,\mathrm{cm}} = \frac{s + (m_3^2 - m_4^2)c^4}{2\,\sqrt{s}}, \qquad
+     E_{4,\mathrm{cm}} = E_{\mathrm{cm}} - E_{3,\mathrm{cm}},
+
+and the common CM momentum magnitude from :math:`E^2 = p^2 c^2 + m^2 c^4`. Notation (CM frame). Let :math:`\vec{P}_1` denote the three-momentum of product 1 and :math:`\vec{P}_2` that of product 2, then
+
+.. math::
+
+  \vec{P}_1 = -\vec{P}_2,\qquad |\vec{P}_1| = |\vec{P}_2| = \frac{\sqrt{E_{3,\mathrm{cm}}^2 - m_3^2 c^4}}{c} = \frac{\sqrt{E_{4,\mathrm{cm}}^2 - m_4^2 c^4}}{c}.
+
+Finally, boost both product four-momenta back to the lab frame with :math:`\vec{V}_{\mathrm{cm}}`.
+
+Code Flow and Implementation
+----------------------------
+
+Fusion Extension Architecture
+-----------------------------
+
+The fusion extension follows this execution flow::
+
+    simulation.hpp → Fusion.x.cpp → Collider.hpp → WithPeer.hpp
+                                                                 └─ IntraCollision.hpp 
+
+Core Components
+---------------
+- ``Inter/Intra-Collision.hpp``: Combines collision algorithm from the Collision Extension with the Creation Kernel
+- ``FusionFunctor.hpp``: Interfaces between collision framework and physics algorithms
+- ``FusionAlgorithm.hpp``: Implements relativistic fusion physics and product momentum calculation
+
+Algorithm Chain
+---------------
+
+.. code-block:: text
+
+                        uses: (particles/fusion/detail)
+    InterCollision.hpp ──────────────────────────────→ FusionFunctor.hpp → FusionAlgorithm.hpp
+                     └─ Creation.hpp
+
+References
+----------
+
+.. [Wu2021]
+        D. Wu, Z. M. Sheng, W. Yu, S. Fritzsche, and X. T. He.
+        *A pairwise nuclear fusion algorithm for particle-in-cell simulations: Weighted particles at relativistic energies.*
+        AIP Advances 11, 075003 (2021).
+        https://doi.org/10.1063/5.0051178
+
+.. [Higginson2019]
+        D. P. Higginson, A. Link, and A. Schmidt.
+        *A pairwise nuclear fusion algorithm for weighted particle-in-cell plasma simulations.*
+        Journal of Computational Physics 388, 439–453 (2019).
+        https://doi.org/10.1016/j.jcp.2019.03.020
+
+.. [Cannoni2016]
+        M. Cannoni.
+        *Lorentz invariant relative velocity and relativistic binary collisions.*
+        arXiv:1605.00569 [hep-ph] (2016).
+        https://arxiv.org/abs/1605.00569v2
+
+.. [Takizuka1977]
+        T. Takizuka and H. Abe.
+        *A binary collision model for plasma simulation with a particle code.*
+        Journal of Computational Physics 25(3), 205–219 (1977).
+        https://doi.org/10.1016/0021-9991(77)90099-7

docs/source/usage/workflows.rst

@@ -18,5 +18,6 @@ This section contains typical user workflows and best practices.
    workflows/probeParticles
    workflows/tracerParticles
    workflows/particleFilters
+   workflows/fusionReactions.rst
    workflows/synchrotronExtensionDocumentation
    workflows/ionization

docs/source/usage/workflows/fusionReactions.rst

@@ -0,0 +1,147 @@
+.. _fusionReactions:
+
+Fusion Workflow
+===============
+
+This page explains how to configure and run the Fusion Extension in simulations.
+For the physics background, algorithms, and references, see the model page:
+:doc:`../../models/fusion`.
+
+Species Definition Requirements
+-------------------------------
+All species participating in fusion reactions must be defined with ``massRatio`` and ``chargeRatio`` flags. These are mandatory for the fusion algorithms to function correctly.
+
+To simplify defining particle masses relative to the electron mass, you can define an atomic mass unit (``amu``) constant. This constant is the ratio of one atomic mass unit to the electron rest mass:
+
+.. literalinclude:: ../../../../share/picongpu/tests/Fusion/include/picongpu/param/speciesDefinition.param
+   :language: cpp
+   :start-after: doc-include-start: amu-definition
+   :end-before: doc-include-end: amu-definition
+   :dedent:
+
+Example species definition for deuterons:
+
+.. literalinclude:: ../../../../share/picongpu/tests/Fusion/include/picongpu/param/speciesDefinition.param
+   :language: cpp
+   :start-after: doc-include-start: deuteron-definition
+   :end-before: doc-include-end: deuteron-definition
+   :dedent:
+
+Important: Mass ratios must be as precise as possible since fusion energy release is calculated from the mass deficit between reactants and products.
+
+Fusion Reaction Configuration (Summary)
+--------------------------------------
+Configure fusion reactions in ``include/picongpu/param/fusion.param``.
+
+1) Define reaction parameters
+
+.. code-block:: cpp
+
+    struct DT {
+        // Reactant species
+        using reactants = pmacc::mp_list<Pair<PIC_Tritons, PIC_Deuterons>>;
+        // Filter for reactants
+        using FilterPair = OneFilter<filter::All>;
+
+        // Product species (multiplicities deduced automatically)
+        using products = pmacc::mp_list<Pair<PIC_Neutrons, PIC_He4>>;
+
+        // Cross-section parameters (Bosch–Hale parameterization)
+        struct Params {
+            static constexpr float_X BG = 34.3827_X;
+            static constexpr float_X A1 = 6.927e4_X;
+            static constexpr float_X A2 = 7.454e8_X;
+            // ... additional parameters
+        };
+
+        using CrossSectionInterpolator = relativistic::FusionFunctor<Params>;
+    };
+
+.. note::
+   Product multiplicities (stoichiometric coefficients) are deduced automatically
+   from species charge and mass numbers. You only declare the product species;
+   multiplicity caps are inferred at compile time.
+
+   Examples:
+
+   - T + T → ⁴He + 2n:
+
+     .. code-block:: cpp
+
+         using products = pmacc::mp_list<Pair<PIC_Neutrons, PIC_He4>>;
+
+   - B + p → 3×⁴He:
+
+     .. code-block:: cpp
+
+         using products = pmacc::mp_list<Pair<PIC_He4, PIC_He4>>;
+
+Complete Example: D-T Fusion Reaction
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Here is a complete example for the Deuterium-Tritium fusion reaction (D + T → n + He⁴).
+This is the easiest fusion reaction to achieve and is commonly used in fusion energy research.
+
+The cross-section parameters are based on the Bosch-Hale parameterization from experimental data
+(Source: https://hdl.handle.net/11858/00-001M-0000-0027-6535-1, page 29).
+
+.. literalinclude:: ../../../../share/picongpu/tests/Fusion/include/picongpu/param/fusion.param
+   :language: cpp
+   :start-after: doc-include-start: DT-reaction
+   :end-before: doc-include-end: DT-reaction
+   :dedent:
+
+2) Configure the fusion pipeline
+
+.. literalinclude:: ../../../../share/picongpu/tests/Fusion/include/picongpu/param/fusion.param
+   :language: cpp
+   :start-after: doc-include-start: fusion-pipeline
+   :end-before: doc-include-end: fusion-pipeline
+   :dedent:
+
+.. note:: Multiple reactions
+
+   You can add multiple reactions to the pipeline by defining additional reaction structs
+   and adding them to the ``FusionPipeline`` list. The reactions are processed sequentially
+   in the order they appear.
+
+Simulation Parameters
+---------------------
+
+.. literalinclude:: ../../../../share/picongpu/tests/Fusion/include/picongpu/param/fusion.param
+   :language: cpp
+   :start-after: doc-include-start: simulation-parameters
+   :end-before: doc-include-end: simulation-parameters
+   :dedent:
+
+.. note:: Fusion production multiplier (Fmult)
+
+   The ``Fmult`` parameter controls how the consumed reactant weight (the
+   minimum of the two reactant weights in a pair) is split into multiple product
+   macro-particles while keeping total weight conserved. Increasing ``maxFmult`` 
+   allows more, lighter products to be created. Decreasing ``maxFmult`` results 
+   in fewer, heavier products.
+
+   Algorithm (per fused pair):
+
+   1. Start with ``Fmult = maxFmult``
+   2. Calculate the base fusion probability ``P``
+   3. Adjust ``Fmult`` to ensure valid Monte Carlo sampling:
+
+      - If ``P > 1.0``: Set ``Fmult = 1.0`` and ``P = 1.0`` (a warning is issued)
+      - If ``0 < P ≤ 1.0``: Calculate the maximum allowed multiplier as 
+        ``Fmult = min(maxFmult, 0.99/P)`` to ensure ``P * Fmult ≤ 0.99``
+
+   4. Update the fusion probability: ``P = P * Fmult``
+   5. Compute ``productWeighting = minWeighting / Fmult``
+   6. Distribute the ``productWeighting`` across the two reactant sites and 
+      product species using the site fractions ``W₁…W₄``
+
+   This approach ensures that the fusion probability remains below 1.0 (required
+   for proper Monte Carlo sampling) while maximizing the number of product 
+   particles created per fusion event.
+
+See also
+--------
+- Physics and algorithms: :doc:`../../models/fusion`
+- Related: :doc:`../../models/binary_collisions`