Skip to content

Commit 2419daf

Browse files
RemiLehelucafedeli88NeilZaim
authored and
dpgrote
committed
D-T fusion (BLAST-WarpX#3153)
* initial work * fixed bugs and added species * update documentation * delete unused file * Add properties for neutron, hydrogen isotopes, helium isotopes * Update code to be more consistent * Correct typo * Parse deuterium-tritium fusion * Start putting in place the files for deuterium-tritium * Update documentation * Prepare structures for deuterium tritium * Fix typo * Fix compilation * Add neutron * Add correct formula for the cross-section * Correct compilation error * Fix nuclear fusion * Reset benchmarks * Prepare creation functor for 2-product fusion * First implementation of momentum initialization * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Use utility function for fusion * Minor modification of variable names * Fix GPU compilation * Fix single precision compilation * Update types * Use util function in P-B fusion * Correct compilation errors * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Correct errors * Update values of mass and charge * Correct compilation error * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Correct compilation error * Correct compilation error * Correct compilation error * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Reset benchmark * Use helium particle in proton-boron, to avoid resetting benchmark * Fixed proton-boron test * Revert "Fixed proton-boron test" This reverts commit 73c8d9d. * Incorporate Neil's recommendations * Reset benchmarks * Correct compilation errors * Add new deuterium tritium automated test * Correct formula of cross-section * Correct cross-section * Improve analysis script * Add test of energy conservation * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Add test of conservation of momentum * Progress in analysis script * Fix error in the initial energy of the deuterium particles * Add check of isotropy * Clean up the test script * Rewrite p_sq formula in a way to avoids machine-precision negative numbers * Add checksum * Clean up code * Apply suggestions from code review * Update PR according to comments * Update benchmark * Address additional comments * Numerical Literals Co-authored-by: Luca Fedeli <[email protected]> Co-authored-by: Neïl Zaim <[email protected]>
1 parent 69f9c2e commit 2419daf

15 files changed

+868
-8
lines changed

‎Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py‎

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

‎Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_3d‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
#################################
2+
####### GENERAL PARAMETERS ######
3+
#################################
4+
## With these parameters, each cell has a size of exactly 1 by 1 by 1
5+
max_step = 1
6+
amr.n_cell = 8 8 16
7+
amr.max_grid_size = 8
8+
amr.blocking_factor = 8
9+
amr.max_level = 0
10+
geometry.dims = 3
11+
geometry.prob_lo = 0. 0. 0.
12+
geometry.prob_hi = 8. 8. 16.
13+
14+
#################################
15+
###### Boundary Condition #######
16+
#################################
17+
boundary.field_lo = periodic periodic periodic
18+
boundary.field_hi = periodic periodic periodic
19+
20+
#################################
21+
############ NUMERICS ###########
22+
#################################
23+
warpx.verbose = 1
24+
warpx.cfl = 1.0
25+
26+
# Order of particle shape factors
27+
algo.particle_shape = 1
28+
29+
#################################
30+
############ PLASMA #############
31+
#################################
32+
particles.species_names = deuterium1 tritium1 helium1 neutron1 deuterium2 tritium2 helium2 neutron2
33+
34+
my_constants.m_deuterium = 2.01410177812*m_u
35+
my_constants.m_tritium = 3.0160492779*m_u
36+
my_constants.m_reduced = m_deuterium*m_tritium/(m_deuterium+m_tritium)
37+
my_constants.keV_to_J = 1.e3*q_e
38+
my_constants.Energy_step = 22. * keV_to_J
39+
40+
deuterium1.species_type = deuterium
41+
deuterium1.injection_style = "NRandomPerCell"
42+
deuterium1.num_particles_per_cell = 10000
43+
deuterium1.profile = constant
44+
deuterium1.density = 1.
45+
deuterium1.momentum_distribution_type = "parse_momentum_function"
46+
deuterium1.momentum_function_ux(x,y,z) = 0.
47+
deuterium1.momentum_function_uy(x,y,z) = 0.
48+
## Thanks to the floor, all particles in the same cell have the exact same momentum
49+
deuterium1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_deuterium*clight)
50+
deuterium1.do_not_push = 1
51+
deuterium1.do_not_deposit = 1
52+
53+
tritium1.species_type = tritium
54+
tritium1.injection_style = "NRandomPerCell"
55+
tritium1.num_particles_per_cell = 10000
56+
tritium1.profile = constant
57+
tritium1.density = 1.
58+
tritium1.momentum_distribution_type = "parse_momentum_function"
59+
tritium1.momentum_function_ux(x,y,z) = 0.
60+
tritium1.momentum_function_uy(x,y,z) = 0.
61+
## Thanks to the floor, all particles in the same cell have the exact same momentum
62+
tritium1.momentum_function_uz(x,y,z) = -sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_tritium*clight)
63+
tritium1.do_not_push = 1
64+
tritium1.do_not_deposit = 1
65+
66+
helium1.species_type = helium4
67+
helium1.do_not_push = 1
68+
helium1.do_not_deposit = 1
69+
70+
neutron1.species_type = neutron
71+
neutron1.do_not_push = 1
72+
neutron1.do_not_deposit = 1
73+
74+
my_constants.background_dens = 1.e26
75+
my_constants.beam_dens = 1.e20
76+
77+
deuterium2.species_type = deuterium
78+
deuterium2.injection_style = "NRandomPerCell"
79+
deuterium2.num_particles_per_cell = 1000
80+
deuterium2.profile = "parse_density_function"
81+
## A tenth of the macroparticles in each cell is made of immobile high-density background deuteriums.
82+
## The other nine tenths are made of fast low-density beam deuteriums.
83+
deuterium2.density_function(x,y,z) = if(y - floor(y) < 0.1, 10.*background_dens, 10./9.*beam_dens)
84+
deuterium2.momentum_distribution_type = "parse_momentum_function"
85+
deuterium2.momentum_function_ux(x,y,z) = 0.
86+
deuterium2.momentum_function_uy(x,y,z) = 0.
87+
deuterium2.momentum_function_uz(x,y,z) = "if(y - floor(y) < 0.1,
88+
0., sqrt(2*m_deuterium*Energy_step*(floor(z)**2))/(m_deuterium*clight))"
89+
deuterium2.do_not_push = 1
90+
deuterium2.do_not_deposit = 1
91+
92+
tritium2.species_type = tritium
93+
tritium2.injection_style = "NRandomPerCell"
94+
tritium2.num_particles_per_cell = 100
95+
tritium2.profile = constant
96+
tritium2.density = background_dens
97+
tritium2.momentum_distribution_type = "constant"
98+
tritium2.do_not_push = 1
99+
tritium2.do_not_deposit = 1
100+
101+
helium2.species_type = helium4
102+
helium2.do_not_push = 1
103+
helium2.do_not_deposit = 1
104+
105+
neutron2.species_type = neutron
106+
neutron2.do_not_push = 1
107+
neutron2.do_not_deposit = 1
108+
109+
#################################
110+
############ COLLISION ##########
111+
#################################
112+
collisions.collision_names = DTF1 DTF2
113+
114+
DTF1.species = deuterium1 tritium1
115+
DTF1.product_species = helium1 neutron1
116+
DTF1.type = nuclearfusion
117+
DTF1.fusion_multiplier = 1.e50
118+
119+
DTF2.species = deuterium2 tritium2
120+
DTF2.product_species = helium2 neutron2
121+
DTF2.type = nuclearfusion
122+
DTF2.fusion_multiplier = 1.e15
123+
DTF2.fusion_probability_target_value = 0.02
124+
125+
# Diagnostics
126+
diagnostics.diags_names = diag1
127+
diag1.intervals = 1
128+
diag1.diag_type = Full
129+
diag1.fields_to_plot = rho

‎Examples/Modules/nuclear_fusion/inputs_2d‎ renamed to ‎Examples/Modules/nuclear_fusion/inputs_proton_boron_2d‎

File renamed without changes.

‎Examples/Modules/nuclear_fusion/inputs_3d‎ renamed to ‎Examples/Modules/nuclear_fusion/inputs_proton_boron_3d‎

File renamed without changes.

‎Regression/Checksum/benchmarks_json/Deuterium_Tritium_Fusion_3D.json‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
{
2+
"deuterium1": {
3+
"particle_cpu": 5120000.0,
4+
"particle_id": 26214405120000.0,
5+
"particle_momentum_x": 0.0,
6+
"particle_momentum_y": 0.0,
7+
"particle_momentum_z": 2.8875978729147693e-13,
8+
"particle_position_x": 40960140.72983793,
9+
"particle_position_y": 40959772.69310104,
10+
"particle_position_z": 81919021.52308556,
11+
"particle_weight": 1024.000000000021
12+
},
13+
"deuterium2": {
14+
"particle_cpu": 512000.0,
15+
"particle_id": 10747904512000.0,
16+
"particle_momentum_x": 0.0,
17+
"particle_momentum_y": 0.0,
18+
"particle_momentum_z": 3.356807324363973e-14,
19+
"particle_position_x": 4095630.698135355,
20+
"particle_position_y": 4096073.5517983637,
21+
"particle_position_z": 8191737.5566503005,
22+
"particle_weight": 1.0227810240779905e+29
23+
},
24+
"helium1": {
25+
"particle_cpu": 20564.0,
26+
"particle_id": 414455929591.0,
27+
"particle_momentum_x": 1.7519716491839538e-15,
28+
"particle_momentum_y": 1.7523289312260283e-15,
29+
"particle_momentum_z": 1.7480231586369996e-15,
30+
"particle_position_x": 154379.32401483235,
31+
"particle_position_y": 152618.63815943015,
32+
"particle_position_z": 325970.4138010667,
33+
"particle_weight": 4.421535775967805e-28
34+
},
35+
"helium2": {
36+
"particle_cpu": 18146.0,
37+
"particle_id": 370840895571.0,
38+
"particle_momentum_x": 1.5330942227771018e-15,
39+
"particle_momentum_y": 1.5328473121602395e-15,
40+
"particle_momentum_z": 1.7635828326228758e-15,
41+
"particle_position_x": 137011.89739173267,
42+
"particle_position_y": 136605.24328988983,
43+
"particle_position_z": 290143.4673994485,
44+
"particle_weight": 5.756530048087129e+18
45+
},
46+
"lev=0": {
47+
"rho": 0.0
48+
},
49+
"neutron1": {
50+
"particle_cpu": 20564.0,
51+
"particle_id": 415194438443.0,
52+
"particle_momentum_x": 1.7519716491839538e-15,
53+
"particle_momentum_y": 1.7523289312260283e-15,
54+
"particle_momentum_z": 1.7480231586369996e-15,
55+
"particle_position_x": 154379.32401483235,
56+
"particle_position_y": 152618.63815943015,
57+
"particle_position_z": 325970.4138010667,
58+
"particle_weight": 4.421535775967805e-28
59+
},
60+
"neutron2": {
61+
"particle_cpu": 18146.0,
62+
"particle_id": 371427197911.0,
63+
"particle_momentum_x": 1.5330942227771018e-15,
64+
"particle_momentum_y": 1.5328473121602395e-15,
65+
"particle_momentum_z": 1.549297051563983e-15,
66+
"particle_position_x": 137011.89739173267,
67+
"particle_position_y": 136605.24328988983,
68+
"particle_position_z": 290143.4673994485,
69+
"particle_weight": 5.756530048087129e+18
70+
},
71+
"tritium1": {
72+
"particle_cpu": 5120000.0,
73+
"particle_id": 78643205120000.0,
74+
"particle_momentum_x": 0.0,
75+
"particle_momentum_y": 0.0,
76+
"particle_momentum_z": 2.8875978729147693e-13,
77+
"particle_position_x": 40958301.591654316,
78+
"particle_position_y": 40961136.14476712,
79+
"particle_position_z": 81920546.19181262,
80+
"particle_weight": 1024.000000000021
81+
},
82+
"tritium2": {
83+
"particle_cpu": 51200.0,
84+
"particle_id": 1103626291200.0,
85+
"particle_momentum_x": 0.0,
86+
"particle_momentum_y": 0.0,
87+
"particle_momentum_z": 0.0,
88+
"particle_position_x": 409798.0158217681,
89+
"particle_position_y": 409670.9858143465,
90+
"particle_position_z": 819255.8152412223,
91+
"particle_weight": 1.0239999999424347e+29
92+
}
93+
}

‎Regression/WarpX-tests.ini‎

Lines changed: 18 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2170,7 +2170,7 @@ aux1File = Regression/PostProcessingUtils/post_processing_utils.py
21702170

21712171
[Proton_Boron_Fusion_3D]
21722172
buildDir = .
2173-
inputFile = Examples/Modules/nuclear_fusion/inputs_3d
2173+
inputFile = Examples/Modules/nuclear_fusion/inputs_proton_boron_3d
21742174
runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1
21752175
dim = 3
21762176
addToCompileString =
@@ -2186,7 +2186,7 @@ analysisRoutine = Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.p
21862186

21872187
[Proton_Boron_Fusion_2D]
21882188
buildDir = .
2189-
inputFile = Examples/Modules/nuclear_fusion/inputs_2d
2189+
inputFile = Examples/Modules/nuclear_fusion/inputs_proton_boron_2d
21902190
runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1
21912191
dim = 2
21922192
addToCompileString =
@@ -2200,6 +2200,22 @@ compileTest = 0
22002200
doVis = 0
22012201
analysisRoutine = Examples/Modules/nuclear_fusion/analysis_proton_boron_fusion.py
22022202

2203+
[Deuterium_Tritium_Fusion_3D]
2204+
buildDir = .
2205+
inputFile = Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_3d
2206+
runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1
2207+
dim = 3
2208+
addToCompileString =
2209+
cmakeSetupOpts = -DWarpX_DIMS=3
2210+
restartTest = 0
2211+
useMPI = 1
2212+
numprocs = 2
2213+
useOMP = 1
2214+
numthreads = 2
2215+
compileTest = 0
2216+
doVis = 0
2217+
analysisRoutine = Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py
2218+
22032219
[Maxwell_Hybrid_QED_solver]
22042220
buildDir = .
22052221
inputFile = Examples/Tests/Maxwell_Hybrid_QED/inputs_2d

‎Source/Particles/Collision/BinaryCollision/BinaryCollision.H‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -219,6 +219,7 @@ public:
219219
amrex::Vector<ParticleTileType*> tile_products;
220220
amrex::Vector<GetParticlePosition> get_position_products;
221221
amrex::Vector<index_type> products_np;
222+
amrex::Vector<amrex::ParticleReal> products_mass;
222223
constexpr int getpos_offset = 0;
223224
for (int i = 0; i < n_product_species; i++)
224225
{
@@ -227,6 +228,7 @@ public:
227228
get_position_products.push_back(GetParticlePosition(ptile_product,
228229
getpos_offset));
229230
products_np.push_back(ptile_product.numParticles());
231+
products_mass.push_back(product_species_vector[i]->getMass());
230232
}
231233
auto tile_products_data = tile_products.data();
232234

@@ -358,7 +360,7 @@ public:
358360
const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
359361
soa_1, soa_1, tile_products_data,
360362
particle_ptr_1, particle_ptr_1, m1, m1,
361-
p_mask, products_np,
363+
products_mass, p_mask, products_np,
362364
copy_species1, copy_species2,
363365
p_pair_indices_1, p_pair_indices_2,
364366
p_pair_reaction_weight);
@@ -519,7 +521,7 @@ public:
519521
const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
520522
soa_1, soa_2, tile_products_data,
521523
particle_ptr_1, particle_ptr_2, m1, m2,
522-
p_mask, products_np,
524+
products_mass, p_mask, products_np,
523525
copy_species1, copy_species2,
524526
p_pair_indices_1, p_pair_indices_2,
525527
p_pair_reaction_weight);

‎Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -12,9 +12,9 @@
1212

1313
#include "Particles/MultiParticleContainer.H"
1414

15-
enum struct CollisionType { ProtonBoronFusion, Undefined };
15+
enum struct CollisionType { DeuteriumTritiumFusion, ProtonBoronFusion, Undefined };
1616

17-
enum struct NuclearFusionType { ProtonBoron, Undefined };
17+
enum struct NuclearFusionType { DeuteriumTritium, ProtonBoron, Undefined };
1818

1919
namespace BinaryCollisionUtils{
2020

‎Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,14 @@ namespace BinaryCollisionUtils{
2626
auto& species1 = mypc->GetParticleContainerFromName(species_names[0]);
2727
auto& species2 = mypc->GetParticleContainerFromName(species_names[1]);
2828

29-
if ((species1.AmIA<PhysicalSpecies::proton>() && species2.AmIA<PhysicalSpecies::boron11>())
29+
if ((species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::hydrogen3>())
30+
||
31+
(species1.AmIA<PhysicalSpecies::hydrogen3>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
32+
)
33+
{
34+
return NuclearFusionType::DeuteriumTritium;
35+
}
36+
else if ((species1.AmIA<PhysicalSpecies::proton>() && species2.AmIA<PhysicalSpecies::boron11>())
3037
||
3138
(species1.AmIA<PhysicalSpecies::boron11>() && species2.AmIA<PhysicalSpecies::proton>())
3239
)
@@ -56,6 +63,8 @@ namespace BinaryCollisionUtils{
5663

5764
CollisionType nuclear_fusion_type_to_collision_type (const NuclearFusionType fusion_type)
5865
{
66+
if (fusion_type == NuclearFusionType::DeuteriumTritium)
67+
return CollisionType::DeuteriumTritiumFusion;
5968
if (fusion_type == NuclearFusionType::ProtonBoron)
6069
return CollisionType::ProtonBoronFusion;
6170
amrex::Abort("Invalid nuclear fusion type");

0 commit comments

Comments
 (0)