Update nltepop.cc

artis-mcrt · Nov 20, 2024 · cbedde4 · cbedde4
1 parent 6f3fd84
commit cbedde4
Showing 1 changed file with 14 additions and 14 deletions.
diff --git a/nltepop.cc b/nltepop.cc
@@ -15,6 +15,7 @@
 #include <cstdlib>
 #include <ctime>
 #include <ranges>
+#include <span>
 #include <tuple>
 #include <vector>
 
@@ -450,10 +451,8 @@ void nltepop_matrix_add_boundbound(const int nonemptymgi, const int element, con
 
     // de-excitation
     const int ndowntrans = get_ndowntrans(element, ion, level);
-    const auto *const leveldowntranslist = get_downtranslist(element, ion, level);
-    const auto ndowntransindices = std::ranges::iota_view{0, ndowntrans};
-    std::for_each(EXEC_PAR ndowntransindices.begin(), ndowntransindices.end(), [&](const auto i) {
-      const auto &downtransition = leveldowntranslist[i];
+    const auto leveldowntrans = std::span(get_downtranslist(element, ion, level), ndowntrans);
+    std::for_each(EXEC_PAR leveldowntrans.begin(), leveldowntrans.end(), [&](const auto &downtransition) {
       const double A_ul = downtransition.einstein_A;
       const int lower = downtransition.targetlevelindex;
 
@@ -529,7 +528,8 @@ void nltepop_matrix_add_ionisation(const int modelgridindex, const int element,
   const int nionisinglevels = get_nlevels_ionising(element, ion);
   const int maxrecombininglevel = get_maxrecombininglevel(element, ion + 1);
 
-  for (int level = 0; level < nionisinglevels; level++) {
+  const auto levels = std::ranges::iota_view{0, nionisinglevels};
+  std::for_each(EXEC_PAR levels.begin(), levels.end(), [&](const auto level) {
     const int lower_index = get_nlte_vector_index(element, ion, level);
 
     // thermal collisional ionization, photoionisation and recombination processes
@@ -548,10 +548,10 @@ void nltepop_matrix_add_ionisation(const int modelgridindex, const int element,
       const double C_ionisation =
           col_ionization_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans);
 
-      *gsl_matrix_ptr(rate_matrix_rad_bf, lower_index, lower_index) -= R_ionisation * s_renorm[level];
-      *gsl_matrix_ptr(rate_matrix_rad_bf, upper_index, lower_index) += R_ionisation * s_renorm[level];
-      *gsl_matrix_ptr(rate_matrix_coll_bf, lower_index, lower_index) -= C_ionisation * s_renorm[level];
-      *gsl_matrix_ptr(rate_matrix_coll_bf, upper_index, lower_index) += C_ionisation * s_renorm[level];
+      atomicadd(*gsl_matrix_ptr(rate_matrix_rad_bf, lower_index, lower_index), -R_ionisation * s_renorm[level]);
+      atomicadd(*gsl_matrix_ptr(rate_matrix_rad_bf, upper_index, lower_index), R_ionisation * s_renorm[level]);
+      atomicadd(*gsl_matrix_ptr(rate_matrix_coll_bf, lower_index, lower_index), -C_ionisation * s_renorm[level]);
+      atomicadd(*gsl_matrix_ptr(rate_matrix_coll_bf, upper_index, lower_index), C_ionisation * s_renorm[level]);
 
       if ((R_ionisation < 0) || (C_ionisation < 0)) {
         printout("  WARNING: Negative ionization rate from ionstage %d level %d phixstargetindex %d\n",
@@ -564,18 +564,18 @@ void nltepop_matrix_add_ionisation(const int modelgridindex, const int element,
         const double R_recomb = rad_recombination_ratecoeff(T_e, nne, element, ion + 1, upper, level, nonemptymgi);
         const double C_recomb = col_recombination_ratecoeff(T_e, nne, element, ion + 1, upper, level, epsilon_trans);
 
-        *gsl_matrix_ptr(rate_matrix_rad_bf, upper_index, upper_index) -= R_recomb * s_renorm[upper];
-        *gsl_matrix_ptr(rate_matrix_rad_bf, lower_index, upper_index) += R_recomb * s_renorm[upper];
-        *gsl_matrix_ptr(rate_matrix_coll_bf, upper_index, upper_index) -= C_recomb * s_renorm[upper];
-        *gsl_matrix_ptr(rate_matrix_coll_bf, lower_index, upper_index) += C_recomb * s_renorm[upper];
+        atomicadd(*gsl_matrix_ptr(rate_matrix_rad_bf, upper_index, upper_index), -R_recomb * s_renorm[upper]);
+        atomicadd(*gsl_matrix_ptr(rate_matrix_rad_bf, lower_index, upper_index), R_recomb * s_renorm[upper]);
+        atomicadd(*gsl_matrix_ptr(rate_matrix_coll_bf, upper_index, upper_index), -C_recomb * s_renorm[upper]);
+        atomicadd(*gsl_matrix_ptr(rate_matrix_coll_bf, lower_index, upper_index), C_recomb * s_renorm[upper]);
 
         if ((R_recomb < 0) || (C_recomb < 0)) {
           printout("  WARNING: Negative recombination rate to ionstage %d level %d phixstargetindex %d\n",
                    get_ionstage(element, ion), level, phixstargetindex);
         }
       }
     }
-  }
+  });
 }
 
 void nltepop_matrix_add_nt_ionisation(const int nonemptymgi, const int element, const int ion,