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>

9 files changed, 236 insertions, 6 deletions

diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
index fd33c0ff4..b1c1ebd93 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
@@ -219,6 +219,7 @@ public:
         amrex::Vector<ParticleTileType*> tile_products;
         amrex::Vector<GetParticlePosition> get_position_products;
         amrex::Vector<index_type> products_np;
+        amrex::Vector<amrex::ParticleReal> products_mass;
         constexpr int getpos_offset = 0;
         for (int i = 0; i < n_product_species; i++)
         {
@@ -227,6 +228,7 @@ public:
             get_position_products.push_back(GetParticlePosition(ptile_product,
                                                                       getpos_offset));
             products_np.push_back(ptile_product.numParticles());
+            products_mass.push_back(product_species_vector[i]->getMass());
         }
         auto tile_products_data = tile_products.data();
 
@@ -358,7 +360,7 @@ public:
             const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
                                                     soa_1, soa_1, tile_products_data,
                                                     particle_ptr_1, particle_ptr_1, m1, m1,
-                                                    p_mask, products_np,
+                                                    products_mass, p_mask, products_np,
                                                     copy_species1, copy_species2,
                                                     p_pair_indices_1, p_pair_indices_2,
                                                     p_pair_reaction_weight);
@@ -519,7 +521,7 @@ public:
             const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
                                                     soa_1, soa_2, tile_products_data,
                                                     particle_ptr_1, particle_ptr_2, m1, m2,
-                                                    p_mask, products_np,
+                                                    products_mass, p_mask, products_np,
                                                     copy_species1, copy_species2,
                                                     p_pair_indices_1, p_pair_indices_2,
                                                     p_pair_reaction_weight);
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
index c20e9716c..2333e1f41 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H
@@ -12,9 +12,9 @@
 
 #include "Particles/MultiParticleContainer.H"
 
-enum struct CollisionType { ProtonBoronFusion, Undefined };
+enum struct CollisionType { DeuteriumTritiumFusion, ProtonBoronFusion, Undefined };
 
-enum struct NuclearFusionType { ProtonBoron, Undefined };
+enum struct NuclearFusionType { DeuteriumTritium, ProtonBoron, Undefined };
 
 namespace BinaryCollisionUtils{
 
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
index cc0c70fce..4bbdb4bc4 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
@@ -26,7 +26,14 @@ namespace BinaryCollisionUtils{
             auto& species1 = mypc->GetParticleContainerFromName(species_names[0]);
             auto& species2 = mypc->GetParticleContainerFromName(species_names[1]);
 
-            if ((species1.AmIA<PhysicalSpecies::proton>() && species2.AmIA<PhysicalSpecies::boron11>())
+            if ((species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::hydrogen3>())
+                ||
+                (species1.AmIA<PhysicalSpecies::hydrogen3>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
+                )
+            {
+                return NuclearFusionType::DeuteriumTritium;
+            }
+            else if ((species1.AmIA<PhysicalSpecies::proton>() && species2.AmIA<PhysicalSpecies::boron11>())
                 ||
                 (species1.AmIA<PhysicalSpecies::boron11>() && species2.AmIA<PhysicalSpecies::proton>())
                 )
@@ -56,6 +63,8 @@ namespace BinaryCollisionUtils{
 
     CollisionType nuclear_fusion_type_to_collision_type (const NuclearFusionType fusion_type)
         {
+            if (fusion_type == NuclearFusionType::DeuteriumTritium)
+                return CollisionType::DeuteriumTritiumFusion;
             if (fusion_type == NuclearFusionType::ProtonBoron)
                 return CollisionType::ProtonBoronFusion;
             amrex::Abort("Invalid nuclear fusion type");
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionCrossSection.H
new file mode 100644
index 000000000..58aaa986b
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionCrossSection.H
@@ -0,0 +1,64 @@
+/* Copyright 2022 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef DEUTERIUM_TRITIUM_FUSION_CROSS_SECTION_H
+#define DEUTERIUM_TRITIUM_FUSION_CROSS_SECTION_H
+
+#include "Utils/WarpXConst.H"
+
+#include <AMReX_REAL.H>
+
+#include <cmath>
+
+/**
+ * \brief Computes the total deuterium-tritium fusion cross section, using
+ * the analytical fits given in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611
+ *
+ * @param[in] E_kin_star the kinetic energy of the deuterium-tritium pair in its center of mass frame, in SI units.
+ * @return The total cross section in SI units (square meters).
+ */
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+amrex::ParticleReal DeuteriumTritiumFusionCrossSection (const amrex::ParticleReal& E_kin_star)
+{
+    using namespace amrex::literals;
+
+    constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e;
+    const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV;
+
+    // If kinetic energy is 0, return a 0 cross section and avoid later division by 0.
+    if (E_keV == 0._prt) {return 0._prt;}
+
+    // Compute the Gamow constant B_G (in keV^{1/2})
+    // (See Eq. 3 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
+    constexpr amrex::ParticleReal m_D = 2.01410177812 * PhysConst::m_u;
+    constexpr amrex::ParticleReal m_T = 3.0160492779 * PhysConst::m_u;
+    constexpr amrex::ParticleReal m_reduced = m_D / (1._prt + m_D/m_T);
+    amrex::ParticleReal B_G = MathConst::pi * PhysConst::alpha * 1._prt * 1._prt
+         * std::sqrt( 2._prt*m_reduced*PhysConst::c*PhysConst::c * joule_to_keV );
+
+    // Compute astrophysical_factor
+    // (See Eq. 9 and Table IV in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
+    constexpr amrex::ParticleReal A1 = 6.927e4;
+    constexpr amrex::ParticleReal A2 = 7.454e8;
+    constexpr amrex::ParticleReal A3 = 2.050e6;
+    constexpr amrex::ParticleReal A4 = 5.2002e4;
+    constexpr amrex::ParticleReal B1 = 6.38e1;
+    constexpr amrex::ParticleReal B2 = -9.95e-1;
+    constexpr amrex::ParticleReal B3 = 6.981e-5;
+    constexpr amrex::ParticleReal B4 = 1.728e-4;
+
+    amrex::ParticleReal astrophysical_factor =
+        (A1 + E_keV*(A2 + E_keV*(A3 + E_keV*A4))) /
+        (1_prt + E_keV*(B1 + E_keV*(B2 + E_keV*(B3 + E_keV*B4))));
+
+    // Compute cross-section in SI units
+    // (See Eq. 8 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
+    constexpr amrex::ParticleReal millibarn_to_sqm = 1.e-31_prt;
+    return millibarn_to_sqm * astrophysical_factor/E_keV * std::exp(-B_G/std::sqrt(E_keV));
+}
+
+#endif // DEUTERIUM_TRITIUM_TRITIUM_FUSION_CROSS_SECTION_H
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
index be1d9622e..2da560ea3 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
@@ -68,6 +68,19 @@ public:
         amrex::Vector<std::string> product_species_name;
         pp_collision_name.getarr("product_species", product_species_name);
 
+        if (m_fusion_type == NuclearFusionType::DeuteriumTritium)
+        {
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                product_species_name.size() == 2u,
+                "ERROR: Deuterium-tritium fusion must contain exactly two product species");
+            auto& product_species1 = mypc->GetParticleContainerFromName(product_species_name[0]);
+            auto& product_species2 = mypc->GetParticleContainerFromName(product_species_name[1]);
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                (product_species1.AmIA<PhysicalSpecies::helium4>() && product_species2.AmIA<PhysicalSpecies::neutron>())
+                ||
+                (product_species1.AmIA<PhysicalSpecies::neutron>() && product_species2.AmIA<PhysicalSpecies::helium4>()),
+                "ERROR: Product species of deuterium-tritium fusion must be of type neutron and helium4");
+        }
         if (m_fusion_type == NuclearFusionType::ProtonBoron)
         {
             WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
index 8b3bd625b..4df0c583e 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H
@@ -8,6 +8,7 @@
 #ifndef SINGLE_NUCLEAR_FUSION_EVENT_H_
 #define SINGLE_NUCLEAR_FUSION_EVENT_H_
 
+#include "DeuteriumTritiumFusionCrossSection.H"
 #include "ProtonBoronFusionCrossSection.H"
 
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
@@ -111,7 +112,11 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part
 
     // Compute fusion cross section as a function of kinetic energy in the center of mass frame
     auto fusion_cross_section = amrex::ParticleReal(0.);
-    if (fusion_type == NuclearFusionType::ProtonBoron)
+    if (fusion_type == NuclearFusionType::DeuteriumTritium)
+    {
+        fusion_cross_section = DeuteriumTritiumFusionCrossSection(E_kin_star);
+    }
+    else if (fusion_type == NuclearFusionType::ProtonBoron)
     {
         fusion_cross_section = ProtonBoronFusionCrossSection(E_kin_star);
     }
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H
new file mode 100644
index 000000000..be3f5b2d9
--- /dev/null
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H
@@ -0,0 +1,97 @@
+/* Copyright 2022 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef TWO_PRODUCT_FUSION_INITIALIZE_MOMENTUM_H
+#define TWO_PRODUCT_FUSION_INITIALIZE_MOMENTUM_H
+
+#include "TwoProductFusionUtil.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/ParticleUtils.H"
+#include "Utils/WarpXConst.H"
+
+#include <AMReX_DenseBins.H>
+#include <AMReX_Random.H>
+#include <AMReX_REAL.H>
+
+#include <cmath>
+#include <limits>
+
+namespace {
+    // Define shortcuts for frequently-used type names
+    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+    using ParticleType = WarpXParticleContainer::ParticleType;
+    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using index_type = ParticleBins::index_type;
+
+    /**
+     * \brief This function initializes the momentum of the product particles,
+     * in a fusion event where only two products are produced.
+     * (In this case, conservation of energy and momentum determines
+     * the amplitude of the momentum of the particles exactly.)
+     * We assume that the emission of the product is isotropic in the center-of-mass frame
+     *
+     * @param[in] soa1_in struct of array data of the first colliding species
+     * @param[in] soa2_in struct of array data of the second colliding species
+     * @param[out] soa1_out struct of array data of the first product species
+     * @param[out] soa2_out struct of array data of the first product species
+     * @param[in] idx1_in index of first colliding macroparticle
+     * @param[in] idx2_in index of second colliding macroparticle
+     * @param[in] idx1_out_start index of first product macroparticle
+     * @param[in] idx2_out_start index of second product macroparticle
+     * @param[in] m1_in mass of first colliding species
+     * @param[in] m2_in mass of second colliding species
+     * @param[in] m1_out mass of first product species
+     * @param[in] m2_out mass of second product species
+     * @param[in] E_fusion energy released in the fusion reaction
+     * @param[in] engine the random engine
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    void TwoProductFusionInitializeMomentum (
+                            const SoaData_type& soa1_in, const SoaData_type& soa2_in,
+                            SoaData_type& soa1_out, SoaData_type& soa2_out,
+                            const index_type& idx1_in, const index_type& idx2_in,
+                            const index_type& idx1_out_start, const index_type& idx2_out_start,
+                            const amrex::ParticleReal& m1_in, const amrex::ParticleReal& m2_in,
+                            const amrex::ParticleReal& m1_out, const amrex::ParticleReal& m2_out,
+                            const amrex::ParticleReal& E_fusion,
+                            const amrex::RandomEngine& engine)
+    {
+        using namespace amrex::literals;
+
+        amrex::ParticleReal ux1_out = 0.0_prt, uy1_out = 0.0_prt, uz1_out = 0.0_prt;
+        amrex::ParticleReal ux2_out = 0.0_prt, uy2_out = 0.0_prt, uz2_out = 0.0_prt;
+
+        TwoProductFusionComputeProductMomenta(
+            soa1_in.m_rdata[PIdx::ux][idx1_in],
+            soa1_in.m_rdata[PIdx::uy][idx1_in],
+            soa1_in.m_rdata[PIdx::uz][idx1_in], m1_in,
+            soa2_in.m_rdata[PIdx::ux][idx2_in],
+            soa2_in.m_rdata[PIdx::uy][idx2_in],
+            soa2_in.m_rdata[PIdx::uz][idx2_in], m2_in,
+            ux1_out, uy1_out, uz1_out, m1_out,
+            ux2_out, uy2_out, uz2_out, m2_out,
+            E_fusion,
+            engine);
+
+        // Fill momentum of product species (note that we actually
+        // create 4 products, 2 at the position of each incident particle)
+        soa1_out.m_rdata[PIdx::ux][idx1_out_start] = ux1_out;
+        soa1_out.m_rdata[PIdx::uy][idx1_out_start] = uy1_out;
+        soa1_out.m_rdata[PIdx::uz][idx1_out_start] = uz1_out;
+        soa1_out.m_rdata[PIdx::ux][idx1_out_start + 1] = ux1_out;
+        soa1_out.m_rdata[PIdx::uy][idx1_out_start + 1] = uy1_out;
+        soa1_out.m_rdata[PIdx::uz][idx1_out_start + 1] = uz1_out;
+        soa2_out.m_rdata[PIdx::ux][idx2_out_start] = ux2_out;
+        soa2_out.m_rdata[PIdx::uy][idx2_out_start] = uy2_out;
+        soa2_out.m_rdata[PIdx::uz][idx2_out_start] = uz2_out;
+        soa2_out.m_rdata[PIdx::ux][idx2_out_start + 1] = ux2_out;
+        soa2_out.m_rdata[PIdx::uy][idx2_out_start + 1] = uy2_out;
+        soa2_out.m_rdata[PIdx::uz][idx2_out_start + 1] = uz2_out;
+    }
+}
+
+#endif // TWO_PRODUCT_FUSION_INITIALIZE_MOMENTUM_H
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
index 6a8387eef..510c59094 100644
--- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
@@ -10,6 +10,7 @@
 
 #include "BinaryCollisionUtils.H"
 
+#include "Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H"
 #include "Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H"
 #include "Particles/ParticleCreation/SmartCopy.H"
 #include "Particles/MultiParticleContainer.H"
@@ -77,6 +78,7 @@ public:
      *             reaches 0.
      * @param[in] m1 mass of the first colliding particle species
      * @param[in] m2 mass of the second colliding particle species
+     * @param[in] products_mass array storing the mass of product particles
      * @param[in] p_mask a mask that is 1 if binary collision has resulted in particle creation
      *            event, 0 otherwise.
      * @param[in] products_np array storing the number of existing product particles in that tile
@@ -102,6 +104,7 @@ public:
                     ParticleTileType** AMREX_RESTRICT tile_products,
                     ParticleType* particle_ptr_1, ParticleType*  particle_ptr_2,
                     const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
+                    const amrex::Vector<amrex::ParticleReal>& products_mass,
                     const index_type* AMREX_RESTRICT p_mask,
                     const amrex::Vector<index_type>& products_np,
                     const SmartCopy* AMREX_RESTRICT copy_species1,
@@ -142,18 +145,24 @@ public:
 #ifdef AMREX_USE_GPU
         amrex::Gpu::DeviceVector<SoaData_type> device_soa_products(m_num_product_species);
         amrex::Gpu::DeviceVector<index_type> device_products_np(m_num_product_species);
+        amrex::Gpu::DeviceVector<amrex::ParticleReal> device_products_mass(m_num_product_species);
         amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, soa_products.begin(),
                               soa_products.end(),
                               device_soa_products.begin());
         amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, products_np.begin(),
                               products_np.end(),
                               device_products_np.begin());
+       amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, products_mass.begin(),
+                              products_mass.end(),
+                              device_products_mass.begin());
         amrex::Gpu::streamSynchronize();
         SoaData_type* AMREX_RESTRICT soa_products_data = device_soa_products.data();
         const index_type* AMREX_RESTRICT products_np_data = device_products_np.data();
+        const amrex::ParticleReal* AMREX_RESTRICT products_mass_data = device_products_mass.data();
 #else
         SoaData_type* AMREX_RESTRICT soa_products_data = soa_products.data();
         const index_type* AMREX_RESTRICT products_np_data = products_np.data();
+        const amrex::ParticleReal* AMREX_RESTRICT products_mass_data = products_mass.data();
 #endif
 
         const int t_num_product_species = m_num_product_species;
@@ -216,6 +225,20 @@ public:
                                                         p_pair_indices_1[i], p_pair_indices_2[i],
                                                         product_start_index, m1, m2, engine);
                 }
+                else if (t_collision_type == CollisionType::DeuteriumTritiumFusion)
+                {
+                    amrex::ParticleReal fusion_energy = 0.0_prt;
+                    if (t_collision_type == CollisionType::DeuteriumTritiumFusion) {
+                        fusion_energy = 17.5893e6_prt * PhysConst::q_e; // 17.6 MeV
+                    }
+                    TwoProductFusionInitializeMomentum(soa_1, soa_2,
+                        soa_products_data[0], soa_products_data[1],
+                        p_pair_indices_1[i], p_pair_indices_2[i],
+                        products_np_data[0] + 2*p_offsets[i]*p_num_products_device[0],
+                        products_np_data[1] + 2*p_offsets[i]*p_num_products_device[1],
+                        m1, m2, products_mass_data[0], products_mass_data[1], fusion_energy, engine);
+                }
+
             }
         });
 
@@ -260,6 +283,7 @@ public:
                     ParticleTileType** /*tile_products*/,
                     ParticleType* /*particle_ptr_1*/, ParticleType* /*particle_ptr_2*/,
                     const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/,
+                    const amrex::Vector<amrex::ParticleReal>& /*products_mass*/,
                     const index_type* /*p_mask*/, const amrex::Vector<index_type>& /*products_np*/,
                     const SmartCopy* /*copy_species1*/, const SmartCopy* /*copy_species2*/,
                     const index_type* /*p_pair_indices_1*/, const index_type* /*p_pair_indices_2*/,
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
index ba6b20f1b..bab92e486 100644
--- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.cpp
@@ -34,6 +34,22 @@ ParticleCreationFunc::ParticleCreationFunc (const std::string collision_name,
                 m_num_products_device.push_back(3);
 #endif
             }
+        else if (m_collision_type == CollisionType::DeuteriumTritiumFusion)
+        {
+            m_num_product_species = 2;
+            m_num_products_host.push_back(1);
+            m_num_products_host.push_back(1);
+#ifndef AMREX_USE_GPU
+            // On CPU, the device vector can be filled immediatly
+            m_num_products_device.push_back(1);
+            m_num_products_device.push_back(1);
+#endif
+        }
+        else
+        {
+          amrex::Abort("Unknown collision type in ParticleCreationFunc");
+        }
+
 #ifdef AMREX_USE_GPU
         m_num_products_device.resize(m_num_product_species);
         amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_num_products_host.begin(),