Pass heatingcoolingrates as reference

nonthermal.cc

@@ -2104,40 +2104,39 @@ void init() {
 // set total non-thermal deposition rate from individual gamma/positron/electron/alpha rates. This should be called
 // after packet propagation is finished for this timestep and normalise_deposition_estimators() has been done
 void calculate_deposition_rate_density(const int nonemptymgi, const int timestep,
-                                       HeatingCoolingRates *heatingcoolingrates) {
-  heatingcoolingrates->dep_gamma = globals::dep_estimator_gamma[nonemptymgi];
+                                       HeatingCoolingRates &heatingcoolingrates) {
+  heatingcoolingrates.dep_gamma = globals::dep_estimator_gamma[nonemptymgi];
 
   const double tmid = globals::timesteps[timestep].mid;
   const double rho = grid::get_rho(nonemptymgi);
 
   // if INSTANT_PARTICLE_DEPOSITION, use the analytic rate at t_mid since it will have no Monte Carlo noise (although
   // strictly, it should be an integral from the timestep start to the end)
   // with time-dependent deposition, we don't have an analytic rate, so we use the Monte Carlo rate
-  assert_always(heatingcoolingrates != nullptr);
 
-  heatingcoolingrates->eps_gamma_ana = rho * decay::get_gamma_emission_rate(nonemptymgi, tmid);
+  heatingcoolingrates.eps_gamma_ana = rho * decay::get_gamma_emission_rate(nonemptymgi, tmid);
 
-  heatingcoolingrates->eps_positron_ana =
+  heatingcoolingrates.eps_positron_ana =
       rho * decay::get_particle_injection_rate(nonemptymgi, tmid, decay::DECAYTYPE_BETAPLUS);
 
-  heatingcoolingrates->eps_electron_ana =
+  heatingcoolingrates.eps_electron_ana =
       (rho * decay::get_particle_injection_rate(nonemptymgi, tmid, decay::DECAYTYPE_BETAMINUS));
 
-  heatingcoolingrates->eps_alpha_ana =
+  heatingcoolingrates.eps_alpha_ana =
       rho * decay::get_particle_injection_rate(nonemptymgi, tmid, decay::DECAYTYPE_ALPHA);
 
   if (PARTICLE_THERMALISATION_SCHEME == ThermalisationScheme::INSTANT) {
-    heatingcoolingrates->dep_positron = heatingcoolingrates->eps_positron_ana;
-    heatingcoolingrates->dep_electron = heatingcoolingrates->eps_electron_ana;
-    heatingcoolingrates->dep_alpha = heatingcoolingrates->eps_alpha_ana;
+    heatingcoolingrates.dep_positron = heatingcoolingrates.eps_positron_ana;
+    heatingcoolingrates.dep_electron = heatingcoolingrates.eps_electron_ana;
+    heatingcoolingrates.dep_alpha = heatingcoolingrates.eps_alpha_ana;
   } else {
-    heatingcoolingrates->dep_positron = globals::dep_estimator_positron[nonemptymgi];
-    heatingcoolingrates->dep_electron = globals::dep_estimator_electron[nonemptymgi];
-    heatingcoolingrates->dep_alpha = globals::dep_estimator_alpha[nonemptymgi];
+    heatingcoolingrates.dep_positron = globals::dep_estimator_positron[nonemptymgi];
+    heatingcoolingrates.dep_electron = globals::dep_estimator_electron[nonemptymgi];
+    heatingcoolingrates.dep_alpha = globals::dep_estimator_alpha[nonemptymgi];
   }
 
   deposition_rate_density_all_cells[nonemptymgi] =
-      (heatingcoolingrates->dep_gamma + heatingcoolingrates->dep_positron + heatingcoolingrates->dep_electron);
+      (heatingcoolingrates.dep_gamma + heatingcoolingrates.dep_positron + heatingcoolingrates.dep_electron);
 }
 
 // get non-thermal deposition rate density in erg / s / cm^3 previously stored by calculate_deposition_rate_density()

nonthermal.h

@@ -15,7 +15,7 @@ void solve_spencerfano(int nonemptymgi, int timestep, int iteration);
                                                       bool energyweighted) -> double;
 [[nodiscard]] auto nt_ionisation_maxupperion(int element, int lowerion) -> int;
 [[nodiscard]] auto nt_random_upperion(int nonemptymgi, int element, int lowerion, bool energyweighted) -> int;
-void calculate_deposition_rate_density(int nonemptymgi, int timestep, HeatingCoolingRates *heatingcoolingrates);
+void calculate_deposition_rate_density(int nonemptymgi, int timestep, HeatingCoolingRates &heatingcoolingrates);
 [[nodiscard]] auto get_deposition_rate_density(int nonemptymgi) -> double;
 [[nodiscard]] auto get_nt_frac_heating(int modelgridindex) -> float;
 #pragma omp declare simd

thermalbalance.cc

@@ -132,7 +132,7 @@ auto get_heating_ion_coll_deexc(const int nonemptymgi, const int element, const
 // Calculate the heating rates for a given cell. Results are returned via the elements of the heatingrates data
 // structure.
 void calculate_heating_rates(const int nonemptymgi, const double T_e, const double nne,
-                             HeatingCoolingRates *heatingcoolingrates, const std::vector<double> &bfheatingcoeffs) {
+                             HeatingCoolingRates &heatingcoolingrates, const std::vector<double> &bfheatingcoeffs) {
   double C_deexc = 0.;
 
   // double C_recomb = 0.;
@@ -172,14 +172,14 @@ void calculate_heating_rates(const int nonemptymgi, const double T_e, const doub
   ffheating = globals::ffheatingestimator[nonemptymgi];
 
   if constexpr (DIRECT_COL_HEAT) {
-    heatingcoolingrates->heating_collisional = C_deexc;
+    heatingcoolingrates.heating_collisional = C_deexc;
   } else {
     // Collisional heating (from estimators)
-    heatingcoolingrates->heating_collisional = globals::colheatingestimator[nonemptymgi];  // C_deexc + C_recomb;
+    heatingcoolingrates.heating_collisional = globals::colheatingestimator[nonemptymgi];  // C_deexc + C_recomb;
   }
 
-  heatingcoolingrates->heating_bf = bfheating;
-  heatingcoolingrates->heating_ff = ffheating;
+  heatingcoolingrates.heating_bf = bfheating;
+  heatingcoolingrates.heating_ff = ffheating;
 }
 
 // Thermal balance equation on which we have to iterate to get T_e
@@ -188,7 +188,7 @@ auto T_e_eqn_heating_minus_cooling(const double T_e, void *paras) -> double {
 
   const auto nonemptymgi = params->nonemptymgi;
   const double t_current = params->t_current;
-  auto *const heatingcoolingrates = params->heatingcoolingrates;
+  auto &heatingcoolingrates = *params->heatingcoolingrates;
   if constexpr (!LTEPOP_EXCITATION_USE_TJ) {
     if (std::abs((T_e / grid::get_Te(nonemptymgi)) - 1.) > 0.1) {
       grid::set_Te(nonemptymgi, T_e);
@@ -211,30 +211,30 @@ auto T_e_eqn_heating_minus_cooling(const double T_e, void *paras) -> double {
   const auto nne = grid::get_nne(nonemptymgi);
 
   // Then calculate heating and cooling rates
-  kpkt::calculate_cooling_rates(nonemptymgi, heatingcoolingrates);
+  kpkt::calculate_cooling_rates(nonemptymgi, &heatingcoolingrates);
   calculate_heating_rates(nonemptymgi, T_e, nne, heatingcoolingrates, *params->bfheatingcoeffs);
   const auto modelgridindex = grid::get_mgi_of_nonemptymgi(nonemptymgi);
 
   const auto ntlepton_frac_heating = nonthermal::get_nt_frac_heating(modelgridindex);
   const auto ntlepton_dep = nonthermal::get_deposition_rate_density(nonemptymgi);
   const auto ntalpha_frac_heating = 1.;
-  const auto ntalpha_dep = heatingcoolingrates->dep_alpha;
-  heatingcoolingrates->heating_dep = ntlepton_dep * ntlepton_frac_heating + ntalpha_dep * ntalpha_frac_heating;
-  heatingcoolingrates->dep_frac_heating =
-      (ntalpha_dep > 0) ? heatingcoolingrates->heating_dep / (ntlepton_dep + ntalpha_dep) : ntlepton_frac_heating;
+  const auto ntalpha_dep = heatingcoolingrates.dep_alpha;
+  heatingcoolingrates.heating_dep = ntlepton_dep * ntlepton_frac_heating + ntalpha_dep * ntalpha_frac_heating;
+  heatingcoolingrates.dep_frac_heating =
+      (ntalpha_dep > 0) ? heatingcoolingrates.heating_dep / (ntlepton_dep + ntalpha_dep) : ntlepton_frac_heating;
 
   // Adiabatic cooling term
   const double nntot = get_nnion_tot(nonemptymgi) + nne;
   const double p = nntot * KB * T_e;  // pressure in [erg/cm^3]
   const double volumetmin = grid::get_modelcell_assocvolume_tmin(modelgridindex);
   const double dV = 3 * volumetmin / pow(globals::tmin, 3) * pow(t_current, 2);  // really dV/dt
   const double V = volumetmin * pow(t_current / globals::tmin, 3);
-  heatingcoolingrates->cooling_adiabatic = p * dV / V;
+  heatingcoolingrates.cooling_adiabatic = p * dV / V;
 
-  const double total_heating_rate = heatingcoolingrates->heating_ff + heatingcoolingrates->heating_bf +
-                                    heatingcoolingrates->heating_collisional + heatingcoolingrates->heating_dep;
-  const double total_coolingrate = heatingcoolingrates->cooling_ff + heatingcoolingrates->cooling_fb +
-                                   heatingcoolingrates->cooling_collisional + heatingcoolingrates->cooling_adiabatic;
+  const double total_heating_rate = heatingcoolingrates.heating_ff + heatingcoolingrates.heating_bf +
+                                    heatingcoolingrates.heating_collisional + heatingcoolingrates.heating_dep;
+  const double total_coolingrate = heatingcoolingrates.cooling_ff + heatingcoolingrates.cooling_fb +
+                                   heatingcoolingrates.cooling_collisional + heatingcoolingrates.cooling_adiabatic;
 
   return total_heating_rate - total_coolingrate;
 }
@@ -312,14 +312,14 @@ void calculate_bfheatingcoeffs(int nonemptymgi, std::vector<double> &bfheatingco
 }
 
 void call_T_e_finder(const int nonemptymgi, const double t_current, const double T_min, const double T_max,
-                     HeatingCoolingRates *heatingcoolingrates, const std::vector<double> &bfheatingcoeffs) {
+                     HeatingCoolingRates &heatingcoolingrates, const std::vector<double> &bfheatingcoeffs) {
   const int modelgridindex = grid::get_mgi_of_nonemptymgi(nonemptymgi);
   const double T_e_old = grid::get_Te(nonemptymgi);
   printout("Finding T_e in cell %d at timestep %d...", modelgridindex, globals::timestep);
 
   TeSolutionParams paras = {.t_current = t_current,
                             .nonemptymgi = nonemptymgi,
-                            .heatingcoolingrates = heatingcoolingrates,
+                            .heatingcoolingrates = &heatingcoolingrates,
                             .bfheatingcoeffs = &bfheatingcoeffs};
 
   gsl_function find_T_e_f = {.function = &T_e_eqn_heating_minus_cooling, .params = &paras};

thermalbalance.h

@@ -25,7 +25,7 @@ struct HeatingCoolingRates {
 };
 
 void call_T_e_finder(int nonemptymgi, double t_current, double T_min, double T_max,
-                     HeatingCoolingRates *heatingcoolingrates, const std::vector<double> &bfheatingcoeffs);
+                     HeatingCoolingRates &heatingcoolingrates, const std::vector<double> &bfheatingcoeffs);
 [[nodiscard]] auto get_bfheatingcoeff_ana(int element, int ion, int level, int phixstargetindex, double T_R, double W)
     -> double;
 void calculate_bfheatingcoeffs(int nonemptymgi, std::vector<double> &bfheatingcoeffs);

update_grid.cc

@@ -32,7 +32,7 @@ namespace {
 std::vector<HeatingCoolingRates> heatingcoolingrates_thisrankcells;
 
 void write_to_estimators_file(FILE *estimators_file, const int mgi, const int timestep, const int titer,
-                              const HeatingCoolingRates *heatingcoolingrates) {
+                              const HeatingCoolingRates &heatingcoolingrates) {
   // return; disable for better performance (if estimators files are not needed)
 
   if (grid::get_numpropcells(mgi) < 1) {
@@ -653,17 +653,17 @@ void write_to_estimators_file(FILE *estimators_file, const int mgi, const int ti
   // ana means analytical at t_mid, i.e. the rates calculated from the nuclear abundances and decay data, not from
   // Monte Carlo
   fprintf(estimators_file, "emission_ana: gamma %11.5e positron %11.5e electron %11.5e alpha %11.5e\n",
-          heatingcoolingrates->eps_gamma_ana, heatingcoolingrates->eps_positron_ana,
-          heatingcoolingrates->eps_electron_ana, heatingcoolingrates->eps_alpha_ana);
+          heatingcoolingrates.eps_gamma_ana, heatingcoolingrates.eps_positron_ana, heatingcoolingrates.eps_electron_ana,
+          heatingcoolingrates.eps_alpha_ana);
   fprintf(estimators_file, "deposition: gamma %11.5e positron %11.5e electron %11.5e alpha %11.5e\n",
-          heatingcoolingrates->dep_gamma, heatingcoolingrates->dep_positron, heatingcoolingrates->dep_electron,
-          heatingcoolingrates->dep_alpha);
+          heatingcoolingrates.dep_gamma, heatingcoolingrates.dep_positron, heatingcoolingrates.dep_electron,
+          heatingcoolingrates.dep_alpha);
   fprintf(estimators_file, "heating: ff %11.5e bf %11.5e coll %11.5e       dep %11.5e heating_dep/total_dep %.3f\n",
-          heatingcoolingrates->heating_ff, heatingcoolingrates->heating_bf, heatingcoolingrates->heating_collisional,
-          heatingcoolingrates->heating_dep, heatingcoolingrates->dep_frac_heating);
+          heatingcoolingrates.heating_ff, heatingcoolingrates.heating_bf, heatingcoolingrates.heating_collisional,
+          heatingcoolingrates.heating_dep, heatingcoolingrates.dep_frac_heating);
   fprintf(estimators_file, "cooling: ff %11.5e fb %11.5e coll %11.5e adiabatic %11.5e\n",
-          heatingcoolingrates->cooling_ff, heatingcoolingrates->cooling_fb, heatingcoolingrates->cooling_collisional,
-          heatingcoolingrates->cooling_adiabatic);
+          heatingcoolingrates.cooling_ff, heatingcoolingrates.cooling_fb, heatingcoolingrates.cooling_collisional,
+          heatingcoolingrates.cooling_adiabatic);
 
   fprintf(estimators_file, "\n");
 
@@ -676,7 +676,7 @@ void write_to_estimators_file(FILE *estimators_file, const int mgi, const int ti
 }
 
 void solve_Te_nltepops(const int nonemptymgi, const int nts, const int nts_prev,
-                       HeatingCoolingRates *heatingcoolingrates)
+                       HeatingCoolingRates &heatingcoolingrates)
 // nts is the timestep number
 {
   const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);
@@ -873,7 +873,7 @@ static void titer_average_estimators(const int nonemptymgi) {
 #endif
 
 void update_grid_cell(const int mgi, const int nts, const int nts_prev, const int titer, const double tratmid,
-                      const double deltat, HeatingCoolingRates *heatingcoolingrates) {
+                      const double deltat, HeatingCoolingRates &heatingcoolingrates) {
   const ptrdiff_t nonemptymgi = grid::get_nonemptymgi_of_mgi(mgi);
 
   const double deltaV =
@@ -1127,19 +1127,17 @@ void update_grid(FILE *estimators_file, const int nts, const int nts_prev, const
 #endif
 
     for (int mgi = nstart; mgi < nstart + ndo; mgi++) {
-      // Check if this task should work on the current model grid cell.
-      // If yes, update the cell and write out the estimators
+      auto &heatingcoolingrates = heatingcoolingrates_thisrankcells.at(mgi - nstart);
       if (grid::get_numpropcells(mgi) > 0) {
-        update_grid_cell(mgi, nts, nts_prev, titer, tratmid, deltat,
-                         &heatingcoolingrates_thisrankcells.at(mgi - nstart));
+        update_grid_cell(mgi, nts, nts_prev, titer, tratmid, deltat, heatingcoolingrates);
       }
 
       // maybe want to add omp ordered here if the modelgrid cells should be output in order
 #ifdef _OPENMP
 #pragma omp critical(estimators_file)
 #endif
       {
-        write_to_estimators_file(estimators_file, mgi, nts, titer, &heatingcoolingrates_thisrankcells.at(mgi - nstart));
+        write_to_estimators_file(estimators_file, mgi, nts, titer, heatingcoolingrates);
       }
     }  // end parallel for loop over all modelgrid cells