D-T fusion (#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 73c8d9d0be8417d5cd08a23daeebbc322c984808.

* 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 <luca.fedeli@cea.fr>
Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com>

4 files changed, 521 insertions, 0 deletions

diff --git a/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py b/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py
new file mode 100755
index 000000000..f218df54d
--- /dev/null
+++ b/Examples/Modules/nuclear_fusion/analysis_deuterium_tritium_fusion.py
@@ -0,0 +1,392 @@
+#! /usr/bin/env python
+# Copyright 2022 Neil Zaim, Remi Lehe
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+import os
+import sys
+
+import yt
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+import numpy as np
+import scipy.constants as scc
+
+## This script performs various checks for the fusion module. The simulation
+## that we check is made of 2 different tests, each with different reactant and product species.
+##
+## The first test is performed in the center of mass frame of the reactant. It could correspond to the
+## physical case of two beams colliding with each other. The kinetic energy of the colliding
+## particles depends on the cell number in the z direction and varies in the few keV to few MeV
+## range. All the particles within a cell have the exact same momentum. The reactant species have the same
+## density and number of particles in this test. The number of product species is much smaller than
+## the initial number of reactants
+##
+## The second test is performed in the rest frame of the second reactant. It corresponds to the
+## physical case of a low density beam colliding with a high-density mixed target. The energy of the
+## beam particles is varied in the few keV to few MeV range, depending on the cell number in the z
+## direction. As in the previous case, all the particles within a cell have the exact same
+## momentum. In this test, there are 100 immobile macroparticles for each reactant per cell, as well as 900
+## of the beam reactant macroparticles per cell. The density of the immobile particles is 6 orders of
+## magnitude higher than the number of beam particles, which means that they have a much higher
+## weight. This test is similar to the example given in section 3 of Higginson et al.,
+## Journal of Computation Physics, 388 439–453 (2019), which was found to be sensitive to the way
+## unsampled pairs are accounted for. As before, the number of product particles is much smaller than
+## the initial number of reactants.
+##
+## In all simulations, we check particle number, charge, momentum and energy conservation and
+## perform basic checks regarding the produced particles. When possible, we also compare the number
+## of produced macroparticles, fusion yield and energy of the produced particles to theoretical
+## values.
+##
+## Please be aware that the relative tolerances are often set empirically in this analysis script,
+## so it would not be surprising that some tolerances need to be increased in the future.
+
+default_tol = 1.e-12 # Default relative tolerance
+
+## Define reactants and products
+reactant_species = ['deuterium', 'tritium']
+product_species = ['helium', 'neutron']
+
+mass = {
+    'deuterium': 2.01410177812*scc.m_u,
+    'tritium': 3.0160492779*scc.m_u,
+    'helium': 4.00260325413*scc.m_u,
+    'neutron': 1.0013784193052508*scc.m_p
+}
+m_reduced = np.product([mass[s] for s in reactant_species])/np.sum([mass[s] for s in reactant_species])
+
+## Some physical parameters
+keV_to_Joule = scc.e*1e3
+MeV_to_Joule = scc.e*1e6
+barn_to_square_meter = 1.e-28
+
+E_fusion = 17.5893*MeV_to_Joule # Energy released during the fusion reaction
+
+## Some numerical parameters for this test
+size_x = 8
+size_y = 8
+size_z = 16
+dV_total = size_x*size_y*size_z # Total simulation volume
+# Volume of a slice corresponding to a single cell in the z direction. In tests 1 and 2, all the
+# particles of a given species in the same slice have the exact same momentum
+dV_slice = size_x*size_y
+dt = 1./(scc.c*np.sqrt(3.))
+# In test 1 and 2, the energy in cells number i (in z direction) is typically Energy_step * i**2
+Energy_step = 22.*keV_to_Joule
+
+def is_close(val1, val2, rtol=default_tol, atol=0.):
+    ## Wrapper around numpy.isclose, used to override the default tolerances.
+    return np.isclose(val1, val2, rtol=rtol, atol=atol)
+
+def add_existing_species_to_dict(yt_ad, data_dict, species_name, prefix, suffix):
+    data_dict[prefix+"_px_"+suffix]  = yt_ad[species_name, "particle_momentum_x"].v
+    data_dict[prefix+"_py_"+suffix]  = yt_ad[species_name, "particle_momentum_y"].v
+    data_dict[prefix+"_pz_"+suffix]  = yt_ad[species_name, "particle_momentum_z"].v
+    data_dict[prefix+"_w_"+suffix]   = yt_ad[species_name, "particle_weight"].v
+    data_dict[prefix+"_id_"+suffix]  = yt_ad[species_name, "particle_id"].v
+    data_dict[prefix+"_cpu_"+suffix] = yt_ad[species_name, "particle_cpu"].v
+    data_dict[prefix+"_z_"+suffix]   = yt_ad[species_name, "particle_position_z"].v
+
+def add_empty_species_to_dict(data_dict, species_name, prefix, suffix):
+    data_dict[prefix+"_px_"+suffix]  = np.empty(0)
+    data_dict[prefix+"_py_"+suffix]  = np.empty(0)
+    data_dict[prefix+"_pz_"+suffix]  = np.empty(0)
+    data_dict[prefix+"_w_"+suffix]   = np.empty(0)
+    data_dict[prefix+"_id_"+suffix]  = np.empty(0)
+    data_dict[prefix+"_cpu_"+suffix] = np.empty(0)
+    data_dict[prefix+"_z_"+suffix]   = np.empty(0)
+
+def add_species_to_dict(yt_ad, data_dict, species_name, prefix, suffix):
+    try:
+        ## If species exist, we add its data to the dictionary
+        add_existing_species_to_dict(yt_ad, data_dict, species_name, prefix, suffix)
+    except yt.utilities.exceptions.YTFieldNotFound:
+        ## If species does not exist, we avoid python crash and add empty arrays to the
+        ## dictionnary.
+        add_empty_species_to_dict(data_dict, species_name, prefix, suffix)
+
+def check_particle_number_conservation(data):
+    # Check consumption of reactants
+    total_w_reactant1_start = np.sum(data[reactant_species[0] + "_w_start"])
+    total_w_reactant1_end   = np.sum(data[reactant_species[0] + "_w_end"])
+    total_w_reactant2_start = np.sum(data[reactant_species[1] + "_w_start"])
+    total_w_reactant2_end   = np.sum(data[reactant_species[1] + "_w_end"])
+    consumed_reactant1 = total_w_reactant1_start - total_w_reactant1_end
+    consumed_reactant2 = total_w_reactant2_start - total_w_reactant2_end
+    assert(consumed_reactant1 >= 0.)
+    assert(consumed_reactant2 >= 0.)
+    ## Check that number of consumed reactants are equal
+    assert_scale = max(total_w_reactant1_start, total_w_reactant2_start)
+    assert(is_close(consumed_reactant1, consumed_reactant2, rtol = 0., atol = default_tol*assert_scale))
+
+    # That the number of products corresponds consumed particles
+    for species_name in product_species:
+        created_product = np.sum(data[species_name + "_w_end"])
+        assert(created_product >= 0.)
+        assert(is_close(total_w_reactant1_start, total_w_reactant1_end + created_product))
+        assert(is_close(total_w_reactant2_start, total_w_reactant2_end + created_product))
+
+def compute_energy_array(data, species_name, suffix, m):
+    ## Relativistic computation of kinetic energy for a given species
+    psq_array = data[species_name+'_px_'+suffix]**2 + data[species_name+'_py_'+suffix]**2 + \
+                data[species_name+'_pz_'+suffix]**2
+    rest_energy = m*scc.c**2
+    return np.sqrt(psq_array*scc.c**2 + rest_energy**2) - rest_energy
+
+def check_energy_conservation(data):
+    total_energy_start = 0
+    for species_name in reactant_species:
+        total_energy_start +=  np.sum( data[species_name + "_w_start"] * \
+            compute_energy_array(data, species_name, "start", mass[species_name]) )
+    total_energy_end = 0
+    for species_name in product_species + reactant_species:
+        total_energy_end +=  np.sum( data[species_name + "_w_end"] * \
+            compute_energy_array(data, species_name, "end", mass[species_name]) )
+    n_fusion_reaction = np.sum(data[product_species[0] + "_w_end"])
+    assert(is_close(total_energy_end,
+                    total_energy_start + n_fusion_reaction*E_fusion,
+                    rtol = 1.e-8))
+
+def check_momentum_conservation(data):
+    total_px_start = 0
+    total_py_start = 0
+    total_pz_start = 0
+    for species_name in reactant_species:
+        total_px_start += np.sum(
+            data[species_name+'_px_start'] * data[species_name+'_w_start'])
+        total_py_start += np.sum(
+            data[species_name+'_py_start'] * data[species_name+'_w_start'])
+        total_pz_start += np.sum(
+            data[species_name+'_pz_start'] * data[species_name+'_w_start'])
+    total_px_end = 0
+    total_py_end = 0
+    total_pz_end = 0
+    for species_name in reactant_species + product_species:
+        total_px_end += np.sum(
+            data[species_name+'_px_end'] * data[species_name+'_w_end'])
+        total_py_end += np.sum(
+            data[species_name+'_py_end'] * data[species_name+'_w_end'])
+        total_pz_end += np.sum(
+            data[species_name+'_pz_end'] * data[species_name+'_w_end'])
+
+    ## Absolute tolerance is needed because sometimes the initial momentum is exactly 0
+    assert(is_close(total_px_start, total_px_end, atol=1.e-15))
+    assert(is_close(total_py_start, total_py_end, atol=1.e-15))
+    assert(is_close(total_pz_start, total_pz_end, atol=1.e-15))
+
+def check_id(data):
+    ## Check that all created particles have unique id + cpu identifier (two particles with
+    ## different cpu can have the same id)
+    for species_name in product_species:
+        complex_id = data[species_name + "_id_end"] + 1j*data[species_name + "_cpu_end"]
+        assert(complex_id.shape == np.unique(complex_id).shape)
+
+def generic_check(data):
+    check_particle_number_conservation(data)
+    check_energy_conservation(data)
+    check_momentum_conservation(data)
+    check_id(data)
+
+def check_isotropy(data, relative_tolerance):
+    ## Checks that the product particles are emitted isotropically
+    for species_name in product_species:
+        average_px_sq = np.average(data[species_name+"_px_end"]*data[species_name+"_px_end"])
+        average_py_sq = np.average(data[species_name+"_py_end"]*data[species_name+"_py_end"])
+        average_pz_sq = np.average(data[species_name+"_pz_end"]*data[species_name+"_pz_end"])
+        assert(is_close(average_px_sq, average_py_sq, rtol = relative_tolerance))
+        assert(is_close(average_px_sq, average_pz_sq, rtol = relative_tolerance))
+
+def check_xy_isotropy(data):
+    ## Checks that the product particles are emitted isotropically in x and y
+    for species_name in product_species:
+        average_px_sq = np.average(data[species_name+"_px_end"]*data[species_name+"_px_end"])
+        average_py_sq = np.average(data[species_name+"_py_end"]*data[species_name+"_py_end"])
+        average_pz_sq = np.average(data[species_name+"_pz_end"]*data[species_name+"_pz_end"])
+        assert(is_close(average_px_sq, average_py_sq, rtol = 5.e-2))
+        assert(average_pz_sq > average_px_sq)
+        assert(average_pz_sq > average_py_sq)
+
+def cross_section( E_keV ):
+    ## Returns cross section in b, using the analytical fits given
+    ## in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611
+    joule_to_keV = 1.e-3/scc.e
+    B_G = scc.pi * scc.alpha * np.sqrt( 2.*m_reduced * scc.c**2 * joule_to_keV );
+    A1 = 6.927e4;
+    A2 = 7.454e8;
+    A3 = 2.050e6;
+    A4 = 5.2002e4;
+    B1 = 6.38e1;
+    B2 = -9.95e-1;
+    B3 = 6.981e-5;
+    B4 = 1.728e-4;
+    astrophysical_factor = (A1 + E_keV*(A2 + E_keV*(A3 + E_keV*A4))) / (1 + E_keV*(B1 + E_keV*(B2 + E_keV*(B3 + E_keV*B4))));
+    millibarn_to_barn = 1.e-3;
+    return millibarn_to_barn * astrophysical_factor/E_keV * np.exp(-B_G/np.sqrt(E_keV))
+
+def E_com_to_p_sq_com(m1, m2, E):
+    ## E is the total (kinetic+mass) energy of a two particle (with mass m1 and m2) system in
+    ## its center of mass frame, in J.
+    ## Returns the square norm of the momentum of each particle in that frame.
+    E_ratio = E/((m1+m2)*scc.c**2)
+    return m1*m2*scc.c**2 * (E_ratio**2 - 1) + (m1-m2)**2*scc.c**2/4 * (E_ratio - 1./E_ratio)**2
+
+def compute_relative_v_com(E):
+    ## E is the kinetic energy of reactants in the center of mass frame, in keV
+    ## Returns the relative velocity between reactants in this frame, in m/s
+    m0 = mass[reactant_species[0]]
+    m1 = mass[reactant_species[1]]
+    E_J  = E*keV_to_Joule + (m0 + m1)*scc.c**2
+    p_sq = E_com_to_p_sq_com(m0, m1, E_J)
+    p = np.sqrt(p_sq)
+    gamma0 = np.sqrt(1. + p_sq / (m0*scc.c)**2)
+    gamma1 = np.sqrt(1. + p_sq / (m1*scc.c)**2)
+    v0 = p/(gamma0*m0)
+    v1 = p/(gamma1*m1)
+    return v0+v1
+
+def expected_weight_com(E_com, reactant0_density, reactant1_density, dV, dt):
+    ## Computes expected number of product particles as a function of energy E_com in the
+    ## center of mass frame. E_com is in keV.
+    assert(np.all(E_com>=0))
+    ## Case E_com == 0 is handled manually to avoid division by zero
+    conditions = [E_com == 0, E_com > 0]
+    ## Necessary to avoid division by 0 warning when pb_cross_section is evaluated
+    E_com_never_zero = np.clip(E_com, 1.e-15, None)
+    choices = [0., cross_section(E_com_never_zero)*compute_relative_v_com(E_com_never_zero)]
+    sigma_times_vrel = np.select(conditions, choices)
+    return reactant0_density*reactant1_density*sigma_times_vrel*barn_to_square_meter*dV*dt
+
+def check_macroparticle_number(data, fusion_probability_target_value, num_pair_per_cell):
+    ## Checks that the number of macroparticles is as expected for the first and second tests
+
+    ## The first slice 0 < z < 1 does not contribute to product species creation
+    numcells = dV_total - dV_slice
+    ## In these tests, the fusion_multiplier is so high that the fusion probability per pair is
+    ## equal to the parameter fusion_probability_target_value
+    fusion_probability_per_pair = fusion_probability_target_value
+    expected_fusion_number = numcells*num_pair_per_cell*fusion_probability_per_pair
+    expected_macroparticle_number = 2*expected_fusion_number
+    std_macroparticle_number = 2*np.sqrt(expected_fusion_number)
+    actual_macroparticle_number = data[product_species[0] + "_w_end"].shape[0]
+    # 5 sigma test that has an intrinsic probability to fail of 1 over ~2 millions
+    assert(is_close(actual_macroparticle_number, expected_macroparticle_number, rtol = 0.,
+                    atol = 5.*std_macroparticle_number))
+
+    ## used in subsequent function
+    return expected_fusion_number
+
+def p_sq_reactant1_frame_to_E_COM_frame(p_reactant0_sq):
+    # Takes the reactant0 square norm of the momentum in the reactant1 rest frame and returns the total
+    # kinetic energy in the center of mass frame. Everything is in SI units.
+    m0 = mass[reactant_species[0]]
+    m1 = mass[reactant_species[1]]
+
+    # Total (kinetic + mass) energy in lab frame
+    E_lab = np.sqrt(p_reactant0_sq*scc.c**2 + (m0*scc.c**2)**2) + m1*scc.c**2
+    # Use invariant E**2 - p**2c**2 of 4-momentum norm to compute energy in center of mass frame
+    E_com = np.sqrt(E_lab**2 - p_reactant0_sq*scc.c**2)
+    # Corresponding kinetic energy
+    E_com_kin = E_com - (m1+m0)*scc.c**2
+    return E_com_kin*(p_reactant0_sq>0.)
+
+def p_sq_to_kinetic_energy(p_sq, m):
+    ## Returns the kinetic energy of a particle as a function of its squared momentum.
+    ## Everything is in SI units.
+    return np.sqrt(p_sq*scc.c**2 + (m*scc.c**2)**2) - (m*scc.c**2)
+
+def compute_E_com1(data):
+    ## Computes kinetic energy (in Joule) in the center of frame for the first test
+
+    ## Square norm of the momentum of reactant as a function of cell number in z direction
+    p_sq = 2.*m_reduced*(Energy_step*np.arange(size_z)**2)
+    Ekin = 0
+    for species_name in reactant_species:
+        Ekin += p_sq_to_kinetic_energy( p_sq, mass[species_name] )
+    return Ekin
+
+def compute_E_com2(data):
+    ## Computes kinetic energy (in Joule) in the center of frame for the second test
+
+    ## Square norm of the momentum of reactant0 as a function of cell number in z direction
+    p_reactant0_sq = 2.*mass[reactant_species[0]]*(Energy_step*np.arange(size_z)**2)
+    return p_sq_reactant1_frame_to_E_COM_frame(p_reactant0_sq)
+
+def check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density, reactant1_density):
+    ## Checks that the fusion yield is as expected for the first and second tests.
+    product_weight_theory = expected_weight_com(E_com/keV_to_Joule,
+        reactant0_density, reactant1_density, dV_slice, dt)
+    for species_name in product_species:
+        product_weight_simulation = np.histogram(data[species_name+"_z_end"],
+                bins=size_z, range=(0, size_z), weights = data[species_name+"_w_end"])[0]
+        ## -1 is here because the first slice 0 < z < 1 does not contribute to fusion
+        expected_fusion_number_per_slice = expected_fusion_number/(size_z-1)
+        relative_std_weight = 1./np.sqrt(expected_fusion_number_per_slice)
+
+        # 5 sigma test that has an intrinsic probability to fail of 1 over ~2 millions
+        assert(np.all(is_close(product_weight_theory, product_weight_simulation,
+                               rtol = 5.*relative_std_weight)))
+
+def specific_check1(data):
+    check_isotropy(data, relative_tolerance = 3.e-2)
+    expected_fusion_number = check_macroparticle_number(data,
+                                                        fusion_probability_target_value = 0.002,
+                                                        num_pair_per_cell = 10000)
+    E_com = compute_E_com1(data)
+    check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density = 1.,
+                                                           reactant1_density = 1.)
+
+def specific_check2(data):
+    check_xy_isotropy(data)
+    ## Only 900 particles pairs per cell here because we ignore the 10% of reactants that are at rest
+    expected_fusion_number = check_macroparticle_number(data,
+                                                        fusion_probability_target_value = 0.02,
+                                                        num_pair_per_cell = 900)
+    E_com = compute_E_com2(data)
+    check_fusion_yield(data, expected_fusion_number, E_com, reactant0_density = 1.e20,
+                                                           reactant1_density = 1.e26)
+
+def check_charge_conservation(rho_start, rho_end):
+    assert(np.all(is_close(rho_start, rho_end, rtol=2.e-11)))
+
+def main():
+    filename_end = sys.argv[1]
+    filename_start = filename_end[:-4] + '0000'
+    ds_end = yt.load(filename_end)
+    ds_start = yt.load(filename_start)
+    ad_end = ds_end.all_data()
+    ad_start = ds_start.all_data()
+    field_data_end = ds_end.covering_grid(level=0, left_edge=ds_end.domain_left_edge,
+                                          dims=ds_end.domain_dimensions)
+    field_data_start = ds_start.covering_grid(level=0, left_edge=ds_start.domain_left_edge,
+                                              dims=ds_start.domain_dimensions)
+
+    ntests = 2
+    for i in range(1, ntests+1):
+        data = {}
+
+        for species_name in reactant_species:
+            add_species_to_dict(ad_start, data, species_name+str(i), species_name, "start")
+            add_species_to_dict(ad_end, data, species_name+str(i), species_name, "end")
+
+        for species_name in product_species:
+            add_species_to_dict(ad_end, data, species_name+str(i), species_name, "end")
+
+        # General checks that are performed for all tests
+        generic_check(data)
+
+        # Checks that are specific to test number i
+        eval("specific_check"+str(i)+"(data)")
+
+    rho_start = field_data_start["rho"].to_ndarray()
+    rho_end = field_data_end["rho"].to_ndarray()
+    check_charge_conservation(rho_start, rho_end)
+
+    test_name = os.path.split(os.getcwd())[1]
+    checksumAPI.evaluate_checksum(test_name, filename_end)
+
+if __name__ == "__main__":
+    main()
diff --git a/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_3d b/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_3d
new file mode 100644
index 000000000..575e9c428
--- /dev/null
+++ b/Examples/Modules/nuclear_fusion/inputs_deuterium_tritium_3d
@@ -0,0 +1,129 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+## With these parameters, each cell has a size of exactly 1 by 1 by 1
+max_step = 1
+amr.n_cell = 8 8 16
+amr.max_grid_size = 8
+amr.blocking_factor = 8
+amr.max_level = 0
+geometry.dims = 3
+geometry.prob_lo     = 0.    0.    0.
+geometry.prob_hi     = 8.    8.   16.
+
+#################################
+###### Boundary Condition #######
+#################################
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.verbose = 1
+warpx.cfl = 1.0
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = deuterium1 tritium1 helium1 neutron1 deuterium2 tritium2 helium2 neutron2
+
+my_constants.m_deuterium = 2.01410177812*m_u
+my_constants.m_tritium = 3.0160492779*m_u
+my_constants.m_reduced = m_deuterium*m_tritium/(m_deuterium+m_tritium)
+my_constants.keV_to_J = 1.e3*q_e
+my_constants.Energy_step = 22. * keV_to_J
+
+deuterium1.species_type = deuterium
+deuterium1.injection_style = "NRandomPerCell"
+deuterium1.num_particles_per_cell = 10000
+deuterium1.profile = constant
+deuterium1.density = 1.
+deuterium1.momentum_distribution_type = "parse_momentum_function"
+deuterium1.momentum_function_ux(x,y,z) = 0.
+deuterium1.momentum_function_uy(x,y,z) = 0.
+## Thanks to the floor, all particles in the same cell have the exact same momentum
+deuterium1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_deuterium*clight)
+deuterium1.do_not_push = 1
+deuterium1.do_not_deposit = 1
+
+tritium1.species_type = tritium
+tritium1.injection_style = "NRandomPerCell"
+tritium1.num_particles_per_cell = 10000
+tritium1.profile = constant
+tritium1.density = 1.
+tritium1.momentum_distribution_type = "parse_momentum_function"
+tritium1.momentum_function_ux(x,y,z) = 0.
+tritium1.momentum_function_uy(x,y,z) = 0.
+## Thanks to the floor, all particles in the same cell have the exact same momentum
+tritium1.momentum_function_uz(x,y,z) = -sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_tritium*clight)
+tritium1.do_not_push = 1
+tritium1.do_not_deposit = 1
+
+helium1.species_type = helium4
+helium1.do_not_push = 1
+helium1.do_not_deposit = 1
+
+neutron1.species_type = neutron
+neutron1.do_not_push = 1
+neutron1.do_not_deposit = 1
+
+my_constants.background_dens = 1.e26
+my_constants.beam_dens = 1.e20
+
+deuterium2.species_type = deuterium
+deuterium2.injection_style = "NRandomPerCell"
+deuterium2.num_particles_per_cell = 1000
+deuterium2.profile = "parse_density_function"
+## A tenth of the macroparticles in each cell is made of immobile high-density background deuteriums.
+## The other nine tenths are made of fast low-density beam deuteriums.
+deuterium2.density_function(x,y,z) = if(y - floor(y) < 0.1, 10.*background_dens, 10./9.*beam_dens)
+deuterium2.momentum_distribution_type = "parse_momentum_function"
+deuterium2.momentum_function_ux(x,y,z) = 0.
+deuterium2.momentum_function_uy(x,y,z) = 0.
+deuterium2.momentum_function_uz(x,y,z) = "if(y - floor(y) < 0.1,
+                                          0., sqrt(2*m_deuterium*Energy_step*(floor(z)**2))/(m_deuterium*clight))"
+deuterium2.do_not_push = 1
+deuterium2.do_not_deposit = 1
+
+tritium2.species_type = tritium
+tritium2.injection_style = "NRandomPerCell"
+tritium2.num_particles_per_cell = 100
+tritium2.profile = constant
+tritium2.density = background_dens
+tritium2.momentum_distribution_type = "constant"
+tritium2.do_not_push = 1
+tritium2.do_not_deposit = 1
+
+helium2.species_type = helium4
+helium2.do_not_push = 1
+helium2.do_not_deposit = 1
+
+neutron2.species_type = neutron
+neutron2.do_not_push = 1
+neutron2.do_not_deposit = 1
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = DTF1 DTF2
+
+DTF1.species = deuterium1 tritium1
+DTF1.product_species = helium1 neutron1
+DTF1.type = nuclearfusion
+DTF1.fusion_multiplier = 1.e50
+
+DTF2.species = deuterium2 tritium2
+DTF2.product_species = helium2 neutron2
+DTF2.type = nuclearfusion
+DTF2.fusion_multiplier = 1.e15
+DTF2.fusion_probability_target_value = 0.02
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = rho
diff --git a/Examples/Modules/nuclear_fusion/inputs_2d b/Examples/Modules/nuclear_fusion/inputs_proton_boron_2d
index 592f44578..592f44578 100644
--- a/Examples/Modules/nuclear_fusion/inputs_2d
+++ b/Examples/Modules/nuclear_fusion/inputs_proton_boron_2d
diff --git a/Examples/Modules/nuclear_fusion/inputs_3d b/Examples/Modules/nuclear_fusion/inputs_proton_boron_3d
index 73f0d52a5..73f0d52a5 100644
--- a/Examples/Modules/nuclear_fusion/inputs_3d
+++ b/Examples/Modules/nuclear_fusion/inputs_proton_boron_3d