Skip to content

Commit 31cc402

Browse files
FilipO28555optolowicz1
FilipO28555
authored and
optolowicz1
committed
Add fusion implementation DOI:/10.1063/5.0051178
Features: - binary fusion reactions between two different species - choosing product weights depending on particles masses and charge. - documentation
1 parent ed2545b commit 31cc402

25 files changed

+3673
-1
lines changed

‎docs/source/usage/workflows.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,5 +18,6 @@ This section contains typical user workflows and best practices.
1818
workflows/probeParticles
1919
workflows/tracerParticles
2020
workflows/particleFilters
21+
workflows/fusionReactions.rst
2122
workflows/synchrotronExtensionDocumentation
2223
workflows/ionization

‎docs/source/usage/workflows/fusionReactions.rst‎

Lines changed: 454 additions & 0 deletions
Large diffs are not rendered by default.

‎include/picongpu/param/fusion.param‎

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
/* Copyright 2025 Filip Optolowicz
2+
*
3+
* This file is part of PIConGPU.
4+
*
5+
* PIConGPU is free software: you can redistribute it and/or modify
6+
* it under the terms of the GNU General Public License as published by
7+
* the Free Software Foundation, either version 3 of the License, or
8+
* (at your option) any later version.
9+
*
10+
* PIConGPU is distributed in the hope that it will be useful,
11+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
* GNU General Public License for more details.
14+
*
15+
* You should have received a copy of the GNU General Public License
16+
* along with PIConGPU.
17+
* If not, see <http://www.gnu.org/licenses/>.
18+
*/
19+
20+
#pragma once
21+
22+
#include "picongpu/param/particleFilters.param"
23+
#include "picongpu/particles/fusion/fusion.def"
24+
#include "picongpu/plugins/misc/SpeciesFilter.hpp"
25+
26+
#include <pmacc/meta/conversion/MakeSeq.hpp>
27+
28+
namespace picongpu
29+
{
30+
namespace particles
31+
{
32+
namespace fusion
33+
{
34+
//! ===================================================================
35+
//! PRECISION SETTINGS
36+
//! ===================================================================
37+
namespace precision
38+
{
39+
//! Use double precision for fusion collision calculations
40+
using float_COLL = float_64;
41+
} // namespace precision
42+
43+
using namespace picongpu::plugins;
44+
45+
//! ===================================================================
46+
//! FUSION REACTION DEFINITIONS
47+
//! ===================================================================
48+
49+
/** Deuterium-Tritium Fusion Reaction
50+
*
51+
* This defines the D + T -> n + He4 fusion reaction
52+
*
53+
* USAGE: Uncomment this struct to enable D-T fusion in your simulation
54+
*
55+
* The cross-section parameters are based on the Bosch-Hale
56+
* parameterization from experimental data.
57+
* Source: https://hdl.handle.net/11858/00-001M-0000-0027-6535-1 (page 29)
58+
*/
59+
// struct DT
60+
// {
61+
// //! Reactant species: Deuterium and Tritium
62+
// using reactants = pmacc::mp_list<Pair<PIC_Tritons, PIC_Deuterons>>;
63+
//
64+
// //! Product species: Neutron and Helium-4
65+
// using products = pmacc::mp_list<Pair<PIC_Neutrons, PIC_He4>>;
66+
//
67+
// /** Cross-section parameters for D-T fusion */
68+
// struct Params
69+
// {
70+
// //! Gamov constant: BG = pi*alpha*Z1*Z2*sqrt(2*mu*c^2) in keV^(1/2)
71+
// //! where alpha is fine structure constant, Z1,Z2 are charges, mu is reduced mass
72+
// static constexpr float_X BG = 34.3827_X;
73+
//
74+
// //! Bosch-Hale parameterization coefficients for D-T cross-section
75+
// //! These coefficients fit experimental cross-section data
76+
// static constexpr float_X A1 = 6.927e4_X;
77+
// static constexpr float_X A2 = 7.454e8_X;
78+
// static constexpr float_X A3 = 2.050e6_X;
79+
// static constexpr float_X A4 = 5.2002e4_X;
80+
// static constexpr float_X A5 = 0.0_X;
81+
//
82+
// //! Additional Bosch-Hale coefficients
83+
// static constexpr float_X B1 = 6.38e1_X;
84+
// static constexpr float_X B2 = -9.95e-1_X;
85+
// static constexpr float_X B3 = 6.981e-5_X;
86+
// static constexpr float_X B4 = 1.728e-4_X;
87+
// };
88+
// //! For now it's the only cross-section interpolator available
89+
// //! Cross-section calculation using Bosch-Hale parameterization.
90+
// //! Returns cross-section in milli-barns.
91+
// using CrossSectionInterpolator = relativistic::FusionFunctor<Params>;
92+
// };
93+
94+
// empty pipeline - no fusion
95+
using FusionPipeline = pmacc::mp_list<>;
96+
97+
//! ===================================================================
98+
//! FUSION PIPELINE CONFIGURATION
99+
//! ===================================================================
100+
101+
/** FusionPipeline: Defines the sequence of fusion reactions
102+
*
103+
* This controls which fusion reactions occur and in what order.
104+
* Each reaction is processed sequentially from first to last.
105+
*
106+
* USAGE EXAMPLES:
107+
*
108+
* 1. Single D-T reaction:
109+
* using FusionPipeline = pmacc::mp_list<ColliderFromStruct<DT>>;
110+
*
111+
* 2. Multiple reactions (processed in order):
112+
* using FusionPipeline = pmacc::mp_list<
113+
* ColliderFromStruct<DT>,
114+
* ColliderFromStruct<DD>,
115+
* ColliderFromStruct<DHe3>
116+
* >;
117+
*
118+
* 3. Manual collider specification:
119+
* using FusionPipeline = pmacc::mp_list<
120+
* Collider<DT::CrossSectionInterpolator,
121+
* DT::reactants,
122+
* DT::products,
123+
* OneFilter<filter::All>>
124+
* >;
125+
*/
126+
127+
128+
//! Method 1: Using ColliderFromStruct (recommended - cleaner syntax)
129+
// using FusionPipeline = pmacc::mp_list<ColliderFromStruct<DT>>;
130+
131+
//! Method 2: Manual collider specification
132+
// using FusionPipeline = pmacc::mp_list<
133+
// Collider<DT::CrossSectionInterpolator, DT::reactants, DT::products, OneFilter<filter::All>>
134+
// >;
135+
136+
// empty pipeline - no fusion
137+
using FusionPipeline = pmacc::mp_list<>;
138+
139+
//! ===================================================================
140+
//! SIMULATION PARAMETERS
141+
//! ===================================================================
142+
143+
/** Memory Management Parameters */
144+
145+
/** Cell list chunk size for memory allocation
146+
*
147+
* Controls memory allocation strategy for collision algorithm.
148+
* The algorithm allocates multiples of this value to store
149+
* particle ID lists, reducing memory fragmentation on GPUs.
150+
*
151+
* Larger values = less fragmentation but more memory overhead
152+
*/
153+
constexpr uint32_t cellListChunkSize = TYPICAL_PARTICLES_PER_CELL;
154+
155+
/** Product Weighting Parameters */
156+
157+
/** Minimum weighting threshold for fusion products
158+
*
159+
* When fusion product weighting <= productMinWeighting, the
160+
* multiplication factor is set to:
161+
* Fmult = max(1, productWeighting/productMinWeighting)
162+
*
163+
* Fmult is then used to scale the fusion probability:
164+
* larger Fmult means we produce more products with smaller weighting.
165+
*
166+
* note: This does not ensure that the product weighting is always above this threshold - see docs "Fusion
167+
* Particle Creation Algorithm"
168+
*/
169+
constexpr float_X productMinWeighting = 16.1;
170+
171+
/** Maximum multiplication factor for products
172+
* the algorithm tries to use this multiplier and if the product weighting is smaller than
173+
* productMinWeighting, Fusion multiplier will be reduced.
174+
*
175+
* Usually set to a value ~10~100x less than the typical particle weighting in a cell.
176+
*/
177+
constexpr uint32_t maxFmult = 1e6;
178+
179+
/** Debug and Control Flags */
180+
181+
/** Enable detailed fusion debugging output
182+
* When true: Prints detailed information about fusion events - warning: large output!
183+
*/
184+
constexpr bool debugFusion = false;
185+
186+
/** Force fusion probability to 100%
187+
* All collision attempts result in fusion
188+
*/
189+
constexpr bool alwaysFuseQ = false;
190+
} // namespace fusion
191+
} // namespace particles
192+
} // namespace picongpu

‎include/picongpu/particles/fusion/Collider.def‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
/* Copyright 2019-2024 Rene Widera, Pawel Ordyna
2+
*
3+
* This file is part of PIConGPU.
4+
*
5+
* PIConGPU is free software: you can redistribute it and/or modify
6+
* it under the terms of the GNU General Public License as published by
7+
* the Free Software Foundation, either version 3 of the License, or
8+
* (at your option) any later version.
9+
*
10+
* PIConGPU is distributed in the hope that it will be useful,
11+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
* GNU General Public License for more details.
14+
*
15+
* You should have received a copy of the GNU General Public License
16+
* along with PIConGPU.
17+
* If not, see <http://www.gnu.org/licenses/>.
18+
*/
19+
20+
#pragma once
21+
22+
#include "picongpu/particles/filter/All.hpp"
23+
24+
#include <pmacc/meta/Pair.hpp>
25+
26+
namespace picongpu::particles::fusion
27+
{
28+
/* A pair of colliding species.
29+
*
30+
* @tparam T_Species0
31+
* @tparam T_Species1
32+
*/
33+
template<typename T_Species0, typename T_Species1>
34+
using Pair = pmacc::meta::Pair<T_Species0, T_Species1>;
35+
36+
/* A set of particle filters for a pair of colliding species.
37+
*
38+
* @tparam T_Filter0 Filter applied on the 1st species in a pair.
39+
* @tparam T_Filter1 Filter applied on the 2nd species in a pair.
40+
*/
41+
template<typename T_Filter0, typename T_Filter1>
42+
using FilterPair = pmacc::meta::Pair<T_Filter0, T_Filter1>;
43+
44+
/* Sets the same particle filter for both species in a pair.
45+
*
46+
* @tparam T_Filter A common particle filter for both colliding species.
47+
*/
48+
template<typename T_Filter>
49+
using OneFilter = pmacc::meta::Pair<T_Filter, T_Filter>;
50+
51+
/* Defines a set of binary collisions with common parameters.
52+
*
53+
* @tparam T_CrossSectionInterpolator A binary particle functor defining a single
54+
* macro particle collision in the binary-collision algorithm.
55+
* @tparam T_SpeciesPairs A sequence of pairs of colliding species.
56+
* @tparam T_Params A struct defining `coulombLog` for the collisions.
57+
* @tparam T_FilterPair A pair of particle filters. Each for every species
58+
* in a pair of colliding species. This pair of filters will be
59+
* aplied to all pairs in T_SpeciesPairs.
60+
*/
61+
template<
62+
typename T_CrossSectionInterpolator,
63+
typename T_SpeciesPairsReactants,
64+
typename T_SpeciesPairsProducts,
65+
typename T_FilterPair = OneFilter<filter::All>>
66+
struct Collider
67+
{
68+
using Functor = T_CrossSectionInterpolator;
69+
using SpeciesPairsReactants = T_SpeciesPairsReactants;
70+
using SpeciesPairsProducts = T_SpeciesPairsProducts;
71+
using FilterPair = T_FilterPair;
72+
};
73+
74+
template<typename T_FusionStruct>
75+
struct ColliderFromStruct
76+
{
77+
using Functor = typename T_FusionStruct::CrossSectionInterpolator;
78+
using SpeciesPairsReactants = typename T_FusionStruct::reactants;
79+
using SpeciesPairsProducts = typename T_FusionStruct::products;
80+
using FilterPair = typename T_FusionStruct::FilterPair;
81+
};
82+
} // namespace picongpu::particles::fusion

‎include/picongpu/particles/fusion/Collider.hpp‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
/* Copyright 2019-2024 Rene Widera, Pawel Ordyna, Filip Optolowicz
2+
*
3+
* This file is part of PIConGPU.
4+
*
5+
* PIConGPU is free software: you can redistribute it and/or modify
6+
* it under the terms of the GNU General Public License as published by
7+
* the Free Software Foundation, either version 3 of the License, or
8+
* (at your option) any later version.
9+
*
10+
* PIConGPU is distributed in the hope that it will be useful,
11+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
* GNU General Public License for more details.
14+
*
15+
* You should have received a copy of the GNU General Public License
16+
* along with PIConGPU.
17+
* If not, see <http://www.gnu.org/licenses/>.
18+
*/
19+
20+
#pragma once
21+
22+
#include "picongpu/defines.hpp"
23+
#include "picongpu/particles/fusion/Collider.def"
24+
#include "picongpu/particles/fusion/WithPeer.hpp"
25+
26+
#include <pmacc/meta/ForEach.hpp>
27+
#include <pmacc/meta/accessors/First.hpp>
28+
#include <pmacc/meta/accessors/Second.hpp>
29+
#include <pmacc/meta/conversion/ApplyGuard.hpp>
30+
#include <pmacc/meta/conversion/ToSeq.hpp>
31+
32+
namespace picongpu::particles::fusion
33+
{
34+
namespace detail
35+
{
36+
// "For each" implementation for calling a collider for each species pair in a list of reactants
37+
template<
38+
typename T_SpeciesPairListReactants,
39+
typename T_SpeciesPairListProducts,
40+
typename T_Collider,
41+
uint32_t colliderId>
42+
struct CallColliderForAPair
43+
{
44+
template<size_t... I>
45+
HINLINE void operator()(
46+
std::index_sequence<I...>,
47+
std::shared_ptr<DeviceHeap> const& deviceHeap,
48+
uint32_t currentStep)
49+
{
50+
(fusion::WithPeer<
51+
typename T_Collider::Functor,
52+
typename pmacc::mp_at_c<T_SpeciesPairListReactants, I>::first,
53+
typename pmacc::mp_at_c<T_SpeciesPairListReactants, I>::second,
54+
typename pmacc::mp_at_c<T_SpeciesPairListProducts, I>::first,
55+
typename pmacc::mp_at_c<T_SpeciesPairListProducts, I>::second,
56+
typename T_Collider::FilterPair,
57+
colliderId,
58+
I>{}(deviceHeap, currentStep),
59+
...);
60+
}
61+
};
62+
} // namespace detail
63+
64+
template<typename T_Collider, uint32_t colliderId>
65+
struct CallCollider
66+
{
67+
void operator()(std::shared_ptr<DeviceHeap> const& deviceHeap, uint32_t currentStep)
68+
{
69+
using SpeciesPairListReactants = pmacc::ToSeq<typename T_Collider::SpeciesPairsReactants>;
70+
using SpeciesPairListProducts = pmacc::ToSeq<typename T_Collider::SpeciesPairsProducts>;
71+
constexpr size_t numPairs = pmacc::mp_size<SpeciesPairListReactants>::value;
72+
std::make_index_sequence<numPairs> index{};
73+
detail::CallColliderForAPair<SpeciesPairListReactants, SpeciesPairListProducts, T_Collider, colliderId>{}(
74+
index,
75+
deviceHeap,
76+
currentStep);
77+
}
78+
};
79+
} // namespace picongpu::particles::fusion

0 commit comments

Comments
 (0)