diff --git a/nltepop.cc b/nltepop.cc index 0bf69ab6e..168518b89 100644 --- a/nltepop.cc +++ b/nltepop.cc @@ -600,9 +600,10 @@ void nltepop_matrix_add_nt_ionisation(const int nonemptymgi, const int element, for (int level = 0; level < nlevels; level++) { const int lower_index = get_nlte_vector_index(element, ion, level); - *gsl_matrix_ptr(rate_matrix_ntcoll_bf, lower_index, lower_index) -= Y_nt_thisupperion * s_renorm[level]; - *gsl_matrix_ptr(rate_matrix_ntcoll_bf, upper_groundstate_index, lower_index) += - Y_nt_thisupperion * s_renorm[level]; + atomicadd(*gsl_matrix_ptr(rate_matrix_ntcoll_bf, lower_index, lower_index), + -Y_nt_thisupperion * s_renorm[level]); + atomicadd(*gsl_matrix_ptr(rate_matrix_ntcoll_bf, upper_groundstate_index, lower_index), + Y_nt_thisupperion * s_renorm[level]); } } } @@ -899,7 +900,8 @@ void solve_nlte_pops_element(const int element, const int nonemptymgi, const int } // printout(" Adding rates for ion stages:"); - for (int ion = 0; ion < nions; ion++) { + const auto ions = std::ranges::iota_view{0, nions}; + std::for_each(EXEC_PAR ions.begin(), ions.end(), [&](const auto ion) { // const int ionstage = get_ionstage(element, ion); // printout(" %d", ionstage); @@ -916,14 +918,14 @@ void solve_nlte_pops_element(const int element, const int nonemptymgi, const int nltepop_matrix_add_boundbound(nonemptymgi, element, ion, t_mid, s_renorm, &rate_matrix_rad_bb, &rate_matrix_coll_bb, &rate_matrix_ntcoll_bb); - if (ion < nions - 1) { + if (ion < (nions - 1)) { // this is the slowest component nltepop_matrix_add_ionisation(modelgridindex, element, ion, s_renorm, &rate_matrix_rad_bf, &rate_matrix_coll_bf); if (NT_ON) { nltepop_matrix_add_nt_ionisation(nonemptymgi, element, ion, s_renorm, &rate_matrix_ntcoll_bf); } } - } + }); // printout("\n"); if (individual_process_matrices) {