Added CI test and updated documentation to referance fusion.param file in it

Author: FilipO28555

docs/source/usage/workflows/fusionReactions.rst

@@ -11,28 +11,21 @@ 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, PIConGPU provides a constant for the atomic mass unit (``amu``). This constant is the ratio of one atomic mass unit to the electron rest mass. You can use it in ``speciesDefinition.param`` as follows:
+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:
 
-.. code-block:: cpp
-
-    #include "picongpu/unitless/simulation.unitless"
-
-    // Atomic mass unit in electron masses (amu/me)
-    constexpr float_X amu = 1.0_X / sim.fusion.electronMassAMU;
-
-    // Deuteron species
-    value_identifier(float_X, MassRatioDeuterons, 2.013553212745*amu);
-    value_identifier(float_X, ChargeRatioDeuterons, 1.0);
+.. 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:
 
-    using ParticleFlagsDeuterons = MakeSeq_t<
-        particlePusher<UsedParticlePusher>,
-        shape<UsedParticleShape>,
-        interpolation<UsedField2Particle>,
-        massRatio<MassRatioDeuterons>,        // Required for fusion
-        chargeRatio<ChargeRatioDeuterons>     // Required for fusion
-    >;
+Example species definition for deuterons:
 
-    using PIC_Deuterons = Particles<PMACC_CSTRING("d"), ParticleFlagsDeuterons, DefaultParticleAttributes>;
+.. 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.
 
@@ -92,92 +85,34 @@ This is the easiest fusion reaction to achieve and is commonly used in fusion en
 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).
 
-.. code-block:: cpp
-
-    struct DT
-    {
-        //! Reactant species: Deuterium and Tritium
-        using reactants = pmacc::mp_list<Pair<PIC_Tritons, PIC_Deuterons>>;
-        
-        //! Filter for reactants (mandatory for ColliderFromStruct)
-        //! Use filter::All to include all particles, or define custom filters
-        using FilterPair = OneFilter<filter::All>;
-        
-        //! Product species: Neutron and Helium-4
-        using products = pmacc::mp_list<Pair<PIC_Neutrons, PIC_He4>>;
-        
-        /** Cross-section parameters for D-T fusion */
-        struct Params
-        {
-            //! Gamow constant: BG = π·α·Z₁·Z₂·√(2·μ·c²) in keV^(1/2)
-            //! where α is the fine structure constant, Z₁,Z₂ are charges, μ is reduced mass
-            static constexpr float_X BG = 34.3827_X;
-            
-            //! Bosch-Hale parameterization coefficients for D-T cross-section
-            //! These coefficients fit experimental cross-section data
-            static constexpr float_X A1 = 6.927e4_X;
-            static constexpr float_X A2 = 7.454e8_X;
-            static constexpr float_X A3 = 2.050e6_X;
-            static constexpr float_X A4 = 5.2002e4_X;
-            static constexpr float_X A5 = 0.0_X;
-            
-            //! Additional Bosch-Hale coefficients
-            static constexpr float_X B1 = 6.38e1_X;
-            static constexpr float_X B2 = -9.95e-1_X;
-            static constexpr float_X B3 = 6.981e-5_X;
-            static constexpr float_X B4 = 1.728e-4_X;
-        };
-        
-        //! Cross-section calculation using Bosch-Hale parameterization
-        //! Returns cross-section in milli-barns
-        using CrossSectionInterpolator = relativistic::FusionFunctor<Params>;
-    };
+.. 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
 
-.. code-block:: cpp
-
-    // Single reaction
-    using FusionPipeline = pmacc::mp_list<ColliderFromStruct<DT>>;
-    
-    // Multiple reactions (processed sequentially)
-    using FusionPipeline = pmacc::mp_list<
-        ColliderFromStruct<DT>,
-        ColliderFromStruct<DD_branch1>,
-        ColliderFromStruct<DD_branch2>
-    >;
-
-.. note:: Alternative explicit form
+.. 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:
 
-   You can also define the fusion pipeline using the explicit ``Collider`` template:
+.. note:: Multiple reactions
 
-   .. code-block:: cpp
-
-      using FusionPipeline = pmacc::mp_list<
-          Collider<DT_He4n::CrossSectionInterpolator, DT_He4n::reactants, DT_He4n::products, OneFilter<filter::All>>,
-          Collider<DD_Tp::CrossSectionInterpolator, DD_Tp::reactants, DD_Tp::products, OneFilter<filter::All>>,
-          Collider<DD_He3n::CrossSectionInterpolator, DD_He3n::reactants, DD_He3n::products, OneFilter<filter::All>>,
-          Collider<He3D_He4p::CrossSectionInterpolator, He3D_He4p::reactants, He3D_He4p::products, OneFilter<filter::All>>,
-          Collider<He3T_He4D::CrossSectionInterpolator, He3T_He4D::reactants, He3T_He4D::products, OneFilter<filter::All>>,
-          Collider<TT_He4nn::CrossSectionInterpolator, TT_He4nn::reactants, TT_He4nn::products, OneFilter<filter::All>>,
-          Collider<Bp_He4He4::CrossSectionInterpolator, Bp_He4He4::reactants, Bp_He4He4::products, OneFilter<filter::All>>
-          >;
+   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
 ---------------------
 
-.. code-block:: cpp
-
-    // Memory allocation control
-    constexpr uint32_t cellListChunkSize = TYPICAL_PARTICLES_PER_CELL;
-    
-    // Product weighting thresholds
-    constexpr float_X productMinWeighting = 16.1;
-    constexpr uint32_t maxFmult = 1e6;
-    
-    // Debug flags
-    constexpr bool debugFusion = false;
-    constexpr bool alwaysFuseQ = false;  // Force 100% fusion probability
+.. 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)
 

include/picongpu/param/fusion.param

@@ -49,35 +49,12 @@ namespace picongpu
              */
             constexpr uint32_t cellListChunkSize = TYPICAL_PARTICLES_PER_CELL;
 
-            /** Minimum weighting threshold for fusion product macro-particles.
+            /** Maximum allowed fusion reaction rate multiplier per time step.
              *
-             * Sets the minimum allowed weighting (statistical weight) for newly created
-             * product macro-particles. When the fusion production multiplier (Fmult) would
-             * create products lighter than this threshold, Fmult is reduced to ensure
-             * products have at least this weighting. Lower values allow more, lighter products
-             * but increase memory usage and computation time.
+             * This is used to increase the number of product particles while simultaneously decreasing the weight of
+             * each product particle to improve statistics of fusion products.
              */
-            constexpr float_X productMinWeighting = 16.1;
-
-            /** Maximum fusion production multiplier (Fmult).
-             *
-             * Controls how the consumed reactant weight is split into multiple product
-             * macro-particles. Higher values create more, lighter products (down to the
-             * productMinWeighting threshold). Lower values create fewer, heavier products.
-             *
-             * Algorithm per fused pair:
-             * - Start with Fmult = maxFmult
-             * - Compute productWeighting = minWeighting / Fmult
-             * - If productWeighting < productMinWeighting, reduce Fmult to
-             *   max(1, minWeighting / productMinWeighting) and recompute
-             *
-             * The base productWeighting is then distributed across reactant sites and
-             * product species while conserving total weight.
-             *
-             * Note: This distribution can result in creation of product particles with
-             * weighting smaller than productMinWeighting.
-             */
-            constexpr uint32_t maxFmult = 1'000'000;
+            constexpr float_X maxFmult = 1e6;
 
             /** Enable debug output for fusion processes.
              *

include/picongpu/particles/fusion/InterCollision.hpp

@@ -183,11 +183,7 @@ namespace picongpu::particles::fusion
             worker.sync();
             detail::maxArrayDestroy<false>(worker, nppc, numCellsPerSuperCell);
             // now in nppc[0] we have the maximum number of particles in the supercell
-            onlyMaster(
-                [&]()
-                {
-                    maxNumParticlesInCell = nppc[0];
-                });
+            onlyMaster([&]() { maxNumParticlesInCell = nppc[0]; });
             // don't need sync
 
             // --- 4. Shuffle Particle Lists ---
@@ -480,7 +476,7 @@ namespace picongpu::particles::fusion
                 for(uint32_t chunkStart = 0; chunkStart < maxNumParticles; chunkStart += numPairsAtOnce)
                 {
                     // Parallel grid-stride loop over the current chunk
-                    constexpr uint32_t step = std::min(worker.numWorkers(), numPairsAtOnce);
+                    uint32_t step = std::min(worker.numWorkers(), numPairsAtOnce);
                     for(int i = chunkStart + worker.workerIdx();
                         i < chunkStart + numPairsAtOnce && i < maxNumParticles;
                         i += step)
@@ -499,7 +495,7 @@ namespace picongpu::particles::fusion
                         // P = n_min * n_a / n_ba * minWeighting * dt * (Fmult*sigma*v_rel*gamma_cm) <- inside fuse()
                         float_X const probabilityCorrectionFactor
                             = minReactantDensity * correctionFactor[cellIdx] * sim.pic.getDt();
-                        
+
                         // Fusion multiplier - can be changed inside fuse() if probability > 1
                         float_X Fmult = maxFmult;
 
@@ -525,7 +521,7 @@ namespace picongpu::particles::fusion
                         {
                             // because we could change Fmult inside fuser (because the probability might have been >1)
                             float_X productWeighting = minWeighting / Fmult;
-                            
+
                             weightingArray[i] = productWeighting; // no atomic needed because i is unique per thread
 
                             uint32_t freeIndex = alpaka::atomicAdd(
@@ -614,17 +610,25 @@ namespace picongpu::particles::fusion
                         // no fusion happened for this pair
                         if(weightingArray[i] == 0._X)
                             continue;
-                            
+
                         if constexpr(debugFusion)
-                            if(weightingArray[i] < 0._X){
-                                printf("Error: negative weighting in fusion reaction! weightingArray[%d] = %f\n", i, weightingArray[i]);
+                            if(weightingArray[i] < 0._X)
+                            {
+                                printf(
+                                    "Error: negative weighting in fusion reaction! weightingArray[%d] = %f\n",
+                                    i,
+                                    weightingArray[i]);
                                 // print fmult and weightingR1 and weightingR2
                                 float_X weightingR1 = accessor1[i % size1][weighting_];
                                 float_X weightingR2 = accessor2[i % size2][weighting_];
                                 printf("weightingR1 = %f, weightingR2 = %f\n", weightingR1, weightingR2);
                                 printf("Fmult = %f\n", maxFmult);
                                 // print the number of particles in longer list
-                                printf("size1 = %d, size2 = %d, weightingArraySize = %d\n", size1, size2, weightingArraySize);
+                                printf(
+                                    "size1 = %d, size2 = %d, weightingArraySize = %d\n",
+                                    size1,
+                                    size2,
+                                    weightingArraySize);
                                 continue;
                             }
                         float_X const oldWeighting1 = accessor1[i % size1][weighting_];
@@ -642,32 +646,35 @@ namespace picongpu::particles::fusion
                         accessor1[i % size1][multiMask_] = (accessor1[i % size1][weighting_] > 1e-6_X);
                         accessor2[i % size2][multiMask_] = (accessor2[i % size2][weighting_] > 1e-6_X);
 
-                        
+
                         // print i and weighting array
                         if constexpr(debugFusion)
-                            if((((i==0 && alwaysFuseQ) || !alwaysFuseQ) || accessor1[i % size1][multiMask_]==0 || accessor2[i % size2][multiMask_]==0)){
-
-                            printf("worker %d: cell %d, i %d, weightingArray %f, new weighting1 %f, new weighting2 %f, old weighting1 %f, old weighting2 %f, difference 1 %f, difference 2 %f\n",
-                                worker.workerIdx(),
-                                cellIdx,
-                                i,
-                                weightingArray[i],
-                                accessor1[i % size1][weighting_],
-                                accessor2[i % size2][weighting_],
-                                oldWeighting1,
-                                oldWeighting2,
-                                oldWeighting1 - accessor1[i % size1][weighting_],
-                                oldWeighting2 - accessor2[i % size2][weighting_]);
-                                
+                            if((((i == 0 && alwaysFuseQ) || !alwaysFuseQ) || accessor1[i % size1][multiMask_] == 0
+                                || accessor2[i % size2][multiMask_] == 0))
+                            {
+                                printf(
+                                    "worker %d: cell %d, i %d, weightingArray %f, new weighting1 %f, new weighting2 "
+                                    "%f, old weighting1 %f, old weighting2 %f, difference 1 %f, difference 2 %f\n",
+                                    worker.workerIdx(),
+                                    cellIdx,
+                                    i,
+                                    weightingArray[i],
+                                    accessor1[i % size1][weighting_],
+                                    accessor2[i % size2][weighting_],
+                                    oldWeighting1,
+                                    oldWeighting2,
+                                    oldWeighting1 - accessor1[i % size1][weighting_],
+                                    oldWeighting2 - accessor2[i % size2][weighting_]);
+
                                 // print the multimask as well
-                                printf("worker %d: cell %d, i %d, multiMask1 %d, multiMask2 %d\n",
+                                printf(
+                                    "worker %d: cell %d, i %d, multiMask1 %d, multiMask2 %d\n",
                                     worker.workerIdx(),
                                     cellIdx,
                                     i,
                                     accessor1[i % size1][multiMask_],
                                     accessor2[i % size2][multiMask_]);
-                                
-                        }
+                            }
                     }
                     worker.sync();
                 }

include/picongpu/particles/fusion/detail/arrayHelpers.hpp

@@ -52,13 +52,13 @@ namespace picongpu::particles::fusion::detail
 
     template<bool debug = false, typename T_worker, typename T_arr>
     DINLINE void maxArrayDestroy(T_worker const& worker, T_arr& arr, int const& size)
-    {   
+    {
         if(worker.workerIdx() == 0)
         {
             auto maxVal = arr[0];
             for(int i = 1; i < size; ++i)
             {
-                if (arr[i] > maxVal)
+                if(arr[i] > maxVal)
                 {
                     maxVal = arr[i];
                 }

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

@@ -309,18 +309,18 @@ namespace picongpu::particles::fusion::relativistic
 
                 float_X sigma_picArea = crossSection(fusionVar.E_r * convToKeV) * millibarn_to_picArea;
                 float_X P = probabilityFactor * sigma_picArea * fusionVar.V_rel_mag * fusionVar.gamma_cm;
-                
-                if(P > 1.0_COLL){
+
+                if(P > 1.0_COLL)
+                {
                     // print warning
-                    printf(
-                        "Warning: Fusion probability exceeds 1.0 (P: %e) the process will be underestimated.\n",
-                        P);
-                    P= 1.0_COLL;
+                    printf("Warning: Fusion probability exceeds 1.0 (P: %e) the process will be underestimated.\n", P);
+                    P = 1.0_COLL;
                     Fmult = 1.0_COLL;
                 }
-                else if(P>0._COLL){
-                // limit P to not more than 0.99_COLL
-                    float_X const maxFmult = 0.99_COLL/P;
+                else if(P > 0._COLL)
+                {
+                    // limit P to not more than 0.99_COLL
+                    float_X const maxFmult = 0.99_COLL / P;
                     Fmult = std::min(Fmult, maxFmult);
                     P *= Fmult;
                 }
@@ -412,15 +412,14 @@ namespace picongpu::particles::fusion::relativistic
                     fusionVar.P<T_Product0Box, T_Product1Box>(dir);
                     mom0 = fusionVar.P0();
                     mom1 = fusionVar.P1();
-                    
+
                     if constexpr(debugFusion)
                     {
-                        if((mom0 == float3_X{0.0_X, 0.0_X, 0.0_X})
-                        || (mom1 == float3_X{0.0_X, 0.0_X, 0.0_X})){
+                        if((mom0 == float3_X{0.0_X, 0.0_X, 0.0_X}) || (mom1 == float3_X{0.0_X, 0.0_X, 0.0_X}))
+                        {
                             printf("Error: Fusion produced zero momentum products.\n");
                         }
                     }
-                        
                 }
                 else
                 {