Add option to limit number of BF estimators

artisoptions_christinenonthermal.h

@@ -81,6 +81,18 @@ constexpr bool DETAILED_BF_ESTIMATORS_ON = true;
 
 constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13;
 
+constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) {
+  bool custom_bf_estimator = false;
+  if (custom_bf_estimator == true) {
+    if (element_z == 26 && ionstage == 2) {
+      return (level <= 197);
+    }
+    return (level <= 80);
+  } else {
+    return true;  // Always treat (element, ion, level) with bf estimator
+  }
+}
+
 constexpr bool USE_LUT_PHOTOION = false;
 
 constexpr bool USE_LUT_BFHEATING = false;

artisoptions_classic.h

@@ -75,6 +75,18 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false;
 
 constexpr bool DETAILED_BF_ESTIMATORS_ON = false;
 
+constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) {
+  bool custom_bf_estimator = false;
+  if (custom_bf_estimator == true) {
+    if (element_z == 26 && ionstage == 2) {
+      return (level <= 197);
+    }
+    return (level <= 80);
+  } else {
+    return true;  // Always treat (element, ion, level) with bf estimator
+  }
+}
+
 constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13;
 
 constexpr bool USE_LUT_PHOTOION = true;

artisoptions_kilonova_lte.h

@@ -74,6 +74,18 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false;
 
 constexpr bool DETAILED_BF_ESTIMATORS_ON = false;
 
+constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) {
+  bool custom_bf_estimator = false;
+  if (custom_bf_estimator == true) {
+    if (element_z == 26 && ionstage == 2) {
+      return (level <= 197);
+    }
+    return (level <= 80);
+  } else {
+    return true;
+  }
+}
+
 constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13;
 
 constexpr bool USE_LUT_PHOTOION = true;

artisoptions_nltenebular.h

@@ -6,7 +6,7 @@
 
 #include "constants.h"
 
-constexpr int MPKTS = 1000000;
+constexpr int MPKTS = 750000;
 
 constexpr int GRID_TYPE = GRID_CARTESIAN3D;
 constexpr int CUBOID_NCOORDGRID_X = 50;
@@ -82,6 +82,18 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false;
 
 constexpr bool DETAILED_BF_ESTIMATORS_ON = true;
 
+constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) {
+  bool custom_bf_estimator = false;
+  if (custom_bf_estimator == true) {
+    if (element_z == 26 && ionstage == 2) {
+      return (level <= 197);
+    }
+    return (level <= 80);
+  } else {
+    return true;
+  }
+}
+
 constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13;
 
 constexpr bool USE_LUT_PHOTOION = false;
@@ -146,4 +158,4 @@ constexpr bool KEEP_ALL_RESTART_FILES = false;
 constexpr bool BFCOOLING_USELEVELPOPNOTIONPOP = false;
 
 // NOLINTEND(modernize*,misc-unused-parameters)
-#endif  // ARTISOPTIONS_H
+#endif  // ARTISOPTIONS_H

artisoptions_nltewithoutnonthermal.h

@@ -78,6 +78,18 @@ constexpr bool DETAILED_LINE_ESTIMATORS_ON = false;
 
 constexpr bool DETAILED_BF_ESTIMATORS_ON = true;
 
+constexpr bool LEVEL_HAS_BFEST(int element_z, int ionstage, int level) {
+  bool custom_bf_estimator = false;
+  if (custom_bf_estimator == true) {
+    if (element_z == 26 && ionstage == 2) {
+      return (level <= 197);
+    }
+    return (level <= 80);
+  } else {
+    return true;  // Always treat (element, ion, level) with bf estimator
+  }
+}
+
 constexpr int DETAILED_BF_ESTIMATORS_USEFROMTIMESTEP = 13;
 
 constexpr bool USE_LUT_PHOTOION = false;

globals.cc

@@ -12,6 +12,8 @@ std::unique_ptr<struct time[]> timesteps = nullptr;
 
 double *rpkt_emiss = nullptr;  /// Volume estimator for the rpkt emissivity
 
+int BFCounter;
+
 // for USE_LUT_PHOTOION = true
 double *corrphotoionrenorm = nullptr;
 double *gammaestimator = nullptr;

globals.h

@@ -50,6 +50,7 @@ struct fullphixslist {
   float *photoion_xs;
   double probability;
   int index_in_groundphixslist;
+  bool has_bf_estimator;
 };
 
 struct groundphixslist {
@@ -209,6 +210,8 @@ extern std::unique_ptr<struct time[]> timesteps;
 
 extern double *rpkt_emiss;
 
+extern int BFCounter;
+
 // for USE_LUT_PHOTOION = true
 extern double *corrphotoionrenorm;
 extern double *gammaestimator;

input.cc

@@ -1356,6 +1356,13 @@ static void write_bflist_file(int includedphotoiontransitions) {
   }
 }
 
+auto is_bfest(const int element, const int ion, const int level) -> bool
+// Returns true if (element,ion,level) is to be treated with a bf estimator.
+{
+  return LEVEL_HAS_BFEST(get_atomicnumber(element), get_ionstage(element, ion),
+                         level);  // defined in artisoptions.h
+}
+
 static void setup_phixs_list() {
   // set up the photoionisation transition lists
   // and temporary gamma/kappa lists for each thread
@@ -1456,6 +1463,8 @@ static void setup_phixs_list() {
            globals::nbfcontinua * (sizeof(fullphixslist)) / 1024. / 1024.);
   size_t nbftables = 0;
   int allcontindex = 0;
+  int bound_free_levels_counter = 0;
+  int nonbound_free_levels_counter = 0;
   for (int element = 0; element < get_nelements(); element++) {
     const int nions = get_nions(element);
     for (int ion = 0; ion < nions - 1; ion++) {
@@ -1482,6 +1491,15 @@ static void setup_phixs_list() {
           nonconstallcont[allcontindex].probability = get_phixsprobability(element, ion, level, phixstargetindex);
           nonconstallcont[allcontindex].upperlevel = get_phixsupperlevel(element, ion, level, phixstargetindex);
 
+          if (is_bfest(element, ion, level)) {
+            nonconstallcont[allcontindex].has_bf_estimator = true;
+            bound_free_levels_counter += 1;
+
+          } else {
+            nonconstallcont[allcontindex].has_bf_estimator = false;
+            nonbound_free_levels_counter += 1;
+          }
+
           if constexpr (USE_LUT_PHOTOION || USE_LUT_BFHEATING) {
             int index_in_groundlevelcontestimator = 0;
             nonconstallcont[allcontindex].index_in_groundphixslist =
@@ -1495,6 +1513,12 @@ static void setup_phixs_list() {
       }
     }
   }
+  int total_bound_free = nonbound_free_levels_counter + bound_free_levels_counter;
+  printout("[info] Bf estimator flag passed through:    \t%d\n", total_bound_free);
+  printout("[info] Bf estimator flag activated for:     \t%d\n", bound_free_levels_counter);
+  printout("[info] Bf estimator flag not activated for: \t%d\n", nonbound_free_levels_counter);
+  globals::BFCounter = bound_free_levels_counter;
+  printout("[info] Set the BFCounter to:\t%d\n", bound_free_levels_counter);
 
   assert_always(allcontindex == globals::nbfcontinua);
   assert_always(globals::nbfcontinua >= 0);  // was initialised as -1 before startup
@@ -2189,4 +2213,4 @@ void write_timestep_file() {
             globals::timesteps[n].width / DAY);
   }
   fclose(timestepfile);
-}
+}

input.h

@@ -11,6 +11,7 @@ void update_parameterfile(int nts);
 void time_init();
 void write_timestep_file();
 auto get_noncommentline(std::fstream &input, std::string &line) -> bool;
+auto is_bfest(int element, int ion, int level) -> bool;
 
 static inline auto lineiscommentonly(const std::string &line) -> bool
 // return true for whitepace-only lines, and lines that are exclusively whitepace up to a '#' character

radfield.cc

@@ -368,10 +368,10 @@ void init(int my_rank, int ndo_nonempty)
     printout("[info] mem_usage: detailed bf estimators for non-empty cells occupy %.3f MB (node shared memory)\n",
              nonempty_npts_model * globals::nbfcontinua * sizeof(float) / 1024. / 1024.);
 
-    bfrate_raw = static_cast<double *>(malloc(nonempty_npts_model * globals::nbfcontinua * sizeof(double)));
+    bfrate_raw = static_cast<double *>(malloc(nonempty_npts_model * globals::BFCounter * sizeof(double)));
 
     printout("[info] mem_usage: detailed bf estimator acculumators for non-empty cells occupy %.3f MB\n",
-             nonempty_npts_model * globals::nbfcontinua * sizeof(double) / 1024. / 1024.);
+             nonempty_npts_model * globals::BFCounter * sizeof(double) / 1024. / 1024.);
   }
 
   for (int modelgridindex = 0; modelgridindex < grid::get_npts_model(); modelgridindex++) {
@@ -655,8 +655,8 @@ void zero_estimators(int modelgridindex)
   if constexpr (DETAILED_BF_ESTIMATORS_ON) {
     assert_always(bfrate_raw != nullptr);
     if (grid::get_numassociatedcells(modelgridindex) > 0) {
-      for (int i = 0; i < globals::nbfcontinua; i++) {
-        bfrate_raw[nonemptymgi * globals::nbfcontinua + i] = 0.;
+      for (int i = 0; i < globals::BFCounter; i++) {
+        bfrate_raw[nonemptymgi * globals::BFCounter + i] = 0.;
       }
     }
   }
@@ -703,18 +703,42 @@ static void update_bfestimators(const int nonemptymgi, const double distance_e_c
   const int tid = get_thread_num();
   const double distance_e_cmf_over_nu =
       distance_e_cmf / nu_cmf * dopplerfactor;  // TODO: Luke: why did I put a doppler factor here?
+  int detailed_counter = 0;
   for (int allcontindex = 0; allcontindex < nbfcontinua; allcontindex++) {
     const double nu_edge = globals::allcont_nu_edge[allcontindex];
     const double nu_max_phixs = nu_edge * last_phixs_nuovernuedge;  // nu of the uppermost point in the phixs table
 
     if (nu_cmf >= nu_edge && nu_cmf <= nu_max_phixs) {
-      safeadd(bfrate_raw[nonemptymgi * nbfcontinua + allcontindex],
-              globals::phixslist[tid].gamma_contr[allcontindex] * distance_e_cmf_over_nu);
+      if (globals::allcont[allcontindex].has_bf_estimator == true) {
+        if (detailed_counter == globals::BFCounter) {
+          printout("detailed_counter %d\n", detailed_counter);
+          printout("globals::BFCounter \n", globals::BFCounter);
+          printout("allcontindex %d\n", allcontindex);
+        }
+        assert_always(detailed_counter < globals::BFCounter);
+
+        if (nonemptymgi * globals::BFCounter + detailed_counter >=
+            grid::get_nonempty_npts_model() * globals::BFCounter) {
+          printout("detailed_counter %d\n", detailed_counter);
+          printout("globals::BFCounter %d\n", globals::BFCounter);
+          printout("grid::get_nonempty_npts_model() %d\n", grid::get_nonempty_npts_model());
+          printout("nonemptymgi %d\n", nonemptymgi);
+        }
+        assert_always(nonemptymgi * globals::BFCounter + detailed_counter <
+                      grid::get_nonempty_npts_model() * globals::BFCounter)
+
+            safeadd(bfrate_raw[nonemptymgi * globals::BFCounter + detailed_counter],
+                    globals::phixslist[tid].gamma_contr[allcontindex] * distance_e_cmf_over_nu);
+      }
 
     } else if (nu_cmf < nu_edge) {
       // list is sorted by nu_edge, so all remaining will have nu_cmf < nu_edge
       break;
     }
+
+    if (globals::allcont[allcontindex].has_bf_estimator == true) {
+      detailed_counter += 1;
+    }
   }
 }
 
@@ -1161,9 +1185,16 @@ void normalise_bf_estimators(const int modelgridindex, const double estimator_no
     printout("normalise_bf_estimators for cell %d with factor %g\n", modelgridindex, estimator_normfactor_over_H);
     const int nonemptymgi = grid::get_modelcell_nonemptymgi(modelgridindex);
     assert_always(nonemptymgi >= 0);
+    int detailed_counter = 0;
     for (int i = 0; i < globals::nbfcontinua; i++) {
       const int mgibfindex = nonemptymgi * globals::nbfcontinua + i;
-      prev_bfrate_normed[mgibfindex] = bfrate_raw[mgibfindex] * estimator_normfactor_over_H;
+      const int detailed_mgibfindex = nonemptymgi * globals::BFCounter + detailed_counter;
+      if (globals::allcont[i].has_bf_estimator == true) {
+        prev_bfrate_normed[mgibfindex] = bfrate_raw[detailed_mgibfindex] * estimator_normfactor_over_H;
+        detailed_counter += 1;
+      } else if (globals::allcont[i].has_bf_estimator == false) {
+        prev_bfrate_normed[mgibfindex] = 0;
+      }
     }
   }
 }
@@ -1257,7 +1288,7 @@ void reduce_estimators()
   MPI_Allreduce(MPI_IN_PLACE, nuJ.data(), nonempty_npts_model, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 
   if constexpr (DETAILED_BF_ESTIMATORS_ON) {
-    MPI_Allreduce(MPI_IN_PLACE, bfrate_raw, nonempty_npts_model * globals::nbfcontinua, MPI_DOUBLE, MPI_SUM,
+    MPI_Allreduce(MPI_IN_PLACE, bfrate_raw, grid::get_nonempty_npts_model() * globals::BFCounter, MPI_DOUBLE, MPI_SUM,
                   MPI_COMM_WORLD);
   }