Add Barnes Particle Thermalisation (#76)

Add Barnes thermalisation for alpha and beta particles --------- Co-authored-by: Luke Shingles <luke.shingles@gmail.com>
artis-mcrt · Jun 11, 2024 · c374125 · c374125
1 parent d3807f3
commit c374125
Show file tree

Hide file tree

Showing 16 changed files with 451 additions and 37 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -28,6 +28,7 @@ jobs:
                         kilonova_1d_1dgrid,
                         kilonova_1d_3dgrid,
                         kilonova_2d_2dgrid,
+                        kilonova_2d_2dgrid_barnesthermalisation,
                         kilonova_2d_2dgrid_expansionopac,
                         kilonova_2d_2dgrid_xcomgammaphotoion,
                         kilonova_2d_3dgrid,

diff --git a/artisoptions_christinenonthermal.h b/artisoptions_christinenonthermal.h
@@ -133,8 +133,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false;
 
 constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false;
 
-constexpr bool INSTANT_PARTICLE_DEPOSITION = true;
-
 constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC;
 
 constexpr double FIXED_TIMESTEP_WIDTH = -1.;
@@ -151,6 +149,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT;
+
 constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
 // NOLINTEND(modernize*,misc-unused-parameters)

diff --git a/artisoptions_classic.h b/artisoptions_classic.h
@@ -131,8 +131,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false;
 
 constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false;
 
-constexpr bool INSTANT_PARTICLE_DEPOSITION = true;
-
 constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC;
 
 constexpr double FIXED_TIMESTEP_WIDTH = -1.;
@@ -149,6 +147,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT;
+
 constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
 // NOLINTEND(modernize*,misc-unused-parameters)

diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h
@@ -130,8 +130,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = true;
 
 constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false;
 
-constexpr bool INSTANT_PARTICLE_DEPOSITION = false;
-
 constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC;
 
 constexpr double FIXED_TIMESTEP_WIDTH = -1.;
@@ -148,6 +146,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
+
 constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
 // NOLINTEND(modernize*,misc-unused-parameters)

diff --git a/artisoptions_nltenebular.h b/artisoptions_nltenebular.h
@@ -142,8 +142,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false;
 
 constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false;
 
-constexpr bool INSTANT_PARTICLE_DEPOSITION = true;
-
 constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC;
 
 constexpr double FIXED_TIMESTEP_WIDTH = -1.;
@@ -160,6 +158,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT;
+
 constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
 // NOLINTEND(modernize*,misc-unused-parameters)

diff --git a/artisoptions_nltewithoutnonthermal.h b/artisoptions_nltewithoutnonthermal.h
@@ -134,8 +134,6 @@ constexpr bool USE_CALCULATED_MEANATOMICWEIGHT = false;
 
 constexpr bool WRITE_PARTIAL_EMISSIONABSORPTIONSPEC = false;
 
-constexpr bool INSTANT_PARTICLE_DEPOSITION = true;
-
 constexpr enum timestepsizemethods TIMESTEP_SIZE_METHOD = TIMESTEP_SIZES_LOGARITHMIC;
 
 constexpr double FIXED_TIMESTEP_WIDTH = 0.1;
@@ -152,6 +150,8 @@ constexpr float RPKT_BOUNDBOUND_THERMALISATION_PROBABILITY = -1.;
 
 constexpr bool USE_XCOM_GAMMAPHOTOION = false;
 
+constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::INSTANT;
+
 constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED;
 
 // NOLINTEND(modernize*,misc-unused-parameters)

diff --git a/constants.h b/constants.h
@@ -67,6 +67,6 @@ enum timestepsizemethods {
   TIMESTEP_SIZES_CONSTANT_THEN_LOGARITHMIC = 3,
 };
 
-enum class ThermalisationScheme { DETAILED, BARNES_GLOBAL, BARNES_LOCAL };
+enum class ThermalisationScheme { INSTANT, DETAILED, BARNES };
 
 #endif
diff --git a/gammapkt.cc b/gammapkt.cc
@@ -969,35 +969,31 @@ void transport_gamma(Packet &pkt, double t2)
   }
 }
 
-void barnes_thermalization(Packet &pkt, bool local)
+void barnes_thermalisation(Packet &pkt)
 // 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);
+  // 0.1 is an average value to fit the analytic approximations from the paper.
+  // Alternative: Distinguish between low-E (kappa = 1) or high-E (kappa = 0.05)
+  // packets.
+  // constexpr double mean_gamma_opac = 0.1;
+
+  // determine average initial density via kinetic energy
+  const double E_kin = grid::get_ejecta_kinetic_energy();
+  const double v_ej = sqrt(E_kin * 2 / grid::mtot_input);
+  const double t_0 = globals::tmin;
+
+  // const double t_ineff = sqrt(rho_0 * R_0 * pow(t_0, 2) * mean_gamma_opac);
+  const double t_ineff = 1.4 * 86400. * sqrt(grid::mtot_input / (5.e-3 * 1.989 * 1.e33)) * ((0.2 * 29979200000) / v_ej);
   // get current time
   const double t = t_0 + pkt.prop_time;
-  const double tau = pow(tau_ineff / t, 2.);
+  const double tau = pow(t_ineff / t, 2.);
   const double f_gamma = 1. - exp(-tau);
+  assert_always(f_gamma >= 0.);
+  assert_always(f_gamma <= 1.);
 
   // either absorb packet or let it escape
   if (rng_uniform() < f_gamma) {
@@ -1014,10 +1010,8 @@ void barnes_thermalization(Packet &pkt, bool local)
 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 if constexpr (GAMMA_THERMALISATION_SCHEME == ThermalisationScheme::BARNES) {
+    barnes_thermalisation(pkt);
   } else {
     __builtin_unreachable();
   }

diff --git a/grid.cc b/grid.cc
@@ -626,7 +626,7 @@ static void set_elem_stable_abund_from_total(const int mgi, const int element, c
   modelgrid[mgi].composition[element].abundance = isofracsum + massfracstable;
 }
 
-static auto get_cellradialposmid(const int cellindex) -> double
+auto get_cellradialposmid(const int cellindex) -> double
 // get the radial distance from the origin to the centre of the cell at time tmin
 {
   if (GRID_TYPE == GRID_SPHERICAL1D) {

diff --git a/grid.h b/grid.h
@@ -82,6 +82,7 @@ void set_elements_uppermost_ion(int modelgridindex, int element, int newvalue);
 [[nodiscard]] auto get_cellcoordmax(int cellindex, int axis) -> double;
 [[nodiscard]] auto get_cellcoordmin(int cellindex, int axis) -> double;
 [[nodiscard]] auto get_cellcoordpointnum(int cellindex, int axis) -> int;
+[[nodiscard]] auto get_cellradialposmid(int cellindex) -> double;
 [[nodiscard]] auto get_coordcellindexincrement(int axis) -> int;
 [[nodiscard]] auto get_rho_tmin(int modelgridindex) -> float;
 [[nodiscard]] auto get_rho(int modelgridindex) -> float;
@@ -156,6 +157,19 @@ inline void change_cell(Packet &pkt, const int snext)
   }
 }
 
+inline auto get_ejecta_kinetic_energy() {
+  double E_kin = 0.;
+  for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) {
+    const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);
+    const int assoc_cells = grid::get_numassociatedcells(mgi);
+    double M_cell = grid::get_rho_tmin(mgi) * grid::get_modelcell_assocvolume_tmin(mgi);
+    const double radial_pos = grid::modelgrid[mgi].initial_radial_pos_sum / assoc_cells;
+    E_kin += 0.5 * M_cell * std::pow(radial_pos / globals::tmin, 2);
+  }
+
+  return E_kin;
+}
+
 }  // namespace grid
 
 #endif  // GRIDINIT_H
diff --git a/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/phixsdata_v2.txt b/tests/kilonova_2d_2dgrid_barnesthermalisation_inputfiles/phixsdata_v2.txt
@@ -0,0 +1,2 @@
+100
+ 3.0000000e-02