Attempt fix of thread local issues with classic

macroatom.cc

@@ -7,7 +7,9 @@
 #include <cstdio>
 #include <cstdlib>
 #include <iterator>
+#if defined(STDPAR_ON) || defined(_OPENMP_ON)
 #include <mutex>
+#endif
 #include <numeric>
 
 #include "artisoptions.h"
@@ -369,8 +371,10 @@ void do_macroatom(Packet &pkt, const MacroAtomState &pktmastate)
     auto &chlevel = globals::cellcache[cellcacheslotid].chelements[element].chions[ion].chlevels[level];
 
     {
+#if defined(STDPAR_ON) || defined(_OPENMP_ON)
       const auto lock =
           std::lock_guard<std::mutex>(globals::mutex_cellcachemacroatom[get_uniquelevelindex(element, ion, level)]);
+#endif
 
       assert_testmodeonly(globals::cellcache[cellcacheslotid].cellnumber == modelgridindex);
 

rpkt.cc

@@ -95,28 +95,37 @@ auto closest_transition(const double nu_cmf, const int next_trans) -> int
 // find the next transition lineindex redder than nu_cmf
 // return -1 if no transition can be reached
 {
-  if (next_trans > (globals::nlines - 1)) {
-    // packet is tagged as having no more line interactions
-    return -1;
-  }
   /// if nu_cmf is smaller than the lowest frequency in the linelist,
   /// no line interaction is possible: return negative value as a flag
   if (nu_cmf < globals::linelist[globals::nlines - 1].nu) {
     return -1;
   }
+  if (next_trans > (globals::nlines - 1)) {
+    // packet is tagged as having no more line interactions
+    return -1;
+  }
 
-  if (next_trans > 0) [[likely]] {
-    /// if left = pkt.next_trans > 0 we know the next line we should interact with, independent of the packets
+  if (next_trans > 0) {
+    /// if next_trans > 0 we know the next line we should interact with, independent of the packets
     /// current nu_cmf which might be smaller than globals::linelist[left].nu due to propagation errors
     return next_trans;
   }
+  if (nu_cmf >= globals::linelist[0].nu) {
+    /// if nu_cmf is larger than the highest frequency in the the linelist,
+    /// interaction with the first line occurs - no search
+    return 0;
+  }  /// otherwise go through the list until nu_cmf is located between two
+  /// entries in the line list and get the index of the closest line
+  /// to lower frequencies
+
   // will find the highest frequency (lowest index) line with nu_line <= nu_cmf
   // lower_bound matches the first element where the comparison function is false
-  const auto *matchline =
-      std::lower_bound(globals::linelist, globals::linelist + globals::nlines, nu_cmf,
-                       [](const auto &line, const double nu_cmf) -> bool { return line.nu > nu_cmf; });
-  const int matchindex = std::distance(globals::linelist, matchline);
-  if (matchindex >= globals::nlines) [[unlikely]] {
+  const int matchindex =
+      std::distance(globals::linelist,
+                    std::lower_bound(globals::linelist, globals::linelist + globals::nlines, nu_cmf,
+                                     [](const auto &line, const double nu_cmf) -> bool { return line.nu > nu_cmf; }));
+
+  if (matchindex >= globals::nlines) {
     return -1;
   }
 
@@ -747,17 +756,22 @@ static auto do_rpkt_step(Packet &pkt, const double t2) -> bool
   const int mgi = grid::get_cell_modelgridindex(cellindex);
   const int nonemptymgi = (mgi != grid::get_npts_model()) ? grid::get_modelcell_nonemptymgi(mgi) : -1;
 
-  static thread_local struct MacroAtomState pktmastate {};
-  static thread_local struct Phixslist phixslist {
-    .groundcont_gamma_contr = std::vector<double>(globals::nbfcontinua_ground, 0.),
-    .chi_bf_sum = std::vector<double>(globals::nbfcontinua, 0.),
-    .gamma_contr = std::vector<double>(globals::bfestimcount, 0.), .allcontend = 1, .allcontbegin = 0, .bfestimend = 1,
-    .bfestimbegin = 0,
-  };
+  MacroAtomState pktmastate{};
+
+  static Phixslist phixslist{};
+  phixslist.groundcont_gamma_contr.resize(globals::nbfcontinua_ground);
+  std::ranges::fill(phixslist.groundcont_gamma_contr, 0.);
+  phixslist.chi_bf_sum.resize(globals::nbfcontinua);
+  std::ranges::fill(phixslist.chi_bf_sum, 0.);
+  phixslist.gamma_contr.resize(globals::bfestimcount);
+  std::ranges::fill(phixslist.gamma_contr, 0.);
+  phixslist.allcontend = 1;
+  phixslist.allcontbegin = 0;
+  phixslist.bfestimend = 1;
+  phixslist.bfestimbegin = 0;
 
-  static thread_local struct Rpkt_continuum_absorptioncoeffs chi_rpkt_cont {
-    .nu = NAN, .total = NAN, .ffescat = NAN, .ffheat = NAN, .bf = NAN, .modelgridindex = -1, .timestep = -1
-  };
+  Rpkt_continuum_absorptioncoeffs chi_rpkt_cont{
+      .nu = NAN, .total = NAN, .ffescat = NAN, .ffheat = NAN, .bf = NAN, .modelgridindex = -1, .timestep = -1};
 
   // Assign optical depth to next physical event
   const double zrand = rng_uniform_pos();
@@ -1286,4 +1300,4 @@ void calculate_expansion_opacities(const int modelgridindex) {
     }
   }
   printout("took %ld seconds\n", std::time(nullptr) - sys_time_start_calc);
-}
+}