Implement D+D and D+He fusion (#3257)

* Implement D+D -> n+He3 fusion

* Fix logic for fusion reaction

* Check products in a different place

* Correct compilation error

* Implement D+D -> T+p cross-section

* Update example

* Use clearer naming convention for fusion types

* Revert changes to example input script

* Add analysis script

* Progress on tests

* Use 2 species in test

* Correct momentum of colliding species

* Update test

* Update test

* Generalize species names in fusion tests

* Update benchmarks

* Correct typo

* Updated scripts

* Update script so that it works for D+T and D+D

* Update CI

* Add benchmark file

* Correct typo

* Fix compilation on GPU

* Update RZ CI test

* Implement Deuterium-Helium reaction

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

1 files changed, 67 insertions, 6 deletions

diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
index 4bbdb4bc4..81acecd3e 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp
@@ -25,20 +25,75 @@ namespace BinaryCollisionUtils{
             pp_collision_name.getarr("species", species_names);
             auto& species1 = mypc->GetParticleContainerFromName(species_names[0]);
             auto& species2 = mypc->GetParticleContainerFromName(species_names[1]);
+            amrex::Vector<std::string> product_species_name;
+            pp_collision_name.getarr("product_species", product_species_name);
 
             if ((species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::hydrogen3>())
                 ||
                 (species1.AmIA<PhysicalSpecies::hydrogen3>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
                 )
             {
-                return 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");
+                return NuclearFusionType::DeuteriumTritiumToNeutronHelium;
+            }
+            else if (species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
+            {
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 2u,
+                    "ERROR: Deuterium-deuterium 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]);
+                if (
+                    (product_species1.AmIA<PhysicalSpecies::helium3>() && product_species2.AmIA<PhysicalSpecies::neutron>())
+                  ||(product_species1.AmIA<PhysicalSpecies::neutron>() && product_species2.AmIA<PhysicalSpecies::helium3>())){
+                    return NuclearFusionType::DeuteriumDeuteriumToNeutronHelium;
+                } else if (
+                    (product_species1.AmIA<PhysicalSpecies::hydrogen3>() && product_species2.AmIA<PhysicalSpecies::proton>())
+                  ||(product_species1.AmIA<PhysicalSpecies::proton>() && product_species2.AmIA<PhysicalSpecies::hydrogen3>())){
+                    return NuclearFusionType::DeuteriumDeuteriumToProtonTritium;
+                } else {
+                    amrex::Abort("ERROR: Product species of proton-boron fusion must be of type helium3 and neutron, or tritium and proton");
+                }
+            }
+            else if ((species1.AmIA<PhysicalSpecies::hydrogen2>() && species2.AmIA<PhysicalSpecies::helium3>())
+                ||
+                (species1.AmIA<PhysicalSpecies::helium3>() && species2.AmIA<PhysicalSpecies::hydrogen2>())
+                )
+            {
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 2u,
+                    "ERROR: Deuterium-helium 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::proton>())
+                    ||
+                    (product_species1.AmIA<PhysicalSpecies::proton>() && product_species2.AmIA<PhysicalSpecies::helium4>()),
+                    "ERROR: Product species of deuterium-helium fusion must be of type proton and helium4");
+                return NuclearFusionType::DeuteriumHeliumToProtonHelium;
             }
             else if ((species1.AmIA<PhysicalSpecies::proton>() && species2.AmIA<PhysicalSpecies::boron11>())
                 ||
                 (species1.AmIA<PhysicalSpecies::boron11>() && species2.AmIA<PhysicalSpecies::proton>())
                 )
             {
-                return NuclearFusionType::ProtonBoron;
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species_name.size() == 1,
+                    "ERROR: Proton-boron must contain exactly one product species");
+                auto& product_species = mypc->GetParticleContainerFromName(product_species_name[0]);
+                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                    product_species.AmIA<PhysicalSpecies::helium>(),
+                    "ERROR: Product species of proton-boron fusion must be of type alpha");
+                return NuclearFusionType::ProtonBoronToAlphas;
             }
             amrex::Abort("Binary nuclear fusion not implemented between species " +
                          species_names[0] + " of type " + species1.getSpeciesTypeName() +
@@ -63,10 +118,16 @@ 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;
+            if (fusion_type == NuclearFusionType::DeuteriumTritiumToNeutronHelium)
+                return CollisionType::DeuteriumTritiumToNeutronHeliumFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumDeuteriumToProtonTritium)
+                return CollisionType::DeuteriumDeuteriumToProtonTritiumFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumDeuteriumToNeutronHelium)
+                return CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion;
+            if (fusion_type == NuclearFusionType::DeuteriumHeliumToProtonHelium)
+                return CollisionType::DeuteriumHeliumToProtonHeliumFusion;
+            if (fusion_type == NuclearFusionType::ProtonBoronToAlphas)
+                return CollisionType::ProtonBoronToAlphasFusion;
             amrex::Abort("Invalid nuclear fusion type");
             return CollisionType::Undefined;
         }