Add microclumping

artisoptions_christinenonthermal.h

@@ -159,5 +159,11 @@ constexpr auto GAMMA_THERMALISATION_SCHEME = GammaThermalisationScheme::FREQUENC
 
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+constexpr bool USE_MICROCLUMPING = false;
+
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE = false;
+
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_classic.h

@@ -154,5 +154,11 @@ constexpr auto GAMMA_THERMALISATION_SCHEME = GammaThermalisationScheme::FREQUENC
 
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+constexpr bool USE_MICROCLUMPING = false;
+
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE = false;
+
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_doc.md

@@ -261,4 +261,17 @@ constexpr bool USE_XCOM_GAMMAPHOTOION;
 // Override frequency-dependent gamma-ray opacity with a grey opacity [cm^2/g]
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+// Use microclumping, which enhances collisional (de)excitation, collisional ionisation, collisional recombination,
+// radiative recombination, collisional capture, stimulated recombination, free-free heating, free-free cooling by the
+// reciprocal of the volume filling factor (the "clumping factor")
+constexpr bool USE_MICROCLUMPING;
+
+// Read the volume filling factors from a plain text file with (num timesteps) rows * (num cells) cols.
+// Entries should be in (0,1] and formatted using %.6e.
+// Has no effect if USE_MICROCLUMPING is set to false.
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE;
+
+// Calculate volume filling factors based on time and radial velocity
+// Will be passed globals::timesteps[nts].mid and grid::get_modelcell_mean_radial_vel(mgi, globals::tmin)
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
 ```

artisoptions_kilonova_lte.h

@@ -154,5 +154,11 @@ constexpr auto GAMMA_THERMALISATION_SCHEME = GammaThermalisationScheme::FREQUENC
 
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+constexpr bool USE_MICROCLUMPING = false;
+
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE = false;
+
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_nltenebular.h

@@ -161,5 +161,11 @@ constexpr auto GAMMA_THERMALISATION_SCHEME = GammaThermalisationScheme::FREQUENC
 
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+constexpr bool USE_MICROCLUMPING = false;
+
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE = false;
+
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_nltephotospheric_dynamic_ion_range.h

@@ -162,5 +162,11 @@ constexpr auto GAMMA_THERMALISATION_SCHEME = GammaThermalisationScheme::FREQUENC
 
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+constexpr bool USE_MICROCLUMPING = false;
+
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE = false;
+
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

artisoptions_nltewithoutnonthermal.h

@@ -157,5 +157,11 @@ constexpr auto GAMMA_THERMALISATION_SCHEME = GammaThermalisationScheme::FREQUENC
 
 constexpr std::optional<double> GAMMA_USE_KAPPA_GREY;
 
+constexpr bool USE_MICROCLUMPING = false;
+
+constexpr bool READ_VOLUME_FILLING_FACTORS_FROM_FILE = false;
+
+constexpr float volume_filling_factor([[maybe_unused]] double tmid, [[maybe_unused]] double rad_vel) { return 1.; }
+
 // NOLINTEND(modernize*,misc-unused-parameters)
 #endif  // ARTISOPTIONS_H

grid.cc

@@ -453,7 +453,11 @@ void allocate_nonemptymodelcells() {
   kappagrey_allcells = MPI_shared_array<float>(nonempty_npts_model, 0.);
   grey_depth_allcells = MPI_shared_array<float>(nonempty_npts_model, 0.);
   thick_allcells = MPI_shared_array<int>(nonempty_npts_model, 0);
-  const auto modelgrid_mem_usage = nonempty_npts_model * ((sizeof(float) * 9) + sizeof(double) + sizeof(int));
+  if constexpr (USE_MICROCLUMPING) {
+    clumpfactor_allcells = MPI_shared_array<float>(nonempty_npts_model, -1.);
+  }
+  const auto modelgrid_mem_usage =
+      nonempty_npts_model * ((sizeof(float) * (USE_MICROCLUMPING ? 10 : 9)) + sizeof(double) + sizeof(int));
   printlnlog(
       "[info] mem_usage: the modelgrid properties (temperatures and electron densities) occupies {:.3f} MB (node "
       "shared memory)",
@@ -1506,6 +1510,18 @@ auto get_rho_tmin(const int modelgridindex) -> float { return modelgrid_input[mo
   return nne_allcells[nonemptymgi];
 }
 
+[[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_clumpfactor(const int nonemptymgi) -> float {
+  assert_testmodeonly(nonemptymgi >= 0);
+  assert_testmodeonly(nonemptymgi < get_nonempty_npts_model());
+
+  if constexpr (USE_MICROCLUMPING) {
+    assert_testmodeonly(std::isfinite(clumpfactor_allcells[nonemptymgi]) && clumpfactor_allcells[nonemptymgi] >= 1.F);
+    return clumpfactor_allcells[nonemptymgi];
+  }
+
+  return 1.;
+}
+
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_nnetot(const int nonemptymgi) -> float {
   assert_testmodeonly(nonemptymgi >= 0);
   assert_testmodeonly(nonemptymgi < get_nonempty_npts_model());
@@ -1522,6 +1538,16 @@ auto get_rho_tmin(const int modelgridindex) -> float { return modelgrid_input[mo
   return modelgrid_input[modelgridindex].initial_radial_pos_sum * tratmid / assoc_cells;
 }
 
+[[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_modelcell_mean_radial_vel(const int modelgridindex, const double tmin)
+    -> double {
+  // TODO: both versions give the same thing. uncommented one should give ever so slightly better performance
+
+  const int assoc_cells = grid::get_numpropcells(modelgridindex);
+  return modelgrid_input[modelgridindex].initial_radial_pos_sum / tmin / assoc_cells;
+  // ---------------
+  // return get_modelcell_mean_radial_pos(modelgridindex, tratmid) / tmid;
+}
+
 // mass fraction of an element (all isotopes combined)
 [[gnu::pure]] [[nodiscard]] auto get_elem_abundance(const std::ptrdiff_t nonemptymgi, const int element) -> float {
   const auto massfrac = elem_massfracs_allcells[(nonemptymgi * get_nelements()) + element];
@@ -1579,6 +1605,12 @@ void set_nne(const int nonemptymgi, const float nne) {
   nne_allcells[nonemptymgi] = nne;
 }
 
+void set_clumpfactor(const int nonemptymgi, const float clumpfactor) {
+  assert_testmodeonly(USE_MICROCLUMPING);
+  assert_always(std::isfinite(clumpfactor) && clumpfactor >= 1);
+  clumpfactor_allcells[nonemptymgi] = clumpfactor;
+}
+
 // Calculate and set the total density of electrons (free and bound) in grid cell. These are targets for Compton
 // scattering of gamma rays
 void set_nnetot(const int nonemptymgi) {
@@ -1678,6 +1710,18 @@ DEVICE_FUNC auto get_modelgridtype() -> GridType {
   return nonemptymgi;
 }
 
+// Check if cell with index `mgi` is a non-empty cell
+// If the cell is non-empty will return `true` and `nonemptymgi` will be set accordingly, otherwise returns `false` and
+// `nonemptymgi` is set to -1
+[[nodiscard]] DEVICE_FUNC auto check_mgi_is_nonempty(const int mgi, int& nonemptymgi) -> bool {
+  assert_testmodeonly(get_nonempty_npts_model() > 0);
+  assert_testmodeonly(mgi >= 0);
+  assert_testmodeonly(mgi < get_npts_model());
+
+  nonemptymgi = nonemptymgi_of_mgi[mgi];
+  return nonemptymgi >= 0;
+}
+
 // get the index in the list of non-empty cells for a given model grid cell
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_mgi_of_nonemptymgi(const ptrdiff_t nonemptymgi) -> int {
   assert_testmodeonly(get_nonempty_npts_model() > 0);

grid.h

@@ -25,6 +25,7 @@ inline MPI_shared_array<float> kappagrey_allcells;
 inline MPI_shared_array<float> grey_depth_allcells;  // Grey optical depth to surface of the modelgridcell
 inline MPI_shared_array<int>
     thick_allcells;  // whether the cell is optically thick (1) or not (0), or (2) thick for vpkts only
+inline MPI_shared_array<float> clumpfactor_allcells;
 
 inline ptrdiff_t ngrid{0};
 
@@ -52,6 +53,7 @@ void set_elements_uppermost_ion(int nonemptymgi, int element, int uppermost_ion)
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_ffegrp(int modelgridindex) -> float;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_modelcell_mean_radial_pos(int modelgridindex, double tratmid)
     -> double;
+[[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_modelcell_mean_radial_vel(int modelgridindex, double tmin) -> double;
 void set_elem_abundance(std::ptrdiff_t nonemptymgi, int element, float newabundance);
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_elem_numberdens(std::ptrdiff_t nonemptymgi, int element) -> double;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_initenergyq(int modelgridindex) -> double;
@@ -78,6 +80,7 @@ void set_element_meanweight(std::ptrdiff_t nonemptymgi, int element, float meanw
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_numpropcells(int modelgridindex) -> int;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_nonemptymgi_of_mgi(int mgi) -> int;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_mgi_of_nonemptymgi(std::ptrdiff_t nonemptymgi) -> int;
+[[nodiscard]] DEVICE_FUNC auto check_mgi_is_nonempty(int mgi, int& nonemptymgi) -> bool;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_modelgridtype() -> GridType;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_npts_model() -> int;
 [[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_nonempty_npts_model() -> int;
@@ -95,6 +98,8 @@ void write_grid_restart_data(int timestep);
 [[nodiscard]] auto get_propcell_random_position_tmin(int cellindex) -> Vec3d;
 [[nodiscard]] DEVICE_FUNC auto boundary_distance(const Vec3d& dir, const Vec3d& pos, double tstart, int cellindex)
     -> std::tuple<double, int>;
+void set_clumpfactor(int nonemptymgi, float clumpfactor);
+[[gnu::pure]] [[nodiscard]] DEVICE_FUNC auto get_clumpfactor(int nonemptymgi) -> float;
 
 [[nodiscard]] auto calculate_cell_kappagrey(int nonemptymgi) -> float;
 

kpkt.cc

@@ -52,6 +52,7 @@ auto calculate_cooling_rates_ion(const int nonemptymgi, const int element, const
                                  double* const C_exc, double* const C_ionisation) -> double {
   const auto nne = grid::get_nne(nonemptymgi);
   const auto T_e = grid::get_Te(nonemptymgi);
+  const float clumpfactor = grid::get_clumpfactor(nonemptymgi);
 
   double C_ion = 0.;
   [[maybe_unused]] int i = 0;  // NOLINT(misc-const-correctness)
@@ -63,7 +64,7 @@ auto calculate_cooling_rates_ion(const int nonemptymgi, const int element, const
   const int ioncharge = get_ionstage(element, ion) - 1;
   if (ioncharge > 0) {
     const double C_ff_ion = 1.426e-27 * sqrt(T_e) * pow(ioncharge, 2) * nncurrention * nne;
-    C_ion += C_ff_ion;
+    C_ion += C_ff_ion * clumpfactor;
 
     if constexpr (update_cellcache_contribs) {
       ion_contribs[i] = C_ion;
@@ -93,9 +94,10 @@ auto calculate_cooling_rates_ion(const int nonemptymgi, const int element, const
       const int upper = globals::alltrans.targetlevelindex[alltransindex];
       const double epsilon_trans = epsilon(ionuniquelevelindexstart + upper) - epsilon_current;
       const auto upper_statweight = stat_weight(ionuniquelevelindexstart + upper);
-      const double C = nnlevel *
-                       col_excitation_ratecoeff(T_e, nne, upper_statweight, alltransindex, epsilon_trans, statweight) *
-                       epsilon_trans;
+      const double C =
+          nnlevel *
+          col_excitation_ratecoeff(T_e, nne, upper_statweight, alltransindex, epsilon_trans, statweight, clumpfactor) *
+          epsilon_trans;
       C_ion += C;
       if constexpr (!update_cellcache_contribs) {
         *C_exc += C;
@@ -126,9 +128,10 @@ auto calculate_cooling_rates_ion(const int nonemptymgi, const int element, const
         const int upper = get_phixsupperlevel(uniquelevelindex, phixstargetindex);
         const double epsilon_upper = epsilon(element, ion + 1, upper);
         const double epsilon_trans = epsilon_upper - epsilon_current;
-        const double C = nnlevel *
-                         col_ionisation_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans) *
-                         epsilon_trans;
+        const double C =
+            nnlevel *
+            col_ionisation_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans, clumpfactor) *
+            epsilon_trans;
 
         C_ion += C;
         if constexpr (update_cellcache_contribs) {
@@ -159,7 +162,7 @@ auto calculate_cooling_rates_ion(const int nonemptymgi, const int element, const
           return nnupperion;
         }();
         const double C = get_bfcoolingcoeff(element, ion, level, phixstargetindex, T_e) * pop * nne;
-        C_ion += C;
+        C_ion += C * clumpfactor;
 
         if constexpr (update_cellcache_contribs) {
           ion_contribs[i] = C_ion;
@@ -479,6 +482,7 @@ DEVICE_FUNC void do_kpkt(Packet& pkt, const double t2, const int nts) {
   } else if (rndcoolingtype == CoolingType::COLLEXC) {
     // the k-packet activates a macro-atom due to collisional excitation
     const float nne = grid::get_nne(nonemptymgi);
+    const float clumpfactor = grid::get_clumpfactor(nonemptymgi);
 
     // if the previous entry belongs to the same ion, then pick up the cumulative sum from
     // the previous entry, otherwise start from zero
@@ -500,9 +504,10 @@ DEVICE_FUNC void do_kpkt(Packet& pkt, const double t2, const int nts) {
       const auto upperuniquelevelindex = ionuniquelevelindexstart + tmpupper;
       const double epsilon_trans = epsilon(upperuniquelevelindex) - epsilon_current;
       const auto upper_statweight = stat_weight(upperuniquelevelindex);
-      const double C = nnlevel *
-                       col_excitation_ratecoeff(T_e, nne, upper_statweight, alltransindex, epsilon_trans, statweight) *
-                       epsilon_trans;
+      const double C =
+          nnlevel *
+          col_excitation_ratecoeff(T_e, nne, upper_statweight, alltransindex, epsilon_trans, statweight, clumpfactor) *
+          epsilon_trans;
       contrib += C;
       if (contrib > rndcool_ion_process) {
         upper = tmpupper;