Add option for Barnes thermalisation of gamma-ray energy (#56)

Co-authored-by: Gerrit Leck <gleck@lxbuild07.gsi.de>
Co-authored-by: Luke Shingles <luke.shingles@gmail.com>

artisoptions_christinenonthermal.h

@@ -151,5 +151,7 @@ constexpr bool EXPANSION_OPAC_SAMPLE_KAPPAPLANCK = false;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_classic.h

@@ -147,5 +147,7 @@ constexpr bool EXPANSION_OPAC_SAMPLE_KAPPAPLANCK = false;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
+
 // NOLINTEND(modernize*,misc-unused-parameters)
-#endif  // ARTISOPTIONS_H
+#endif  // ARTISOPTIONS_H

artisoptions_kilonova_lte.h

@@ -146,5 +146,7 @@ constexpr bool EXPANSION_OPAC_SAMPLE_KAPPAPLANCK = false;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_nltenebular.h

@@ -158,5 +158,7 @@ constexpr bool EXPANSION_OPAC_SAMPLE_KAPPAPLANCK = false;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_nltewithoutnonthermal.h

@@ -150,5 +150,7 @@ constexpr bool EXPANSION_OPAC_SAMPLE_KAPPAPLANCK = false;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

constants.h

@@ -67,4 +67,6 @@ enum timestepsizemethods {
   TIMESTEP_SIZES_CONSTANT_THEN_LOGARITHMIC = 3,
 };
 
-#endif
+enum class ThermalisationScheme { DETAILED, BARNES_GLOBAL, BARNES_LOCAL };
+
+#endif

gammapkt.cc

@@ -823,7 +823,7 @@ void pair_prod(Packet &pkt) {
   }
 }
 
-void do_gamma(Packet &pkt, double t2)
+void transport_gamma(Packet &pkt, double t2)
 // Now routine for moving a gamma packet. Idea is that we have as input
 // a gamma packet with known properties at time t1 and we want to follow it
 // until time t2.
@@ -969,4 +969,58 @@ void do_gamma(Packet &pkt, double t2)
   }
 }
 
-}  // namespace gammapkt
+void barnes_thermalization(Packet &pkt, bool local)
+// Barnes treatment: packet is either getting absorbed immediately and locally
+// creating a k-packet or it escapes. The absorption probability matches the
+// Barnes thermalization efficiency, for expressions see the original paper:
+// https://ui.adsabs.harvard.edu/abs/2016ApJ...829..110B
+{
+  // compute thermalization efficiency (= absorption probability)
+  constexpr double mean_gamma_opac = 0.1;
+
+  // determine average initial density
+  double V_0 = 0.;
+  for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) {
+    const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);
+    V_0 += grid::get_modelcell_assocvolume_tmin(mgi);
+  }
+  double rho_0 = 0.;
+  if (!local) {
+    rho_0 = grid::mtot_input / V_0;
+  } else {
+    rho_0 = grid::get_rho_tmin(grid::get_cell_modelgridindex(pkt.where));
+  }
+
+  const double R_0 = pow(3 * V_0 / (4 * PI), 1 / 3.);
+  const double t_0 = grid::get_t_model();
+  const double tau_ineff = sqrt(rho_0 * R_0 * pow(t_0, 2) * mean_gamma_opac);
+  // get current time
+  const double t = t_0 + pkt.prop_time;
+  const double tau = pow(tau_ineff / t, 2.);
+  const double f_gamma = 1. - exp(-tau);
+
+  // either absorb packet or let it escape
+  if (rng_uniform() < f_gamma) {
+    // packet is absorbed and contributes to the heating as a k-packet
+    pkt.type = TYPE_NTLEPTON;
+    pkt.absorptiontype = -4;
+  } else {
+    // let packet escape, i.e. make it inactive
+    pkt.type = TYPE_ESCAPE;
+    grid::change_cell(pkt, -99);
+  }
+}
+
+void do_gamma(Packet &pkt, double t2) {
+  if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::DETAILED) {
+    transport_gamma(pkt, t2);
+  } else if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::BARNES_GLOBAL) {
+    barnes_thermalization(pkt, false);
+  } else if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::BARNES_LOCAL) {
+    barnes_thermalization(pkt, true);
+  } else {
+    __builtin_unreachable();
+  }
+}
+
+}  // namespace gammapkt

grid.cc

@@ -55,7 +55,6 @@ std::array<int, 3> ncoord_model{0};  // the model.txt input grid dimensions
 
 double min_den;  // minimum model density
 
-double mtot_input;
 double mfeg;  /// Total mass of Fe group elements in ejecta
 
 int first_cellindex = -1;  // auto-dermine first cell index in model.txt (usually 1 or 0)

grid.h

@@ -75,6 +75,8 @@ inline ModelGridCell *modelgrid{};
 
 inline int ngrid{0};
 
+inline double mtot_input;
+
 [[nodiscard]] auto get_elements_uppermost_ion(int modelgridindex, int element) -> int;
 void set_elements_uppermost_ion(int modelgridindex, int element, int newvalue);
 [[nodiscard]] auto wid_init(int cellindex, int axis) -> double;