Update nonthermal.cc

nonthermal.cc

@@ -18,6 +18,7 @@
 #include <functional>
 #include <ios>
 #include <numeric>
+#include <ranges>
 #include <span>
 #include <sstream>
 #include <string>
@@ -1800,7 +1801,8 @@ void sfmatrix_add_excitation(std::vector<double> &sfmatrixuppertri, const int mo
   const int nlevels_all = get_nlevels(element, ion);
   const int nlevels = (nlevels_all > NTEXCITATION_MAXNLEVELS_LOWER) ? NTEXCITATION_MAXNLEVELS_LOWER : nlevels_all;
 
-  for (int lower = 0; lower < nlevels; lower++) {
+  const auto lowers = std::ranges::iota_view{0, nlevels};
+  std::for_each(lowers.begin(), lowers.end(), [&](const int lower) {
     const double statweight_lower = stat_weight(element, ion, lower);
     const double nnlevel = get_levelpop(nonemptymgi, element, ion, lower);
     const double epsilon_lower = epsilon(element, ion, lower);
@@ -1822,24 +1824,25 @@ void sfmatrix_add_excitation(std::vector<double> &sfmatrixuppertri, const int mo
       if (xsstartindex >= 0) {
         cblas_dscal(SFPTS - xsstartindex, DELTA_E, vec_xs_excitation_deltae.data() + xsstartindex, 1);
 
-        for (int i = 0; i < SFPTS; i++) {
+        const auto arri = std::ranges::iota_view{0, SFPTS};
+        std::for_each(EXEC_PAR arri.begin(), arri.end(), [&](const int i) {
           const int rowoffset = uppertriangular(i, 0);
           const double en = engrid(i);
           const int stopindex = get_energyindex_ev_lteq(en + epsilon_trans_ev);
 
           const int startindex = std::max(i, xsstartindex);
           for (int j = startindex; j < stopindex; j++) {
-            sfmatrixuppertri[rowoffset + j] += nnlevel * vec_xs_excitation_deltae[j];
+            atomicadd(sfmatrixuppertri[rowoffset + j], nnlevel * vec_xs_excitation_deltae[j]);
           }
 
           // do the last bit separately because we're not using the full delta_e interval
           const double delta_en_actual = (en + epsilon_trans_ev - engrid(stopindex));
-          sfmatrixuppertri[rowoffset + stopindex] +=
-              nnlevel * vec_xs_excitation_deltae[stopindex] * delta_en_actual / DELTA_E;
-        }
+          atomicadd(sfmatrixuppertri[rowoffset + stopindex],
+                    nnlevel * vec_xs_excitation_deltae[stopindex] * delta_en_actual / DELTA_E);
+        });
       }
     }
-  }
+  });
 }
 
 void sfmatrix_add_ionization(std::vector<double> &sfmatrixuppertri, const int Z, const int ionstage, const double nnion)