Add thermalisation scheme DETAILEDWITHGAMMAPRODUCTS

constants.h

@@ -59,6 +59,6 @@ constexpr size_t GSLWSIZE = 16384;  // GSL integration workspace size
 
 enum class TimeStepSizeMethod { LOGARITHMIC, CONSTANT, LOGARITHMIC_THEN_CONSTANT, CONSTANT_THEN_LOGARITHMIC };
 
-enum class ThermalisationScheme { INSTANT, DETAILED, BARNES, WOLLAEGER, GUTTMAN };
+enum class ThermalisationScheme { INSTANT, DETAILED, DETAILEDWITHGAMMAPRODUCTS, BARNES, WOLLAEGER, GUTTMAN };
 
 #endif

gammapkt.cc

@@ -466,8 +466,13 @@ void compton_scatter(Packet &pkt) {
 
     pkt.last_cross = BOUNDARY_NONE;  // allow it to re-cross a boundary
   } else {
-    // It's converted to an e-minus packet.
-    pkt.type = TYPE_NTLEPTON_DEPOSITED;
+    // energy loss of the gamma becomes energy of the electron (needed to calculate time-dependent thermalisation rate)
+    if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
+      pkt.nu_cmf = pkt.nu_cmf * (1 - 1 / f);
+      pkt.type = TYPE_NONTHERMAL_PREDEPOSIT_BETAMINUS;
+    } else {
+      pkt.type = TYPE_NTLEPTON_DEPOSITED;
+    }
     pkt.absorptiontype = -3;
     stats::increment(stats::COUNTER_NT_STAT_FROM_GAMMA);
   }
@@ -657,6 +662,9 @@ void update_gamma_dep(const Packet &pkt, const double dist, const int mgi, const
   if (!(dist > 0)) {
     return;
   }
+  if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
+    return;  // don't instantly deposit energy from gamma rays, handle the particles they produce instead
+  }
 
   const double doppler_sq = doppler_squared_nucmf_on_nurf(pkt.pos, pkt.dir, pkt.prop_time);
 
@@ -693,11 +701,16 @@ void pair_prod(Packet &pkt) {
     printout("prob_gamma < 0. pair_prod. Abort. %g\n", prob_gamma);
     std::abort();
   }
+  const auto zrand = rng_uniform();
+  if (zrand > prob_gamma) {
+    if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
+      // Convert it to an e-minus or positron kinetic energy packet
+      pkt.type = (rng_uniform() > 0.5) ? TYPE_NONTHERMAL_PREDEPOSIT_BETAMINUS : TYPE_NONTHERMAL_PREDEPOSIT_BETAPLUS;
+    } else {
+      pkt.type = TYPE_NTLEPTON_DEPOSITED;
+    }
 
-  if (rng_uniform() > prob_gamma) {
-    // Convert it to an e-minus packet - actually it could be positron EK too, but this works
-    // for consistency with compton_scatter.
-    pkt.type = TYPE_NTLEPTON_DEPOSITED;
+    // nu_cmf stays the same as the gamma energy becomes the kinetic energy of the electron
     pkt.absorptiontype = -5;
     stats::increment(stats::COUNTER_NT_STAT_FROM_GAMMA);
   } else {
@@ -837,29 +850,23 @@ void transport_gamma(Packet &pkt, const double t2) {
       // Compton scattering.
       compton_scatter(pkt);
     } else if ((chi_compton + chi_photo_electric) > chi_rnd) {
-      // Photo electric effect - makes it a k-packet for sure.
-      pkt.type = TYPE_NTLEPTON_DEPOSITED;
+      // Photo electric effect
+      if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
+        pkt.type = TYPE_NONTHERMAL_PREDEPOSIT_BETAMINUS;
+        // nu_cmf stays the same as the gamma-ray energy becomes the kinetic energy of the electron (minus ionisation
+        // energy but this is neglected here)
+      } else {
+        pkt.type = TYPE_NTLEPTON_DEPOSITED;
+      }
+
       pkt.absorptiontype = -4;
       stats::increment(stats::COUNTER_NT_STAT_FROM_GAMMA);
-    } else if ((chi_compton + chi_photo_electric + chi_pair_prod) > chi_rnd) {
+    } else {
       // It's a pair production
       pair_prod(pkt);
-    } else {
-      printout(
-          "Failed to identify event. Gamma (1). chi_compton %g chi_photo_electric %g chi_tot %g chi_rnd %g Abort.\n",
-          chi_compton, chi_photo_electric, chi_tot, chi_rnd);
-      const int cellindex = pkt.where;
-      printout(
-          " globals::cell[pkt.where].rho %g pkt.nu_cmf %g pkt.dir[0] %g pkt.dir[1] %g "
-          "pkt.dir[2] %g pkt.pos[0] %g pkt.pos[1] %g pkt.pos[2] %g \n",
-          grid::get_rho(grid::get_cell_modelgridindex(cellindex)), pkt.nu_cmf, pkt.dir[0], pkt.dir[0], pkt.dir[1],
-          pkt.dir[2], pkt.pos[1], pkt.pos[2]);
-
-      std::abort();
     }
   } else {
-    printout("Failed to identify event. Gamma (2). edist %g, sdist %g, tdist %g Abort.\n", edist, sdist, tdist);
-    std::abort();
+    assert_always(false);
   }
 }
 
@@ -1095,9 +1102,12 @@ __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) {
   }
 
   if (pkt.type != TYPE_GAMMA && pkt.type != TYPE_ESCAPE) {
-    atomicadd(globals::timesteps[nts].gamma_dep_discrete, pkt.e_cmf);
+    if constexpr (PARTICLE_THERMALISATION_SCHEME != ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
+      atomicadd(globals::timesteps[nts].gamma_dep_discrete, pkt.e_cmf);
+    }
 
-    if constexpr (GAMMA_THERMALISATION_SCHEME != ThermalisationScheme::DETAILED) {
+    if constexpr (GAMMA_THERMALISATION_SCHEME != ThermalisationScheme::DETAILED &&
+                  GAMMA_THERMALISATION_SCHEME != ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
       // no transport, so the path-based gamma deposition estimator won't get updated unless we do it here
       const int mgi = grid::get_cell_modelgridindex(pkt.where);
       const int nonemptymgi = grid::get_modelcell_nonemptymgi(mgi);

update_packets.cc

@@ -71,6 +71,7 @@ void do_nonthermal_predeposit(Packet &pkt, const int nts, const double t2) {
       grid::change_cell(pkt, -99);
     }
   } else {
+    // ThermalisationScheme::DETAILED or ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS
     // local, detailed absorption following Shingles+2023
     const double rho = grid::get_rho(grid::get_cell_modelgridindex(pkt.where));
 
@@ -109,20 +110,29 @@ void do_nonthermal_predeposit(Packet &pkt, const int nts, const double t2) {
     pkt.prop_time = t_new;
   }
 
-  if (priortype == TYPE_NONTHERMAL_PREDEPOSIT_BETAMINUS) {
-    atomicadd(globals::dep_estimator_electron[nonemptymgi], en_deposited);
-    if (pkt.type == TYPE_NTLEPTON_DEPOSITED) {
-      atomicadd(globals::timesteps[nts].electron_dep_discrete, pkt.e_cmf);
-    }
-  } else if (priortype == TYPE_NONTHERMAL_PREDEPOSIT_BETAPLUS) {
-    atomicadd(globals::dep_estimator_positron[nonemptymgi], en_deposited);
-    if (pkt.type == TYPE_NTLEPTON_DEPOSITED) {
-      atomicadd(globals::timesteps[nts].positron_dep_discrete, pkt.e_cmf);
+  // for DETAILEDWITHGAMMAPRODUCTS, gamma-ray deposition will lead to predeposit beta particles, but they will count
+  // toward "gamma deposition" not particle deposition
+  if (pkt.originated_from_particlenotgamma) {
+    if (priortype == TYPE_NONTHERMAL_PREDEPOSIT_BETAMINUS) {
+      atomicadd(globals::dep_estimator_electron[nonemptymgi], en_deposited);
+      if (pkt.type == TYPE_NTLEPTON_DEPOSITED) {
+        atomicadd(globals::timesteps[nts].electron_dep_discrete, pkt.e_cmf);
+      }
+    } else if (priortype == TYPE_NONTHERMAL_PREDEPOSIT_BETAPLUS) {
+      atomicadd(globals::dep_estimator_positron[nonemptymgi], en_deposited);
+      if (pkt.type == TYPE_NTLEPTON_DEPOSITED) {
+        atomicadd(globals::timesteps[nts].positron_dep_discrete, pkt.e_cmf);
+      }
+    } else if (priortype == TYPE_NONTHERMAL_PREDEPOSIT_ALPHA) {
+      atomicadd(globals::dep_estimator_alpha[nonemptymgi], en_deposited);
+      if (pkt.type == TYPE_NTLEPTON_DEPOSITED) {
+        atomicadd(globals::timesteps[nts].alpha_dep_discrete, pkt.e_cmf);
+      }
     }
-  } else if (priortype == TYPE_NONTHERMAL_PREDEPOSIT_ALPHA) {
-    atomicadd(globals::dep_estimator_alpha[nonemptymgi], en_deposited);
+  } else if constexpr (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::DETAILEDWITHGAMMAPRODUCTS) {
+    atomicadd(globals::dep_estimator_gamma[nonemptymgi], en_deposited);
     if (pkt.type == TYPE_NTLEPTON_DEPOSITED) {
-      atomicadd(globals::timesteps[nts].alpha_dep_discrete, pkt.e_cmf);
+      atomicadd(globals::timesteps[nts].gamma_dep_discrete, pkt.e_cmf);
     }
   }
 }