D-T fusion

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{
 

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::hydrogen>() && 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::hydrogen>() && species2.AmIA<PhysicalSpecies::boron11>())
                 ||
                 (species1.AmIA<PhysicalSpecies::boron11>() && species2.AmIA<PhysicalSpecies::hydrogen>())
                 )
@@ -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");

Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumCrossSection.H

@@ -0,0 +1,34 @@
+/* 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.
+ *
+ * @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;
+
+    // Convert cross section to SI units: barn to square meter
+    constexpr auto barn_to_sqm = 1.e-28_prt;
+    return 0.*barn_to_sqm;
+}
+
+#endif // DEUTERIUM_TRITIUM_TRITIUM_FUSION_CROSS_SECTION_H

Source/Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumInitializeMomentum.H

@@ -0,0 +1,64 @@
+/* Copyright 2022 Remi Lehe
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef DEUTERIUM_TRITIUM_FUSION_INITIALIZE_MOMENTUM_H
+#define DEUTERIUM_TRITIUM_FUSION_INITIALIZE_MOMENTUM_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 particles produced from
+     * deuterium-tritium fusion.
+     *
+     * @param[in] soa_1 struct of array data of the first colliding species (can be either deuterium
+     * or tritium)
+     * @param[in] soa_2 struct of array data of the second colliding species (can be either deuterium
+     * or tritium)
+     * @param[out] ...
+     * @param[in] idx_1 index of first colliding macroparticle
+     * @param[in] idx_2 index of second colliding macroparticle
+     * @param[in]...
+     * @param[in] m1 mass of first colliding species
+     * @param[in] m2 mass of second colliding species
+     * @param[in] engine the random engine
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    void DeuteriumTritiumFusionInitializeMomentum (
+                            const SoaData_type& soa_1, const SoaData_type& soa_2,
+                            SoaData_type& soa_alpha,
+                            const index_type& idx_1, const index_type& idx_2,
+                            const index_type& idx_alpha_start,
+                            const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
+                            const amrex::RandomEngine& engine)
+    {
+        // General notations in this function:
+        //     x_sq denotes the square of x
+        //     x_star denotes the value of x in the proton+boron center of mass frame
+        //     x_Bestar denotes the value of x in the Beryllium rest frame
+
+        using namespace amrex::literals;
+    }
+
+}
+
+#endif // DEUTERIUM_TRITIUM_FUSION_INITIALIZE_MOMENTUM_H

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() == 2,
+                "ERROR: Deuterium-tritium 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_species.AmIA<PhysicalSpecies::helium4>() && product_species.AmIA<PhysicalSpecies::neutron>())
+                ||
+                (product_species.AmIA<PhysicalSpecies::neutron>() && product_species.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(

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);
     }

Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H

@@ -10,6 +10,7 @@
 
 #include "BinaryCollisionUtils.H"
 
+#include "Particles/Collision/BinaryCollision/NuclearFusion/DeuteriumTritiumFusionInitializeMomentum.H"
 #include "Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H"
 #include "Particles/ParticleCreation/SmartCopy.H"
 #include "Particles/MultiParticleContainer.H"
@@ -208,6 +209,10 @@ public:
 
                 // Initialize the product particles' momentum, using a function depending on the
                 // specific collision type
+                if (t_collision_type == CollisionType::DeuteriumTritiumFusion)
+                {
+                    // TODO DeuteriumTritiumFusionInitializeMomentum
+                }
                 if (t_collision_type == CollisionType::ProtonBoronFusion)
                 {
                     const index_type product_start_index = products_np_data[0] + 2*p_offsets[i]*

Source/Particles/SpeciesPhysicalProperties.H

@@ -17,8 +17,9 @@
 #include <map>
 #include <string>
 
-enum struct PhysicalSpecies{unspecified=0, electron, positron, photon, hydrogen, helium, boron,
-                            boron10, boron11, carbon, nitrogen, oxygen, argon, copper, xenon};
+enum struct PhysicalSpecies{unspecified=0, electron, positron, photon, neutron,
+    hydrogen, hydrogen2, hydrogen3, helium, helium3, helium4, boron, boron10,
+    boron11, carbon, nitrogen, oxygen, argon, copper, xenon};
 
 namespace species
 {
@@ -33,12 +34,20 @@ namespace species
             return PhysicalSpecies::positron;
         if( species=="photon" )
             return PhysicalSpecies::photon;
-        if( species=="hydrogen" )
-            return PhysicalSpecies::hydrogen;
-        if( species=="proton" )
+        if( species=="neutron" )
+            return PhysicalSpecies::neutron;
+        if( (species=="hydrogen") || (species=="proton") )
             return PhysicalSpecies::hydrogen;
+        if( (species=="hydrogen2") || (species=="deuterium") )
+            return PhysicalSpecies::hydrogen2;
+        if (species=="hydrogen3") || (species=="tritium") )
+            return PhysicalSpecies::hydrogen3;
         if( species=="helium" )
             return PhysicalSpecies::helium;
+        if( species=="helium3" )
+            return PhysicalSpecies::helium3;
+        if( species=="helium4" )
+            return PhysicalSpecies::helium4;
         if( species=="alpha" )
             return PhysicalSpecies::helium;
         if( species=="boron" )
@@ -75,10 +84,20 @@ namespace species
             return PhysConst::q_e;
         case PhysicalSpecies::photon:
             return 0.;
+        case PhysicalSpecies::neutron:
+            return 0.;
         case PhysicalSpecies::hydrogen:
             return PhysConst::q_e;
+        case PhysicalSpecies::hydrogen2:
+            return PhysConst::q_e;
+        case PhysicalSpecies::hydrogen3:
+            return PhysConst::q_e;
         case PhysicalSpecies::helium:
             return PhysConst::q_e * amrex::Real(2.0);
+        case PhysicalSpecies::helium3:
+            return PhysConst::q_e * amrex::Real(2.0);
+        case PhysicalSpecies::helium4:
+            return PhysConst::q_e * amrex::Real(2.0);
         case PhysicalSpecies::boron:
             return PhysConst::q_e * amrex::Real(5.0);
         case PhysicalSpecies::boron10:
@@ -115,10 +134,20 @@ namespace species
             return PhysConst::m_e;
         case PhysicalSpecies::photon:
             return 0.;
+        case PhysicalSpecies::neutron:
+            return PhysConst::m_p * amrex::Real(1.0013784193052508);
         case PhysicalSpecies::hydrogen:
             return PhysConst::m_p;
+        case PhysicalSpecies::hydrogen2:
+            return PhysConst::m_p * amrex::Real(1.99901);
+        case PhysicalSpecies::hydrogen3:
+            return PhysConst::m_p * amrex::Real(2.99372);
         case PhysicalSpecies::helium:
             return PhysConst::m_p * amrex::Real(3.97369);
+        case PhysicalSpecies::helium3:
+            return PhysConst::m_p * amrex::Real(2.99369);
+        case PhysicalSpecies::helium4:
+            return PhysConst::m_p * amrex::Real(3.97314);
         case PhysicalSpecies::boron:
             return PhysConst::m_p * amrex::Real(10.7319);
         case PhysicalSpecies::boron10:
@@ -157,8 +186,16 @@ namespace species
             return "photon";
         case PhysicalSpecies::hydrogen:
             return "hydrogen";
+        case PhysicalSpecies::hydrogen2:
+            return "hydrogen2";
+        case PhysicalSpecies::hydrogen3:
+            return "hydrogen3";
         case PhysicalSpecies::helium:
             return "helium";
+        case PhysicalSpecies::helium3:
+            return "helium3";
+        case PhysicalSpecies::helium4:
+            return "helium4";
         case PhysicalSpecies::boron:
             return "boron";
         case PhysicalSpecies::boron10: