Merge branch 'ComputationalRadiationPhysics:dev' into Fusion_interSpecies

Author: FilipO28555

docs/source/postprocessing/python.rst

@@ -3,7 +3,10 @@
 Python
 ======
 
-.. sectionauthor:: Axel Huebl, Klaus Steiniger
+.. sectionauthor:: Axel Huebl, Klaus Steiniger, Richard Pausch
+
+If you are familiar with Python and would like to install the Python libraries that come with PIConGPU, please refer to the install instructions in our :doc:`python utilities section <../usage/python_utils>`.
+
 
 If you are new to python, get your hands on the tutorials of the following important libraries to get started.
 
@@ -17,7 +20,7 @@ https://mamba.readthedocs.io/en/latest/installation/micromamba-installation.html
 
 and set up a Micromamba environment with the help of this Micromamba environment file
 
-https://gist.github.com/steindev/d19263d41b0964bcecdfb1f47e18a86e
+https://gist.github.com/steindev/c1d678375fd0f974f4b9cb1480b3da03
 
 (see documentation within the file).
 

docs/source/usage/plugins/radiation.rst

@@ -301,8 +301,8 @@ Command line option                                      Description
                                                          Default: ``{}`` (no JSON/TOML configuration used)
 ``--<species>_radiation.distributedAmplitude``           Activate the optional output of distributed amplitudes per MPI rank.  in the openPMD output.
                                                          Default: ``0`` (deactivated/no additional output)
-					  
-========================================= ==============================================================================================================================
+
+======================================================== ==============================================================================================================================
 
 Memory Complexity
 ^^^^^^^^^^^^^^^^^

docs/source/usage/workflows/fusionReactions.rst

@@ -9,15 +9,15 @@ 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.
 
 **Use Cases:**
-- Inertial confinement fusion (ICF) plasma simulations
-- High-energy density physics with fusion heating
-- Thermonuclear burn studies in laser-plasma interactions
+  - 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
+  - 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
 
 How to Use the Fusion Extension
 ===============================
@@ -90,7 +90,7 @@ Configure fusion reactions in ``include/picongpu/param/fusion.param``:
 
    Examples:
 
-   - T + T → n + ⁴He (also covers T + T → ⁴He + 2n): it suffices to write
+   - T + T → ⁴He + 2n: it suffices to write
 
      .. code-block:: cpp
 
@@ -154,23 +154,21 @@ Configure fusion reactions in ``include/picongpu/param/fusion.param``:
 
 .. note:: Fusion production multiplier (Fmult)
 
+   Increasing ``maxFmult`` allows more, lighter products to be created, down to the ``productMinWeighting`` threshold. Decreasing ``maxFmult`` results in fewer, heavier products.
    ``Fmult`` 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.
 
    Algorithm (per fused pair):
 
    - Start with ``Fmult = maxFmult``
    - Compute ``productWeighting = minWeighting / Fmult``
    - If ``productWeighting < productMinWeighting``, reduce ``Fmult`` to ``max(1, minWeighting / productMinWeighting)`` and recompute ``productWeighting``
-   - Guarantees: ``productWeighting ≥ productMinWeighting`` and ``Fmult ≥ 1`` while conserving the total produced weight ``minWeighting``
 
-   Increasing ``maxFmult`` allows more, lighter products to be created, down to the ``productMinWeighting`` threshold. Decreasing ``maxFmult`` results in fewer, heavier products.
-   
    The base ``productWeighting`` is then distributed across the two reactant sites and product species using the site fractions ``W₁…W₄`` (see Particle Creation Algorithm). In symmetric cases this often results in roughly half of a species' weight at each reactant site.
 
 **Parameter Guidelines:**
-- ``productMinWeighting``: Minimum allowed weight per two created product macro-particles; 16 let's the particle split in half 4 times before hitting the threshold. (see: ''Fusion Particle Creation Algorithm'')
-- ``maxFmult``: Upper bound for how many product macroparticles can be spawned per fused pair; higher values allow more (lighter) products down to ``productMinWeighting``.
-- ``cellListChunkSize``: Use ``TYPICAL_PARTICLES_PER_CELL`` for optimal memory usage.
+  - ``productMinWeighting``: Minimum allowed weight per two created product macro-particles; 16 let's the particle split in half 4 times before beeing less than 1. (see: ''Fusion Particle Creation Algorithm'')
+  - ``maxFmult``: Upper bound for how many product macroparticles can be spawned per fused pair; higher values allow more (lighter) products down to ``productMinWeighting``. Usually 10~100 times less than typical particle weighting.
+  - ``cellListChunkSize``: Use ``TYPICAL_PARTICLES_PER_CELL`` for optimal memory usage.
 
 Fusion Particle Creation Algorithm
 ==================================
@@ -188,19 +186,19 @@ Problem Definition
 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.
 
 **Parameter Definitions:**
-- q₁ = charge of reactant 1
-- q₂ = charge of reactant 2  
-- q₃ = charge of product species 1
-- q₄ = charge of product species 2
-- m₁ = mass of reactant 1
-- m₂ = mass of reactant 2
-- m₃ = mass of product species 1
-- m₄ = mass of product species 2
-- W_p = total weight of products (fuel consumption amount)
+  - q₁ = charge of reactant 1
+  - q₂ = charge of reactant 2  
+  - q₃ = charge of product species 1
+  - q₄ = charge of product species 2
+  - m₁ = mass of reactant 1
+  - m₂ = mass of reactant 2
+  - m₃ = mass of product species 1
+  - m₄ = mass of product species 2
+  - W_p = total weight of products (fuel consumption amount)
 
 **Output:** Four product particles with final weights W₁×W_p, W₂×W_p, W₃×W_p, W₄×W_p
-- Product species 1: weights W₁×W_p at reactant 1 position, W₃×W_p at reactant 2 position
-- Product species 2: weights W₂×W_p at reactant 1 position, W₄×W_p at reactant 2 position
+  - Product species 1: weights W₁×W_p at reactant 1 position, W₃×W_p at reactant 2 position
+  - Product species 2: weights W₂×W_p at reactant 1 position, W₄×W_p at reactant 2 position
 
 Stoichiometric Multiplicity Limits and Invariants
 ----------------------------------
@@ -224,10 +222,10 @@ We derive compile-time multiplicity limits (c₃, c₄) from species charges (Z)
      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ᵢ (i ∈ {product 1, product 2} per site)
-- 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₄
+  - Weight bounds: 0 ≤ Wᵢ ≤ cᵢ (i ∈ {product 1, product 2} per site)
+  - 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 Implementation
 -----------------------------
@@ -261,7 +259,7 @@ Algorithm 1 is considered failed when:
 
 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
@@ -278,7 +276,7 @@ Fusion Extension Architecture
 
 The fusion extension follows this execution flow::
 
-    simulation.hpp → Fusion.x.cpp → Collider.hpp → WithPeer.hpp → InterCollision.hpp
+    simulation.hpp → Fusion.x.cpp → Collider.hpp → WithPeer.hpp
                                                                  └─ IntraCollision.hpp 
 
 **Core Components:**
@@ -298,139 +296,99 @@ The fusion extension follows this execution flow::
 Main Physical Algorithms
 ========================
 
-1. Fusion Process - Binary Particle Collisions
------------------------------------------------
+Algorithmic Overview (per cell)
+--------------------------------
 
-**Physical Model:** Two reactant macro-particles undergo fusion, producing up to 4 product macro-particles
-- Reactant selection within computational cells
-- Energy-dependent fusion probability: σ(E) using empirical fits
-- Statistical weighting to handle macro-particle representation
+.. code-block:: text
 
-2. Particle Creation Algorithm
-------------------------------
+    1) Randomly pair reactants in the cell (Takizuka–Abe)
+    2) Compute relativistic kinematics and sample fusion (Cannoni)
+    3) Create product macroparticles with local conservation (new)
 
-**Physical Model:** Local charge (and when possible local mass) conservation during particle creation under stoichiometric caps
+Each stage is summarized below with references and links to implementation.
 
-**Algorithm 1 - Mass-Charge (under caps):**
-2×2 local solve for W₁,W₂; set W₃,W₄ from caps; accept if caps and site-2 checks pass
+1. Binary Pairing of Reactants (Takizuka–Abe, 1977)
+---------------------------------------------------
 
-**Algorithm 2 - Charge-Only (under caps):**
-Neutral-aware, macro-minimizing case analysis; always respects caps and local charge
+Physical idea
+  Within each cell we form unbiased binary encounters by random pairing of two reactant lists following Takizuka–Abe [Takizuka1977]_. If the lists have different sizes, particles from the shorter list are deterministically reused so that all particles participate; see :ref:`model-binaryCollisions::details:duplication` for details. For a broader description of the pairing algorithm, see :ref:`model-binaryCollisions`.
 
-**Key Physics:**
-- Local charge conservation: q₁ = W₁q₃ + W₂q₄, q₂ = W₃q₃ + W₄q₄
-- Cap constraints: 0 ≤ Wᵢ ≤ cᵢ, W₁ + W₃ = c₃, W₂ + W₄ = c₄
-- Local mass conservation when Algorithm 1 is valid
+Algorithm (per cell)
+  - Build two filtered lists A and B of reactant indices (intra-collision: A=B; inter-collision: A≠B).
+  - Apply random permutations to the longer list.
+  - Pair i-th entries up to N_p = min(|A|, |B|): pairs (A[i], B[i]).
+  - If |A| ≠ |B|, handle duplicates of the shorter list so each long-list entry has a partner (The weight of a particle from the shorter list is divided; see implementation note below).
 
-3. Relativistic Framework
--------------------------
+Implementation
+  - Shuffling and pairing are performed by the collision kernels: ``particles/fusion/InterCollision.hpp`` and ``IntraCollision.hpp``.
+  - Deterministic duplication for |A|≠|B|: ``duplicationCorrection()``; randomized ordering via ``detail::shuffle``.
+  - Number-density weighting used later is computed per cell to keep correct rates under unequal weights and counts.
 
-**Physical Model:** Fully relativistic treatment using Lorentz invariants
+2. Relativistic Kinematics and Fusion Sampling (Cannoni, 2016)
+--------------------------------------------------------------
 
-**Implementation in ``FusionAlgorithm.hpp``:**
+Physical model
+  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)`.
 
-**Mandelstam Variable Approach:**
-[Cannoni2016]_
-- Uses Lorentz invariant s = E²_total - (pc)² instead of double transformations
-- Avoids numerical errors in highly relativistic regime
-- Calculates relative velocity through invariant: γᵣ = (s - m₁²c⁴ - m₂²c⁴)/(2m₁m₂c⁴)
+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):
 
-**Lorentz Transformations:**
-- Lab frame → Center-of-mass frame for isotropic product generation
-- CM frame → Lab frame for final product momenta
-- Boost formulas: **p**_lab = **p**_cm + [(**V**_cm·**p**_cm)γ_cm/(γ_cm+1) + γ_cm E_cm/c]**V**_cm/c
+    .. math::
 
-**Key Physics:**
-- Energy-momentum conservation in all reference frames
-- Relativistic energy: E = √[(pc)² + (mc²)²]
-- Isotropic angular distribution in CM frame
+       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.
 
-Common Fusion Reaction Examples  
---------------------------------
+  - 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
+
+    where the factors account for unequal list sizes, macro-weights, and multiplicity splitting and are computed in the collision kernel.
+  - 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}}`.
+
+Implementation
+  - Accumulators and boosts: ``particles/fusion/relativistic/FusionAlgorithm.hpp`` (``acc::Variables`` and ``FusionAlg::fuse``).
+  - Unit conversions (e.g., eV/keV, millibarn→area) are handled internally before σ and P are evaluated.
 
-.. list-table:: Weight Assignments for Common Fusion Reactions (under multiplicity limits)
-   :header-rows: 1
-   :widths: 25 20 10 10 10 10
-
-   * - Reaction
-     - Charges (q₁,q₂,q₃,q₄)
-     - W₁
-     - W₂  
-     - W₃
-     - W₄
-   * - D + T → ⁴He + n
-     - (1,1,2,0)
-     - 0.5
-     - 0
-     - 0.5
-     - 1
-   * - D + D → T + p
-     - (1,1,1,1)
-     - 0.5
-     - 0.5
-     - 0.5
-     - 0.5
-   * - D + D → ³He + n
-     - (1,1,2,0)
-     - 0.5
-     - 0.5
-     - 0.5
-     - 0.5
-   * - ³He + D → ⁴He + p
-     - (2,1,2,1)
-     - 0.5
-     - 1
-     - 0.5
-     - 0
-   * - ³He + T → ⁴He + D
-     - (2,1,2,1)
-     - 1
-     - 0
-     - 0
-     - 1
-   * - T + T → ⁴He + 2n
-     - (1,1,2,0)
-     - 0.5
-     - 1
-     - 0.5
-     - 1
-   * - ³He + ³He → ⁴He + 2p
-     - (2,2,2,1)
-     - 0.5
-     - 1
-     - 0.5
-     - 1
-   * - p + ⁶Li → ³He + ⁴He
-     - (1,3,2,2)
-     - 0.5
-     - 0
-     - 0.5
-     - 1
-   * - D + ⁶Li → ⁴He + ⁴He
-     - (1,3,2,2)
-     - 0.5
-     - 0
-     - 0.5
-     - 1
-   * - p + ⁷Li → ⁴He + ⁴He
-     - (1,3,2,2)
-     - 0.5
-     - 0
-     - 0.5
-     - 1
-   * - p + ¹¹B → 3×⁴He
-     - (1,5,2,2)
-     - 0.5
-     - 0
-     - 1
-     - 1.5
-
-**Weight Interpretation:**
-- W₁, W₃: weights of product 1 at reactant positions 1, 2
-- W₂, W₄: weights of product 2 at reactant positions 1, 2
-- Multiplicity limits can be fractional and > 1; weights may exceed 1 accordingly
-- If Algorithm 1 is valid, neutrals can split across sites (mass+charge enforced)
-- Otherwise, Algorithm 2 minimizes the number of neutral macroparticles (all neutrals at one site)
 
 References
 ----------
@@ -452,3 +410,9 @@ References
         *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

include/picongpu/particles/fusion/relativistic/FusionAlgorithm.hpp

@@ -170,8 +170,7 @@ namespace picongpu::particles::fusion::relativistic
                 }
 
                 // Relativistic formula for the TOTAL energy of product 0 in the CM frame
-                float_COLL const E_p0_tot_cm
-                    = (E_cm_tot * E_cm_tot + (mP0 * mP0 - mP1 * mP1) * c4) / (2.0_COLL * E_cm_tot);
+                float_COLL const E_p0_tot_cm = (E_cm_tot_sq + (mP0 * mP0 - mP1 * mP1) * c4) / (2.0_COLL * E_cm_tot);
                 // TOTAL energy of product 1
                 float_COLL const E_p1_tot_cm = E_cm_tot - E_p0_tot_cm;