Add fusion particle production diagnostic

Docs/source/usage/parameters.rst

@@ -3064,6 +3064,23 @@ Details about the collision models can be found in the :ref:`theory section <mul
     With ``forward``, the scattering angle is set to zero, i.e. the products are emitted in the same direction as the reactant (in the center of mass frame).
     With ``backward``, the scattering angle is set to :math:`\pi`, i.e. the products are emitted in the opposite direction of the reactant (in the center of mass frame).
 
+.. pp:param:: <collision_name>.create_products
+    :type: ``bool``
+    :default: ``1``
+    :optional:
+
+    Only for ``nuclearfusion``. When true, the product particles are created, otherwise not.
+
+.. pp:param:: <collision_name>.save_particle_production
+    :type: ``bool``
+    :default: ``0``
+    :optional:
+
+    Only for ``nuclearfusion``.
+    When true, the integrated product particle density is saved in a MultiFab with the name ``<collision_name>_particle_production``.
+    The data can be written out by adding that name to the ``<diag_name>.fields_to_plot`` input parameter.
+    The option can be used in conjunction with ``<collision_name>.create_products`` to save only the product density and not create particles.
+
 .. pp:param:: <collision_name>.background_density
     :type: ``float``
 

Examples/Tests/nuclear_fusion/analysis_two_product_fusion.py

@@ -65,13 +65,15 @@
     reaction_type = "DT"
     reactant_species = ["deuterium", "tritium"]
     product_species = ["helium4", "neutron"]
+    collision_names = ["DTF1", "DTF2"]
     ntests = 2
     E_fusion = 17.58929696 * MeV_to_Joule
 else:
     # else, this is the D+D test
     reaction_type = "DD"
     reactant_species = ["deuterium", "hydrogen2"]
     product_species = ["helium3", "neutron"]
+    collision_names = ["DDNHeF1"]
     ntests = 1
     E_fusion = 3.26891111e6 * MeV_to_Joule
 
@@ -399,6 +401,19 @@ def check_macroparticle_number(
         atol=5.0 * std_macroparticle_number,
     )
 
+    if "particle_production" in data:
+        w_sum = data[product_species[0] + "_w_end"].sum()
+        n_sum = data["particle_production"].sum()
+        tolerance = 0.02
+        print(
+            f"Check particle production diagnostic for collision {data['collision_name']}:"
+        )
+        print(f"from particles  = {w_sum}")
+        print(f"from diagnostic = {n_sum}")
+        print(f"error = {np.abs(w_sum - n_sum) / w_sum}")
+        print(f"tolerance = {tolerance}")
+        assert is_close(w_sum, n_sum, rtol=tolerance)
+
     ## used in subsequent function
     return expected_fusion_number
 
@@ -545,6 +560,13 @@ def main():
         # General checks that are performed for all tests
         generic_check(data)
 
+        product_production_name = f"{collision_names[i - 1]}_particle_production"
+        if ("boxlib", product_production_name) in ds_end.field_list:
+            data["collision_name"] = collision_names[i - 1]
+            data["particle_production"] = field_data_end[
+                "boxlib", product_production_name
+            ].to_ndarray()
+
         # Checks that are specific to test number i
         eval("specific_check" + str(i) + "(data, dt)")
 

Examples/Tests/nuclear_fusion/inputs_test_3d_deuterium_tritium_fusion

@@ -2,7 +2,9 @@
 ####### GENERAL PARAMETERS ######
 #################################
 ## With these parameters, each cell has a size of exactly 1 by 1 by 1
-max_step = 1
+# Run two steps. The first step to check the fusion collisions, the
+# second step to provide a comparison for the rerun.
+max_step = 2
 amr.n_cell = 8 8 16
 amr.max_grid_size = 8
 amr.blocking_factor = 8
@@ -121,15 +123,24 @@ DTF1.species = deuterium_1 tritium_1
 DTF1.product_species = helium4_1 neutron_1
 DTF1.type = nuclearfusion
 DTF1.event_multiplier = 1.e50
+DTF1.create_products = 1
+DTF1.save_particle_production = 1
 
 DTF2.species = deuterium_2 tritium_2
 DTF2.product_species = helium4_2 neutron_2
 DTF2.type = nuclearfusion
 DTF2.event_multiplier = 1.e15
 DTF2.probability_target_value = 0.02
+DTF2.create_products = 1
+DTF2.save_particle_production = 1
 
 # Diagnostics
-diagnostics.diags_names = diag1
+diagnostics.diags_names = diag1 checkpoint
 diag1.intervals = 1
 diag1.diag_type = Full
-diag1.fields_to_plot = rho
+diag1.fields_to_plot = rho DTF1_particle_production DTF2_particle_production
+
+checkpoint.format = checkpoint
+checkpoint.diag_type = Full
+checkpoint.intervals = 1:1:1
+checkpoint.dump_last_timestep = 0

Regression/Checksum/benchmarks_json/test_3d_deuterium_tritium_fusion.json

@@ -1,24 +1,12 @@
 {
-  "lev=0": {
-    "rho": 0.0
-  },
-  "neutron_2": {
-    "particle_momentum_x": 1.5369063360838545e-15,
-    "particle_momentum_y": 1.5327717119671177e-15,
-    "particle_momentum_z": 1.5632203763888702e-15,
-    "particle_position_x": 136756.17264787608,
-    "particle_position_y": 136453.48037878488,
-    "particle_position_z": 290503.22456411575,
-    "particle_weight": 6.347081228434342e+18
-  },
-  "neutron_1": {
-    "particle_momentum_x": 1.7270063957926637e-15,
-    "particle_momentum_y": 1.7295255445271788e-15,
-    "particle_momentum_z": 1.7442619148907942e-15,
-    "particle_position_x": 151546.7150677448,
-    "particle_position_y": 151695.5086129642,
-    "particle_position_z": 323004.4593236664,
-    "particle_weight": 4.337788155202713e-28
+  "deuterium_1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.8872569136407634e-13,
+    "particle_position_x": 40958427.50992301,
+    "particle_position_y": 40959476.34450768,
+    "particle_position_z": 81921930.27522022,
+    "particle_weight": 1024.0000000000002
   },
   "deuterium_2": {
     "particle_momentum_x": 0.0,
@@ -29,23 +17,14 @@
     "particle_position_z": 8192362.405430986,
     "particle_weight": 1.0240001137714307e+30
   },
-  "tritium_1": {
-    "particle_momentum_x": 0.0,
-    "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.887256913640763e-13,
-    "particle_position_x": 40959200.11081588,
-    "particle_position_y": 40960650.407891415,
-    "particle_position_z": 81920772.7986121,
-    "particle_weight": 1024.0000000000002
-  },
-  "deuterium_1": {
-    "particle_momentum_x": 0.0,
-    "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.8872569136407634e-13,
-    "particle_position_x": 40958427.50992301,
-    "particle_position_y": 40959476.34450768,
-    "particle_position_z": 81921930.27522022,
-    "particle_weight": 1024.0000000000002
+  "helium4_1": {
+    "particle_momentum_x": 1.7270063957926637e-15,
+    "particle_momentum_y": 1.7295255445271788e-15,
+    "particle_momentum_z": 1.7442619148907942e-15,
+    "particle_position_x": 151546.7150677448,
+    "particle_position_y": 151695.5086129642,
+    "particle_position_z": 323004.4593236664,
+    "particle_weight": 4.337788155202713e-28
   },
   "helium4_2": {
     "particle_momentum_x": 1.5369063360838545e-15,
@@ -56,7 +35,12 @@
     "particle_position_z": 290503.22456411575,
     "particle_weight": 6.347081228434342e+18
   },
-  "helium4_1": {
+  "lev=0": {
+    "DTF1_particle_production": 4.415241048343040e-28,
+    "DTF2_particle_production": 6.295151973703336e+18,
+    "rho": 0.0
+  },
+  "neutron_1": {
     "particle_momentum_x": 1.7270063957926637e-15,
     "particle_momentum_y": 1.7295255445271788e-15,
     "particle_momentum_z": 1.7442619148907942e-15,
@@ -65,6 +49,24 @@
     "particle_position_z": 323004.4593236664,
     "particle_weight": 4.337788155202713e-28
   },
+  "neutron_2": {
+    "particle_momentum_x": 1.5369063360838545e-15,
+    "particle_momentum_y": 1.5327717119671177e-15,
+    "particle_momentum_z": 1.5632203763888702e-15,
+    "particle_position_x": 136756.17264787608,
+    "particle_position_y": 136453.48037878488,
+    "particle_position_z": 290503.22456411575,
+    "particle_weight": 6.347081228434342e+18
+  },
+  "tritium_1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.887256913640763e-13,
+    "particle_position_x": 40959200.11081588,
+    "particle_position_y": 40960650.407891415,
+    "particle_position_z": 81920772.7986121,
+    "particle_weight": 1024.0000000000002
+  },
   "tritium_2": {
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
@@ -74,4 +76,4 @@
     "particle_position_z": 819126.8984535292,
     "particle_weight": 1.0239999999365294e+29
   }
-}
+}

Source/Particles/Collision/BinaryCollision/BinaryCollision.H

@@ -177,6 +177,11 @@ public:
     BinaryCollision ( BinaryCollision&& )                  = delete;
     BinaryCollision& operator= ( BinaryCollision&& )       = delete;
 
+    void AllocData () override
+    {
+        m_binary_collision_functor.AllocData();
+    }
+
     /** Perform the collisions
      *
      * @param cur_time Current time
@@ -296,7 +301,7 @@ public:
 
         ABLASTR_PROFILE("BinaryCollision::doCollisionsWithinTile");
 
-        const auto& binary_collision_functor = m_binary_collision_functor.executor();
+        const auto& binary_collision_functor = m_binary_collision_functor.executor(lev, mfi);
         const bool have_product_species = m_have_product_species;
 
         // Store product species data in vectors
@@ -326,6 +331,8 @@ public:
             global_debye_length_data = global_debye_length_fab.dataPtr();
         }
 
+        const amrex::Box tilebox = mfi.tilebox();
+
         amrex::Geometry const& geom_lev = WarpX::GetInstance().Geom(lev);
         // dV is level-specific: cell volume at this refinement level (smaller on fine levels).
         amrex::ParticleReal const dV = AMREX_D_TERM(geom_lev.CellSize(0), *geom_lev.CellSize(1), *geom_lev.CellSize(2));
@@ -652,6 +659,7 @@ public:
 
                     // Use a bisection algorithm to find the index of the cell in which this pair is located
                     const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
+                    const amrex::IntVect global_index = tilebox.atOffset(i_cell);
 
                     // The particles from species1 that are in the cell `i_cell` are
                     // given by the `indices_1[cell_start_1:cell_stop_1]`
@@ -691,7 +699,7 @@ public:
                         soa_1, soa_1, get_position_1, get_position_1,
                         n1, n1, T1, T1, global_lamdb,
                         q1, q1, m1, m1, dt, dV*volume_factor(i_cell), coll_idx,
-                        cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
+                        cell_start_pair, global_index, p_mask, p_pair_indices_1, p_pair_indices_2,
                         p_pair_reaction_weight, p_product_data, engine);
                 }
             );
@@ -1238,6 +1246,7 @@ public:
 
                     // Use a bisection algorithm to find the index of the cell in which this pair is located
                     const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
+                    const amrex::IntVect global_index = tilebox.atOffset(i_cell);
 
                     // The particles from species1 that are in the cell `i_cell` are
                     // given by the `indices_1[cell_start_1:cell_stop_1]`
@@ -1284,7 +1293,7 @@ public:
                         soa_1, soa_2, get_position_1, get_position_2,
                         n1, n2, T1, T2, global_lamdb,
                         q1, q2, m1, m2, dt, dV*volume_factor(i_cell), coll_idx,
-                        cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
+                        cell_start_pair, global_index, p_mask, p_pair_indices_1, p_pair_indices_2,
                         p_pair_reaction_weight, p_product_data, engine);
                 }
             );

Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H

@@ -8,6 +8,7 @@
 #ifndef WARPX_BREMSSTRAHLUNG_FUNC_H_
 #define WARPX_BREMSSTRAHLUNG_FUNC_H_
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/MultiParticleContainer.H"
@@ -31,6 +32,7 @@
  *  effectively created in the particle creation functor.
  */
 class BremsstrahlungFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -92,7 +94,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const m1, amrex::ParticleReal const m2,
             amrex::Real const dt, amrex::Real const /*dV*/, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const /*global_index*/,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* AMREX_RESTRICT p_product_data,
@@ -406,7 +409,7 @@ public:
 
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     [[nodiscard]] bool use_global_debye_length() const { return m_use_global_debye_length; }
 

Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H

@@ -8,6 +8,7 @@
 #ifndef WARPX_PAIRWISE_COULOMB_COLLISION_FUNC_H_
 #define WARPX_PAIRWISE_COULOMB_COLLISION_FUNC_H_
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "ElasticCollisionPerez.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/WarpXParticleContainer.H"
@@ -24,6 +25,7 @@
  *        ElasticCollisionPerez. It also reads and contains the Coulomb logarithm.
  */
 class PairWiseCoulombCollisionFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -120,7 +122,8 @@ public:
             amrex::ParticleReal const  q1, amrex::ParticleReal const  q2,
             amrex::ParticleReal const  m1, amrex::ParticleReal const  m2,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const /*cell_start_pair*/, index_type* /*p_mask*/,
+            index_type const /*cell_start_pair*/, amrex::IntVect const /*global_index*/,
+            index_type* /*p_mask*/,
             index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/,
             amrex::ParticleReal* /*p_pair_reaction_weight*/,
             amrex::ParticleReal* /*p_product_data*/,
@@ -144,7 +147,7 @@ public:
         bool m_isSameSpecies;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     [[nodiscard]] bool use_global_debye_length() const { return m_use_global_debye_length; }
 

Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H

@@ -12,6 +12,7 @@
 
 #include "CollisionFilterFunc.H"
 
+#include "Particles/Collision/CollisionFuncBase.H"
 #include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
 #include "Particles/Collision/BinaryCollision/ParticleShufflers.H"
 #include "Particles/Collision/CollisionBase.H"
@@ -35,6 +36,7 @@
  * used for binary Coulomb collisions and the nuclear fusion module.
  */
 class DSMCFunc
+    : public CollisionFuncBase
 {
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
@@ -103,7 +105,8 @@ public:
             amrex::ParticleReal const  /*q1*/, amrex::ParticleReal const  /*q2*/,
             amrex::ParticleReal const  m1, amrex::ParticleReal const  m2,
             amrex::Real const  dt, amrex::Real const dV, index_type coll_idx,
-            index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
+            index_type const cell_start_pair, amrex::IntVect const /*global_index*/,
+            index_type* AMREX_RESTRICT p_mask,
             index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
             amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
             amrex::ParticleReal* /*p_product_data*/,
@@ -183,7 +186,7 @@ public:
         ScatteringProcess::Executor* m_scattering_processes_data;
     };
 
-    [[nodiscard]] Executor const& executor () const { return m_exe; }
+    [[nodiscard]] Executor const& executor (int const /*lev*/, amrex::MFIter const& /*mfi*/) const { return m_exe; }
 
     bool use_global_debye_length() { return false; }