From eb58112b6bf74cec217f181a1a6fdc341f58779c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 09:36:31 -0500 Subject: [PATCH 01/30] update grackle_rate_functions.h boilerplate (to generate doxygen docs) --- src/include/grackle_rate_functions.h | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/include/grackle_rate_functions.h b/src/include/grackle_rate_functions.h index 9d3647481..c882aad9f 100644 --- a/src/include/grackle_rate_functions.h +++ b/src/include/grackle_rate_functions.h @@ -1,4 +1,15 @@ -// Header file containing all rate function declarations. +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Declares rate functions +/// +//===----------------------------------------------------------------------===// + #include "grackle.h" #ifndef GRACKLE_RATE_FUNCTIONS_H From 75fd181d80b2f13a79e59c682faa8ce65ca9952d Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 09:48:14 -0500 Subject: [PATCH 02/30] add a docstring for regr_rate --- src/clib/rate_functions.c | 1 - src/include/grackle_rate_functions.h | 6 ++++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/clib/rate_functions.c b/src/clib/rate_functions.c index 154b82d31..e2bc10597 100644 --- a/src/clib/rate_functions.c +++ b/src/clib/rate_functions.c @@ -1421,7 +1421,6 @@ double gasGrain2_rate(double T, double units, chemistry_data *my_chemistry) } -//Calculation of regr. double regr_rate(double T, double units, chemistry_data *my_chemistry) { //(Equation 9, Wolfire et al., 1995) diff --git a/src/include/grackle_rate_functions.h b/src/include/grackle_rate_functions.h index c882aad9f..e97f2d552 100644 --- a/src/include/grackle_rate_functions.h +++ b/src/include/grackle_rate_functions.h @@ -101,7 +101,13 @@ double cieco_rate(double T, double units, chemistry_data *my_chemistry); double gasGrain_rate(double T, double units, chemistry_data *my_chemistry); double gasGrain2_rate(double T, double units, chemistry_data *my_chemistry); + +/// Calculate regr (as in GRain REcombination cooling) +/// +/// This is the value of a cluster of variables taken from equation 9 from +/// [Wolfire+95](https://ui.adsabs.harvard.edu/abs/1995ApJ...443..152W/abstract). double regr_rate(double T, double units, chemistry_data *my_chemistry); + double grain_growth_rate(double T, double units, chemistry_data *my_chemistry); double comp_rate(double units, chemistry_data *my_chemistry); From 143181065175a5faaa8ce7fb324dd92c25557d87 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 8 Jan 2026 23:24:19 -0500 Subject: [PATCH 03/30] first stab at extracting photoelectric heating --- src/clib/cool1d_multi_g.cpp | 107 ++++++++++++++++++++---------------- 1 file changed, 59 insertions(+), 48 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 02793cc90..00e3adab2 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -25,6 +25,61 @@ #include "internal_types.hpp" #include "utils-cpp.hpp" +void update_edot_photoelectric_heat( + double* edot, double* tgas, double* dust2gas, double* rhoH, + gr_mask_type* itmask, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, IndexRange idx_range, + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, int& i, + double dom_inv, double pe_eps, double pe_X, std::vector myisrf) { + if (my_chemistry->photoelectric_heating == 1) { + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + if (tgas[i] > 2.e4) { + cool1dmulti_buf.gammaha_eff[i] = 0.; + } else { + cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah; + } + } + } + + // Use eqn. 1 of Wolfire et al. (1995) + } else if (my_chemistry->photoelectric_heating == 2) { + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + if (tgas[i] > 2.e4) { + cool1dmulti_buf.gammaha_eff[i] = 0.; + } else { + // Assume constant epsilon = 0.05. + cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * 0.05 * myisrf[i]; + } + } + } + + // Full calculation of epsilon (eqn. 2 of Wolfire 1995) + } else if (my_chemistry->photoelectric_heating == 3) { + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + pe_X = + myisrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; + pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / + (1. + (pe_X / 5000.))); + cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * pe_eps * myisrf[i]; + } + } + } + + if (my_chemistry->photoelectric_heating > 0) { + for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + edot[i] = edot[i] + cool1dmulti_buf.gammaha_eff[i] * rhoH[i] * dom_inv * + dust2gas[i] / + my_chemistry->local_dust_to_gas_ratio; + } + } + } +} + void grackle::impl::cool1d_multi_g( int imetal, int iter, double* edot, double* tgas, double* mmw, double* p2d, double* tdust, double* metallicity, double* dust2gas, double* rhoH, @@ -1522,53 +1577,9 @@ void grackle::impl::cool1d_multi_g( // Photo-electric heating by UV-irradiated dust - if (my_chemistry->photoelectric_heating == 1) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - if (tgas[i] > 2.e4) { - cool1dmulti_buf.gammaha_eff[i] = 0.; - } else { - cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah; - } - } - } - - // Use eqn. 1 of Wolfire et al. (1995) - } else if (my_chemistry->photoelectric_heating == 2) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - if (tgas[i] > 2.e4) { - cool1dmulti_buf.gammaha_eff[i] = 0.; - } else { - // Assume constant epsilon = 0.05. - cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * 0.05 * myisrf[i]; - } - } - } - - // Full calculation of epsilon (eqn. 2 of Wolfire 1995) - } else if (my_chemistry->photoelectric_heating == 3) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - pe_X = - myisrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; - pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + - ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / - (1. + (pe_X / 5000.))); - cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * pe_eps * myisrf[i]; - } - } - } - - if (my_chemistry->photoelectric_heating > 0) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - edot[i] = edot[i] + cool1dmulti_buf.gammaha_eff[i] * rhoH[i] * dom_inv * - dust2gas[i] / - my_chemistry->local_dust_to_gas_ratio; - } - } - } + update_edot_photoelectric_heat( + edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates, idx_range, + cool1dmulti_buf, i, dom_inv, pe_eps, pe_X, myisrf); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) @@ -1969,4 +1980,4 @@ void grackle::impl::cool1d_multi_g( grackle::impl::drop_GrainSpeciesCollection(&gas_grainsp_heatrate); return; -} \ No newline at end of file +} From bb41cc7055f900fb20b3cb302cadcefb68adfb4c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 8 Jan 2026 23:30:30 -0500 Subject: [PATCH 04/30] a little work to clean up photoelectric heating logic --- src/clib/cool1d_multi_g.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 00e3adab2..3c788f3d4 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -26,13 +26,13 @@ #include "utils-cpp.hpp" void update_edot_photoelectric_heat( - double* edot, double* tgas, double* dust2gas, double* rhoH, - gr_mask_type* itmask, chemistry_data* my_chemistry, + double* edot, const double* tgas, const double* dust2gas, const double* rhoH, + const gr_mask_type* itmask, const chemistry_data* my_chemistry, chemistry_data_storage* my_rates, IndexRange idx_range, - grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, int& i, - double dom_inv, double pe_eps, double pe_X, std::vector myisrf) { + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, + double dom_inv, const double* myisrf) { if (my_chemistry->photoelectric_heating == 1) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { if (tgas[i] > 2.e4) { cool1dmulti_buf.gammaha_eff[i] = 0.; @@ -44,7 +44,7 @@ void update_edot_photoelectric_heat( // Use eqn. 1 of Wolfire et al. (1995) } else if (my_chemistry->photoelectric_heating == 2) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { if (tgas[i] > 2.e4) { cool1dmulti_buf.gammaha_eff[i] = 0.; @@ -57,11 +57,11 @@ void update_edot_photoelectric_heat( // Full calculation of epsilon (eqn. 2 of Wolfire 1995) } else if (my_chemistry->photoelectric_heating == 3) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { - pe_X = + double pe_X = myisrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; - pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + + double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / (1. + (pe_X / 5000.))); cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * pe_eps * myisrf[i]; @@ -70,7 +70,7 @@ void update_edot_photoelectric_heat( } if (my_chemistry->photoelectric_heating > 0) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { edot[i] = edot[i] + cool1dmulti_buf.gammaha_eff[i] * rhoH[i] * dom_inv * dust2gas[i] / @@ -1579,7 +1579,7 @@ void grackle::impl::cool1d_multi_g( update_edot_photoelectric_heat( edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates, idx_range, - cool1dmulti_buf, i, dom_inv, pe_eps, pe_X, myisrf); + cool1dmulti_buf, dom_inv, myisrf.data()); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) From c5e54bd8489494ac8a7369c49f96bc5abfd42d5e Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 8 Jan 2026 23:54:31 -0500 Subject: [PATCH 05/30] additional reformatting --- src/clib/cool1d_multi_g.cpp | 51 ++++++++++++++++++++++++++----------- 1 file changed, 36 insertions(+), 15 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 3c788f3d4..bbca7196a 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -25,19 +25,41 @@ #include "internal_types.hpp" #include "utils-cpp.hpp" +/// update edot, in place, with contributions from Photo-electric heating by +/// UV-irradiated dust +/// +/// Each of the 1D arrays is only valid for the specified @p idx_range +/// +/// @param [in,out] edot 1D array being used to accumulate the net rate of +/// change in thermal energy +/// @param[in] tgas 1D array of gas temperatures +/// @param[in] dust2gas 1D array of ratios between dust & gas densities +/// @param[in] rhoH 1D array holding the Hydrogen mass density +/// @param[in] itmask Specifies the general iteration-mask of the @p idx_range +/// for this calculation. +/// @param[in] my_chemistry holds a number of configuration parameters. +/// @param[in] gammah Parameterizes the calculation +/// @param[in] idx_range Specifies the current index-range +/// @param[in] cool1dmulti_buf Pre-allocated buffers that are used by this +/// function for scratch space (to hold a variety of quantities) +/// @param[in] dom_inv +/// @param[in] isrf 1D array specifying the strength of the interstellar +/// radiation field in Habing units +/// +/// @todo Perhaps we should extract the relevant parameters from my_chemistry? void update_edot_photoelectric_heat( - double* edot, const double* tgas, const double* dust2gas, const double* rhoH, - const gr_mask_type* itmask, const chemistry_data* my_chemistry, - chemistry_data_storage* my_rates, IndexRange idx_range, - grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, - double dom_inv, const double* myisrf) { + double* edot, const double* tgas, const double* dust2gas, + const double* rhoH, const gr_mask_type* itmask, + const chemistry_data* my_chemistry, double gammah, IndexRange idx_range, + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, + const double* isrf) { if (my_chemistry->photoelectric_heating == 1) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { if (tgas[i] > 2.e4) { cool1dmulti_buf.gammaha_eff[i] = 0.; } else { - cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah; + cool1dmulti_buf.gammaha_eff[i] = gammah; } } } @@ -50,7 +72,7 @@ void update_edot_photoelectric_heat( cool1dmulti_buf.gammaha_eff[i] = 0.; } else { // Assume constant epsilon = 0.05. - cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * 0.05 * myisrf[i]; + cool1dmulti_buf.gammaha_eff[i] = gammah * 0.05 * isrf[i]; } } } @@ -60,11 +82,11 @@ void update_edot_photoelectric_heat( for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { double pe_X = - myisrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; + isrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + - ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / - (1. + (pe_X / 5000.))); - cool1dmulti_buf.gammaha_eff[i] = my_rates->gammah * pe_eps * myisrf[i]; + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / + (1. + (pe_X / 5000.))); + cool1dmulti_buf.gammaha_eff[i] = gammah * pe_eps * isrf[i]; } } } @@ -1576,10 +1598,9 @@ void grackle::impl::cool1d_multi_g( } // Photo-electric heating by UV-irradiated dust - - update_edot_photoelectric_heat( - edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates, idx_range, - cool1dmulti_buf, dom_inv, myisrf.data()); + update_edot_photoelectric_heat(edot, tgas, dust2gas, rhoH, itmask, + my_chemistry, my_rates->gammah, idx_range, + cool1dmulti_buf, dom_inv, myisrf.data()); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) From a1edfb8426500b98a2be204690d39f772baa056f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 00:04:42 -0500 Subject: [PATCH 06/30] shift update_edot_photoelectric_heat to a separate file --- src/clib/CMakeLists.txt | 1 + src/clib/cool1d_multi_g.cpp | 85 ++------------------- src/clib/dust/gas_heating_cooling.hpp | 103 ++++++++++++++++++++++++++ 3 files changed, 109 insertions(+), 80 deletions(-) create mode 100644 src/clib/dust/gas_heating_cooling.hpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 7036b28e8..0f7cf58e5 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -109,6 +109,7 @@ add_library(Grackle_Grackle cool1d_multi_g.cpp cool1d_multi_g.hpp cool_multi_time_g.cpp cool_multi_time_g.h dust_props.hpp + dust/gas_heating_cooling.hpp dust/grain_species_info.cpp dust/grain_species_info.hpp dust/lookup_dust_rates1d.hpp init_misc_species_cool_rates.cpp init_misc_species_cool_rates.hpp diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index bbca7196a..d7677f2c2 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -13,6 +13,8 @@ // This file was initially generated automatically during conversion of the // cool1d_multi_g function from FORTRAN to C++ +#include "dust/gas_heating_cooling.hpp" + #include #include #include @@ -25,83 +27,6 @@ #include "internal_types.hpp" #include "utils-cpp.hpp" -/// update edot, in place, with contributions from Photo-electric heating by -/// UV-irradiated dust -/// -/// Each of the 1D arrays is only valid for the specified @p idx_range -/// -/// @param [in,out] edot 1D array being used to accumulate the net rate of -/// change in thermal energy -/// @param[in] tgas 1D array of gas temperatures -/// @param[in] dust2gas 1D array of ratios between dust & gas densities -/// @param[in] rhoH 1D array holding the Hydrogen mass density -/// @param[in] itmask Specifies the general iteration-mask of the @p idx_range -/// for this calculation. -/// @param[in] my_chemistry holds a number of configuration parameters. -/// @param[in] gammah Parameterizes the calculation -/// @param[in] idx_range Specifies the current index-range -/// @param[in] cool1dmulti_buf Pre-allocated buffers that are used by this -/// function for scratch space (to hold a variety of quantities) -/// @param[in] dom_inv -/// @param[in] isrf 1D array specifying the strength of the interstellar -/// radiation field in Habing units -/// -/// @todo Perhaps we should extract the relevant parameters from my_chemistry? -void update_edot_photoelectric_heat( - double* edot, const double* tgas, const double* dust2gas, - const double* rhoH, const gr_mask_type* itmask, - const chemistry_data* my_chemistry, double gammah, IndexRange idx_range, - grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, - const double* isrf) { - if (my_chemistry->photoelectric_heating == 1) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - if (tgas[i] > 2.e4) { - cool1dmulti_buf.gammaha_eff[i] = 0.; - } else { - cool1dmulti_buf.gammaha_eff[i] = gammah; - } - } - } - - // Use eqn. 1 of Wolfire et al. (1995) - } else if (my_chemistry->photoelectric_heating == 2) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - if (tgas[i] > 2.e4) { - cool1dmulti_buf.gammaha_eff[i] = 0.; - } else { - // Assume constant epsilon = 0.05. - cool1dmulti_buf.gammaha_eff[i] = gammah * 0.05 * isrf[i]; - } - } - } - - // Full calculation of epsilon (eqn. 2 of Wolfire 1995) - } else if (my_chemistry->photoelectric_heating == 3) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - double pe_X = - isrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; - double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + - ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / - (1. + (pe_X / 5000.))); - cool1dmulti_buf.gammaha_eff[i] = gammah * pe_eps * isrf[i]; - } - } - } - - if (my_chemistry->photoelectric_heating > 0) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - edot[i] = edot[i] + cool1dmulti_buf.gammaha_eff[i] * rhoH[i] * dom_inv * - dust2gas[i] / - my_chemistry->local_dust_to_gas_ratio; - } - } - } -} - void grackle::impl::cool1d_multi_g( int imetal, int iter, double* edot, double* tgas, double* mmw, double* p2d, double* tdust, double* metallicity, double* dust2gas, double* rhoH, @@ -1598,9 +1523,9 @@ void grackle::impl::cool1d_multi_g( } // Photo-electric heating by UV-irradiated dust - update_edot_photoelectric_heat(edot, tgas, dust2gas, rhoH, itmask, - my_chemistry, my_rates->gammah, idx_range, - cool1dmulti_buf, dom_inv, myisrf.data()); + dust::update_edot_photoelectric_heat( + edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates->gammah, + idx_range, cool1dmulti_buf, dom_inv, myisrf.data()); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp new file mode 100644 index 000000000..e3da28a31 --- /dev/null +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -0,0 +1,103 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// computes contributions to the heating/cooling of gas +/// +//===----------------------------------------------------------------------===// +#ifndef DUST_GAS_HEATING_COOLING_HPP +#define DUST_GAS_HEATING_COOLING_HPP + +#include "fortran_func_decls.h" // gr_mask_type +#include "grackle.h" +#include "index_helper.h" +#include "internal_types.hpp" + +#include + +namespace grackle::impl::dust { + +/// update edot, in place, with contributions from Photo-electric heating by +/// UV-irradiated dust +/// +/// Each of the 1D arrays is only valid for the specified @p idx_range +/// +/// @param [in,out] edot 1D array being used to accumulate the net rate of +/// change in thermal energy +/// @param[in] tgas 1D array of gas temperatures +/// @param[in] dust2gas 1D array of ratios between dust & gas densities +/// @param[in] rhoH 1D array holding the Hydrogen mass density +/// @param[in] itmask Specifies the general iteration-mask of the @p idx_range +/// for this calculation. +/// @param[in] my_chemistry holds a number of configuration parameters. +/// @param[in] gammah Parameterizes the calculation +/// @param[in] idx_range Specifies the current index-range +/// @param[in] cool1dmulti_buf Pre-allocated buffers that are used by this +/// function for scratch space (to hold a variety of quantities) +/// @param[in] dom_inv +/// @param[in] isrf 1D array specifying the strength of the interstellar +/// radiation field in Habing units +/// +/// @todo Perhaps we should extract the relevant parameters from my_chemistry? +void update_edot_photoelectric_heat( + double* edot, const double* tgas, const double* dust2gas, + const double* rhoH, const gr_mask_type* itmask, + const chemistry_data* my_chemistry, double gammah, IndexRange idx_range, + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, + const double* isrf) { + if (my_chemistry->photoelectric_heating == 1) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + if (tgas[i] > 2.e4) { + cool1dmulti_buf.gammaha_eff[i] = 0.; + } else { + cool1dmulti_buf.gammaha_eff[i] = gammah; + } + } + } + + // Use eqn. 1 of Wolfire et al. (1995) + } else if (my_chemistry->photoelectric_heating == 2) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + if (tgas[i] > 2.e4) { + cool1dmulti_buf.gammaha_eff[i] = 0.; + } else { + // Assume constant epsilon = 0.05. + cool1dmulti_buf.gammaha_eff[i] = gammah * 0.05 * isrf[i]; + } + } + } + + // Full calculation of epsilon (eqn. 2 of Wolfire 1995) + } else if (my_chemistry->photoelectric_heating == 3) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + double pe_X = + isrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; + double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / + (1. + (pe_X / 5000.))); + cool1dmulti_buf.gammaha_eff[i] = gammah * pe_eps * isrf[i]; + } + } + } + + if (my_chemistry->photoelectric_heating > 0) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + edot[i] = edot[i] + cool1dmulti_buf.gammaha_eff[i] * rhoH[i] * dom_inv * + dust2gas[i] / + my_chemistry->local_dust_to_gas_ratio; + } + } + } +} + +} // namespace grackle::impl::dust + +#endif // DUST_GAS_HEATING_COOLING_HPP \ No newline at end of file From 5c83471d95ba7b33f24e3096822d63621c0506f8 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 00:17:32 -0500 Subject: [PATCH 07/30] remove cool1dmulti_buf::gammaha_eff (1/2) --- src/clib/dust/gas_heating_cooling.hpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index e3da28a31..a06b9602f 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -43,12 +43,18 @@ namespace grackle::impl::dust { /// radiation field in Habing units /// /// @todo Perhaps we should extract the relevant parameters from my_chemistry? -void update_edot_photoelectric_heat( +inline void update_edot_photoelectric_heat( double* edot, const double* tgas, const double* dust2gas, const double* rhoH, const gr_mask_type* itmask, const chemistry_data* my_chemistry, double gammah, IndexRange idx_range, grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, const double* isrf) { + double local_dust_to_gas_ratio = my_chemistry->local_dust_to_gas_ratio; + auto update_edot_fn = [=](int i, double gammaha_eff) { + edot[i] = edot[i] + gammaha_eff * rhoH[i] * dom_inv * dust2gas[i] / + local_dust_to_gas_ratio; + }; + if (my_chemistry->photoelectric_heating == 1) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { @@ -57,6 +63,7 @@ void update_edot_photoelectric_heat( } else { cool1dmulti_buf.gammaha_eff[i] = gammah; } + update_edot_fn(i, cool1dmulti_buf.gammaha_eff[i]); } } @@ -70,6 +77,7 @@ void update_edot_photoelectric_heat( // Assume constant epsilon = 0.05. cool1dmulti_buf.gammaha_eff[i] = gammah * 0.05 * isrf[i]; } + update_edot_fn(i, cool1dmulti_buf.gammaha_eff[i]); } } @@ -83,16 +91,7 @@ void update_edot_photoelectric_heat( ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / (1. + (pe_X / 5000.))); cool1dmulti_buf.gammaha_eff[i] = gammah * pe_eps * isrf[i]; - } - } - } - - if (my_chemistry->photoelectric_heating > 0) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - edot[i] = edot[i] + cool1dmulti_buf.gammaha_eff[i] * rhoH[i] * dom_inv * - dust2gas[i] / - my_chemistry->local_dust_to_gas_ratio; + update_edot_fn(i, cool1dmulti_buf.gammaha_eff[i]); } } } From 73a6499e7567b2a1cfd52db4a7f0a447d7e9ceea Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 00:21:37 -0500 Subject: [PATCH 08/30] remove cool1dmulti_buf::gammaha_eff (2/2) --- src/clib/dust/gas_heating_cooling.hpp | 22 +++++++--------------- src/clib/internal_types.hpp | 2 -- 2 files changed, 7 insertions(+), 17 deletions(-) diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index a06b9602f..ef40d8645 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -58,12 +58,8 @@ inline void update_edot_photoelectric_heat( if (my_chemistry->photoelectric_heating == 1) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { - if (tgas[i] > 2.e4) { - cool1dmulti_buf.gammaha_eff[i] = 0.; - } else { - cool1dmulti_buf.gammaha_eff[i] = gammah; - } - update_edot_fn(i, cool1dmulti_buf.gammaha_eff[i]); + double gammaha_eff = (tgas[i] > 2.e4) ? 0. : gammah; + update_edot_fn(i, gammaha_eff); } } @@ -71,13 +67,9 @@ inline void update_edot_photoelectric_heat( } else if (my_chemistry->photoelectric_heating == 2) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { - if (tgas[i] > 2.e4) { - cool1dmulti_buf.gammaha_eff[i] = 0.; - } else { - // Assume constant epsilon = 0.05. - cool1dmulti_buf.gammaha_eff[i] = gammah * 0.05 * isrf[i]; - } - update_edot_fn(i, cool1dmulti_buf.gammaha_eff[i]); + // Assume constant epsilon = 0.05. + double gammaha_eff = (tgas[i] > 2.e4) ? 0. : gammah * 0.05 * isrf[i]; + update_edot_fn(i, gammaha_eff); } } @@ -90,8 +82,8 @@ inline void update_edot_photoelectric_heat( double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / (1. + (pe_X / 5000.))); - cool1dmulti_buf.gammaha_eff[i] = gammah * pe_eps * isrf[i]; - update_edot_fn(i, cool1dmulti_buf.gammaha_eff[i]); + double gammaha_eff = gammah * pe_eps * isrf[i]; + update_edot_fn(i, gammaha_eff); } } } diff --git a/src/clib/internal_types.hpp b/src/clib/internal_types.hpp index 7aaa30fa7..b5b150eaa 100644 --- a/src/clib/internal_types.hpp +++ b/src/clib/internal_types.hpp @@ -125,7 +125,6 @@ struct Cool1DMultiScratchBuf { double* tgasold = nullptr; double* mynh = nullptr; double* myde = nullptr; - double* gammaha_eff = nullptr; double* gasgr_tdust = nullptr; double* regr = nullptr; }; @@ -143,7 +142,6 @@ void visit_member_pair( f(VIS_MEMBER_NAME("tgasold"), obj0.tgasold, obj1.tgasold, vis::idx_range_len_multiple(1)); f(VIS_MEMBER_NAME("mynh"), obj0.mynh, obj1.mynh, vis::idx_range_len_multiple(1)); f(VIS_MEMBER_NAME("myde"), obj0.myde, obj1.myde, vis::idx_range_len_multiple(1)); - f(VIS_MEMBER_NAME("gammaha_eff"), obj0.gammaha_eff, obj1.gammaha_eff, vis::idx_range_len_multiple(1)); f(VIS_MEMBER_NAME("gasgr_tdust"), obj0.gasgr_tdust, obj1.gasgr_tdust, vis::idx_range_len_multiple(1)); f(VIS_MEMBER_NAME("regr"), obj0.regr, obj1.regr, vis::idx_range_len_multiple(1)); vis::end_visit(f); From c87c0d194c73c17a156615ba546a2445f9b886e2 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 00:33:31 -0500 Subject: [PATCH 09/30] clean the function --- src/clib/cool1d_multi_g.cpp | 6 +++--- src/clib/dust/gas_heating_cooling.hpp | 16 ++++++++++------ 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index d7677f2c2..1e7f8b14c 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -122,8 +122,8 @@ void grackle::impl::cool1d_multi_g( int i, iZscale, mycmbTfloor; double dom, qq, vibl, logtem0, logtem9, dlogtem, zr, hdlte1, hdlow1, gamma2, x, fudge, gphdl1, dom_inv, tau, ciefudge, coolunit, tbase1, nH2, nother, - nSSh, nratio, nssh_he, nratio_he, fSShHI, fSShHeI, pe_eps, pe_X, grbeta, - ih2cox, min_metallicity; + nSSh, nratio, nssh_he, nratio_he, fSShHI, fSShHeI, grbeta, ih2cox, + min_metallicity; double comp1, comp2; // Performing heap allocations for all of the subsequent buffers within this @@ -1525,7 +1525,7 @@ void grackle::impl::cool1d_multi_g( // Photo-electric heating by UV-irradiated dust dust::update_edot_photoelectric_heat( edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates->gammah, - idx_range, cool1dmulti_buf, dom_inv, myisrf.data()); + idx_range, cool1dmulti_buf.myde, dom_inv, myisrf.data()); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index ef40d8645..c9ad14560 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -36,20 +36,25 @@ namespace grackle::impl::dust { /// @param[in] my_chemistry holds a number of configuration parameters. /// @param[in] gammah Parameterizes the calculation /// @param[in] idx_range Specifies the current index-range -/// @param[in] cool1dmulti_buf Pre-allocated buffers that are used by this -/// function for scratch space (to hold a variety of quantities) +/// @param[in] e_density 1D array of electron mass densities (see below) /// @param[in] dom_inv /// @param[in] isrf 1D array specifying the strength of the interstellar /// radiation field in Habing units /// /// @todo Perhaps we should extract the relevant parameters from my_chemistry? +/// +/// @note +/// For primordial_chemistry > 0, @p e_density is typically just a copy of +/// the electron density field. Otherwise, these values are inferred inline void update_edot_photoelectric_heat( double* edot, const double* tgas, const double* dust2gas, const double* rhoH, const gr_mask_type* itmask, const chemistry_data* my_chemistry, double gammah, IndexRange idx_range, - grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, - const double* isrf) { + const double* e_density, double dom_inv, const double* isrf) { double local_dust_to_gas_ratio = my_chemistry->local_dust_to_gas_ratio; + + // define a lambda function to actually apply the update (once we determine + // the value of gammaha) auto update_edot_fn = [=](int i, double gammaha_eff) { edot[i] = edot[i] + gammaha_eff * rhoH[i] * dom_inv * dust2gas[i] / local_dust_to_gas_ratio; @@ -77,8 +82,7 @@ inline void update_edot_photoelectric_heat( } else if (my_chemistry->photoelectric_heating == 3) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { - double pe_X = - isrf[i] * dom_inv * std::sqrt(tgas[i]) / cool1dmulti_buf.myde[i]; + double pe_X = isrf[i] * dom_inv * std::sqrt(tgas[i]) / e_density[i]; double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / (1. + (pe_X / 5000.))); From eed121f4d7b4d1108d685fbb3ab92f29bf490d87 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 07:51:03 -0500 Subject: [PATCH 10/30] initial attempt to factor out dust recombination cooling --- src/clib/cool1d_multi_g.cpp | 30 ++++-------------------- src/clib/dust/gas_heating_cooling.hpp | 33 +++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 26 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 1e7f8b14c..046f5f911 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -122,7 +122,7 @@ void grackle::impl::cool1d_multi_g( int i, iZscale, mycmbTfloor; double dom, qq, vibl, logtem0, logtem9, dlogtem, zr, hdlte1, hdlow1, gamma2, x, fudge, gphdl1, dom_inv, tau, ciefudge, coolunit, tbase1, nH2, nother, - nSSh, nratio, nssh_he, nratio_he, fSShHI, fSShHeI, grbeta, ih2cox, + nSSh, nratio, nssh_he, nratio_he, fSShHI, fSShHeI, ih2cox, min_metallicity; double comp1, comp2; @@ -1528,31 +1528,9 @@ void grackle::impl::cool1d_multi_g( idx_range, cool1dmulti_buf.myde, dom_inv, myisrf.data()); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) - - if ((my_chemistry->dust_chemistry > 0) || - (my_chemistry->dust_recombination_cooling > 0)) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - cool1dmulti_buf.regr[i] = - my_rates->regr[logTlininterp_buf.indixe[i] - 1] + - logTlininterp_buf.tdef[i] * - (my_rates->regr[logTlininterp_buf.indixe[i] + 1 - 1] - - my_rates->regr[logTlininterp_buf.indixe[i] - 1]); - } - } - - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - grbeta = 0.74 / std::pow(tgas[i], 0.068); - edot[i] = edot[i] - - cool1dmulti_buf.regr[i] * - std::pow((myisrf[i] * dom_inv / cool1dmulti_buf.myde[i]), - grbeta) * - cool1dmulti_buf.myde[i] * rhoH[i] * dust2gas[i] / - my_chemistry->local_dust_to_gas_ratio; - } - } - } + dust::update_edot_dust_recombination( + edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates, idx_range, + logTlininterp_buf, cool1dmulti_buf, dom_inv, myisrf.data()); // Compton cooling or heating and X-ray compton heating diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index c9ad14560..e7265b397 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -93,6 +93,39 @@ inline void update_edot_photoelectric_heat( } } +void update_edot_dust_recombination( + double* edot, double* tgas, double* dust2gas, double* rhoH, + gr_mask_type* itmask, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, IndexRange idx_range, + grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, + grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, + double* isrf) { + if ((my_chemistry->dust_chemistry > 0) || + (my_chemistry->dust_recombination_cooling > 0)) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + cool1dmulti_buf.regr[i] = + my_rates->regr[logTlininterp_buf.indixe[i] - 1] + + logTlininterp_buf.tdef[i] * + (my_rates->regr[logTlininterp_buf.indixe[i] + 1 - 1] - + my_rates->regr[logTlininterp_buf.indixe[i] - 1]); + } + } + + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + double grbeta = 0.74 / std::pow(tgas[i], 0.068); + edot[i] = edot[i] - + cool1dmulti_buf.regr[i] * + std::pow((isrf[i] * dom_inv / cool1dmulti_buf.myde[i]), + grbeta) * + cool1dmulti_buf.myde[i] * rhoH[i] * dust2gas[i] / + my_chemistry->local_dust_to_gas_ratio; + } + } + } +} + } // namespace grackle::impl::dust #endif // DUST_GAS_HEATING_COOLING_HPP \ No newline at end of file From a354f5b91803841403345be67f018d9550060323 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 07:59:38 -0500 Subject: [PATCH 11/30] start cleaning up dust-recombination-cooling --- src/clib/cool1d_multi_g.cpp | 11 +++++-- src/clib/dust/gas_heating_cooling.hpp | 46 +++++++++++++-------------- 2 files changed, 30 insertions(+), 27 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 046f5f911..e38ef67fb 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1528,9 +1528,14 @@ void grackle::impl::cool1d_multi_g( idx_range, cool1dmulti_buf.myde, dom_inv, myisrf.data()); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) - dust::update_edot_dust_recombination( - edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates, idx_range, - logTlininterp_buf, cool1dmulti_buf, dom_inv, myisrf.data()); + + if ((my_chemistry->dust_chemistry > 0) || + (my_chemistry->dust_recombination_cooling > 0)) { + dust::update_edot_dust_recombination( + edot, tgas, dust2gas, rhoH, itmask, + my_chemistry->local_dust_to_gas_ratio, my_rates, idx_range, + logTlininterp_buf, cool1dmulti_buf, dom_inv, myisrf.data()); + } // Compton cooling or heating and X-ray compton heating diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index e7265b397..5db8a858b 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -94,34 +94,32 @@ inline void update_edot_photoelectric_heat( } void update_edot_dust_recombination( - double* edot, double* tgas, double* dust2gas, double* rhoH, - gr_mask_type* itmask, chemistry_data* my_chemistry, - chemistry_data_storage* my_rates, IndexRange idx_range, + double* edot, const double* tgas, const double* dust2gas, + const double* rhoH, const gr_mask_type* itmask, + double local_dust_to_gas_ratio, const chemistry_data_storage* my_rates, + IndexRange idx_range, grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, - double* isrf) { - if ((my_chemistry->dust_chemistry > 0) || - (my_chemistry->dust_recombination_cooling > 0)) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - cool1dmulti_buf.regr[i] = - my_rates->regr[logTlininterp_buf.indixe[i] - 1] + - logTlininterp_buf.tdef[i] * - (my_rates->regr[logTlininterp_buf.indixe[i] + 1 - 1] - - my_rates->regr[logTlininterp_buf.indixe[i] - 1]); - } + const double* isrf) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + cool1dmulti_buf.regr[i] = + my_rates->regr[logTlininterp_buf.indixe[i] - 1] + + logTlininterp_buf.tdef[i] * + (my_rates->regr[logTlininterp_buf.indixe[i] + 1 - 1] - + my_rates->regr[logTlininterp_buf.indixe[i] - 1]); } + } - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { - double grbeta = 0.74 / std::pow(tgas[i], 0.068); - edot[i] = edot[i] - - cool1dmulti_buf.regr[i] * - std::pow((isrf[i] * dom_inv / cool1dmulti_buf.myde[i]), - grbeta) * - cool1dmulti_buf.myde[i] * rhoH[i] * dust2gas[i] / - my_chemistry->local_dust_to_gas_ratio; - } + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask[i] != MASK_FALSE) { + double grbeta = 0.74 / std::pow(tgas[i], 0.068); + edot[i] = + edot[i] - + cool1dmulti_buf.regr[i] * + std::pow((isrf[i] * dom_inv / cool1dmulti_buf.myde[i]), grbeta) * + cool1dmulti_buf.myde[i] * rhoH[i] * dust2gas[i] / + local_dust_to_gas_ratio; } } } From ae3a508a626ed49299088c27414c0cc759f6abb2 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 08:07:16 -0500 Subject: [PATCH 12/30] dust_recombination_cooling: stop tracking an array of local regr values --- src/clib/cool1d_multi_g.cpp | 2 +- src/clib/dust/gas_heating_cooling.hpp | 17 +++++------------ src/clib/internal_types.hpp | 2 -- 3 files changed, 6 insertions(+), 15 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index e38ef67fb..c2bc47284 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1534,7 +1534,7 @@ void grackle::impl::cool1d_multi_g( dust::update_edot_dust_recombination( edot, tgas, dust2gas, rhoH, itmask, my_chemistry->local_dust_to_gas_ratio, my_rates, idx_range, - logTlininterp_buf, cool1dmulti_buf, dom_inv, myisrf.data()); + logTlininterp_buf, cool1dmulti_buf.myde, dom_inv, myisrf.data()); } // Compton cooling or heating and X-ray compton heating diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index 5db8a858b..e41c5be3c 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -93,33 +93,26 @@ inline void update_edot_photoelectric_heat( } } -void update_edot_dust_recombination( +inline void update_edot_dust_recombination( double* edot, const double* tgas, const double* dust2gas, const double* rhoH, const gr_mask_type* itmask, double local_dust_to_gas_ratio, const chemistry_data_storage* my_rates, IndexRange idx_range, grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, - grackle::impl::Cool1DMultiScratchBuf cool1dmulti_buf, double dom_inv, - const double* isrf) { + const double* e_density, double dom_inv, const double* isrf) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { - cool1dmulti_buf.regr[i] = + double cur_regr = my_rates->regr[logTlininterp_buf.indixe[i] - 1] + logTlininterp_buf.tdef[i] * (my_rates->regr[logTlininterp_buf.indixe[i] + 1 - 1] - my_rates->regr[logTlininterp_buf.indixe[i] - 1]); - } - } - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask[i] != MASK_FALSE) { double grbeta = 0.74 / std::pow(tgas[i], 0.068); edot[i] = edot[i] - - cool1dmulti_buf.regr[i] * - std::pow((isrf[i] * dom_inv / cool1dmulti_buf.myde[i]), grbeta) * - cool1dmulti_buf.myde[i] * rhoH[i] * dust2gas[i] / - local_dust_to_gas_ratio; + cur_regr * std::pow((isrf[i] * dom_inv / e_density[i]), grbeta) * + e_density[i] * rhoH[i] * dust2gas[i] / local_dust_to_gas_ratio; } } } diff --git a/src/clib/internal_types.hpp b/src/clib/internal_types.hpp index b5b150eaa..618ab8346 100644 --- a/src/clib/internal_types.hpp +++ b/src/clib/internal_types.hpp @@ -126,7 +126,6 @@ struct Cool1DMultiScratchBuf { double* mynh = nullptr; double* myde = nullptr; double* gasgr_tdust = nullptr; - double* regr = nullptr; }; /// used to help implement the visitor design pattern @@ -143,7 +142,6 @@ void visit_member_pair( f(VIS_MEMBER_NAME("mynh"), obj0.mynh, obj1.mynh, vis::idx_range_len_multiple(1)); f(VIS_MEMBER_NAME("myde"), obj0.myde, obj1.myde, vis::idx_range_len_multiple(1)); f(VIS_MEMBER_NAME("gasgr_tdust"), obj0.gasgr_tdust, obj1.gasgr_tdust, vis::idx_range_len_multiple(1)); - f(VIS_MEMBER_NAME("regr"), obj0.regr, obj1.regr, vis::idx_range_len_multiple(1)); vis::end_visit(f); } From cc4c8c2e2646d0dcec66d20cb98245535236b193 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 08:20:55 -0500 Subject: [PATCH 13/30] cleanup and document update_edot_dust_recombination --- src/clib/cool1d_multi_g.cpp | 9 +++--- src/clib/dust/gas_heating_cooling.hpp | 45 ++++++++++++++++++++++----- 2 files changed, 42 insertions(+), 12 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index c2bc47284..2bde2778f 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1531,10 +1531,11 @@ void grackle::impl::cool1d_multi_g( if ((my_chemistry->dust_chemistry > 0) || (my_chemistry->dust_recombination_cooling > 0)) { - dust::update_edot_dust_recombination( - edot, tgas, dust2gas, rhoH, itmask, - my_chemistry->local_dust_to_gas_ratio, my_rates, idx_range, - logTlininterp_buf, cool1dmulti_buf.myde, dom_inv, myisrf.data()); + dust::update_edot_dust_recombination(edot, tgas, dust2gas, rhoH, itmask, + my_chemistry->local_dust_to_gas_ratio, + logTlininterp_buf, my_rates->regr, + idx_range, cool1dmulti_buf.myde, + dom_inv, myisrf.data()); } // Compton cooling or heating and X-ray compton heating diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index e41c5be3c..2ff70f347 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -93,20 +93,49 @@ inline void update_edot_photoelectric_heat( } } +/// update edot, in place, with contributions from electron recombinations +/// onto dust +/// +/// Each of the 1D arrays (other than @p regr) is only valid for the +/// specified @p idx_range +/// +/// @param [in,out] edot 1D array being used to accumulate the net rate of +/// change in thermal energy +/// @param[in] tgas 1D array of gas temperatures +/// @param[in] dust2gas 1D array of ratios between dust & gas densities +/// @param[in] rhoH 1D array holding the Hydrogen mass density +/// @param[in] itmask Specifies the general iteration-mask of the @p idx_range +/// for this calculation. +/// @param[in] local_dust_to_gas_ratio ratio of total dust mass to gas mass +/// in the local Universe +/// @param[in] logTlininterp_buf Hold pre-computed values for each location +/// in @p idx_range with values that are used to linearly interpolate +/// tables ov values with respect to the natural log of @p tgas1. Used +/// here to interpolate the @p regr table. +/// @param[in] regr 1D table of rate values (precomputed on the shared gas +/// temperature grid) that parameterize the calculation. +/// @param[in] idx_range Specifies the current index-range +/// @param[in] e_density 1D array of electron mass densities (see below) +/// @param[in] dom_inv +/// @param[in] isrf 1D array specifying the strength of the interstellar +/// radiation field in Habing units +/// +/// @note +/// For primordial_chemistry > 0, @p e_density is typically just a copy of +/// the electron density field. Otherwise, these values are inferred inline void update_edot_dust_recombination( double* edot, const double* tgas, const double* dust2gas, const double* rhoH, const gr_mask_type* itmask, - double local_dust_to_gas_ratio, const chemistry_data_storage* my_rates, - IndexRange idx_range, + double local_dust_to_gas_ratio, grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, - const double* e_density, double dom_inv, const double* isrf) { + const double* regr, IndexRange idx_range, const double* e_density, + double dom_inv, const double* isrf) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { - double cur_regr = - my_rates->regr[logTlininterp_buf.indixe[i] - 1] + - logTlininterp_buf.tdef[i] * - (my_rates->regr[logTlininterp_buf.indixe[i] + 1 - 1] - - my_rates->regr[logTlininterp_buf.indixe[i] - 1]); + double cur_regr = regr[logTlininterp_buf.indixe[i] - 1] + + logTlininterp_buf.tdef[i] * + (regr[logTlininterp_buf.indixe[i] + 1 - 1] - + regr[logTlininterp_buf.indixe[i] - 1]); double grbeta = 0.74 / std::pow(tgas[i], 0.068); edot[i] = From 516e37493459a02ddbba2dabc43165acf0554582 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 08:30:28 -0500 Subject: [PATCH 14/30] slightly shuffle the order of arguments in update_edot_photoelectric_heat and update_edot_dust_recombination --- src/clib/cool1d_multi_g.cpp | 13 ++++++------- src/clib/dust/gas_heating_cooling.hpp | 27 +++++++++++++-------------- 2 files changed, 19 insertions(+), 21 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 2bde2778f..fa0c4bf3e 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1524,18 +1524,17 @@ void grackle::impl::cool1d_multi_g( // Photo-electric heating by UV-irradiated dust dust::update_edot_photoelectric_heat( - edot, tgas, dust2gas, rhoH, itmask, my_chemistry, my_rates->gammah, - idx_range, cool1dmulti_buf.myde, dom_inv, myisrf.data()); + edot, tgas, dust2gas, rhoH, cool1dmulti_buf.myde, myisrf.data(), itmask, + my_chemistry, my_rates->gammah, idx_range, dom_inv); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) if ((my_chemistry->dust_chemistry > 0) || (my_chemistry->dust_recombination_cooling > 0)) { - dust::update_edot_dust_recombination(edot, tgas, dust2gas, rhoH, itmask, - my_chemistry->local_dust_to_gas_ratio, - logTlininterp_buf, my_rates->regr, - idx_range, cool1dmulti_buf.myde, - dom_inv, myisrf.data()); + dust::update_edot_dust_recombination( + edot, tgas, dust2gas, rhoH, cool1dmulti_buf.myde, myisrf.data(), itmask, + my_chemistry->local_dust_to_gas_ratio, logTlininterp_buf, + my_rates->regr, idx_range, dom_inv); } // Compton cooling or heating and X-ray compton heating diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index 2ff70f347..3ac57c885 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -31,15 +31,15 @@ namespace grackle::impl::dust { /// @param[in] tgas 1D array of gas temperatures /// @param[in] dust2gas 1D array of ratios between dust & gas densities /// @param[in] rhoH 1D array holding the Hydrogen mass density +/// @param[in] e_density 1D array of electron mass densities (see below) +/// @param[in] isrf 1D array specifying the strength of the interstellar +/// radiation field in Habing units /// @param[in] itmask Specifies the general iteration-mask of the @p idx_range /// for this calculation. /// @param[in] my_chemistry holds a number of configuration parameters. /// @param[in] gammah Parameterizes the calculation /// @param[in] idx_range Specifies the current index-range -/// @param[in] e_density 1D array of electron mass densities (see below) /// @param[in] dom_inv -/// @param[in] isrf 1D array specifying the strength of the interstellar -/// radiation field in Habing units /// /// @todo Perhaps we should extract the relevant parameters from my_chemistry? /// @@ -48,9 +48,9 @@ namespace grackle::impl::dust { /// the electron density field. Otherwise, these values are inferred inline void update_edot_photoelectric_heat( double* edot, const double* tgas, const double* dust2gas, - const double* rhoH, const gr_mask_type* itmask, - const chemistry_data* my_chemistry, double gammah, IndexRange idx_range, - const double* e_density, double dom_inv, const double* isrf) { + const double* rhoH, const double* e_density, const double* isrf, + const gr_mask_type* itmask, const chemistry_data* my_chemistry, + double gammah, IndexRange idx_range, double dom_inv) { double local_dust_to_gas_ratio = my_chemistry->local_dust_to_gas_ratio; // define a lambda function to actually apply the update (once we determine @@ -104,6 +104,9 @@ inline void update_edot_photoelectric_heat( /// @param[in] tgas 1D array of gas temperatures /// @param[in] dust2gas 1D array of ratios between dust & gas densities /// @param[in] rhoH 1D array holding the Hydrogen mass density +/// @param[in] e_density 1D array of electron mass densities (see below) +/// @param[in] isrf 1D array specifying the strength of the interstellar +/// radiation field in Habing units /// @param[in] itmask Specifies the general iteration-mask of the @p idx_range /// for this calculation. /// @param[in] local_dust_to_gas_ratio ratio of total dust mass to gas mass @@ -115,21 +118,17 @@ inline void update_edot_photoelectric_heat( /// @param[in] regr 1D table of rate values (precomputed on the shared gas /// temperature grid) that parameterize the calculation. /// @param[in] idx_range Specifies the current index-range -/// @param[in] e_density 1D array of electron mass densities (see below) /// @param[in] dom_inv -/// @param[in] isrf 1D array specifying the strength of the interstellar -/// radiation field in Habing units /// /// @note /// For primordial_chemistry > 0, @p e_density is typically just a copy of -/// the electron density field. Otherwise, these values are inferred +/// the electron density field. Otherwise, these values are inferred. inline void update_edot_dust_recombination( double* edot, const double* tgas, const double* dust2gas, - const double* rhoH, const gr_mask_type* itmask, - double local_dust_to_gas_ratio, + const double* rhoH, const double* e_density, const double* isrf, + const gr_mask_type* itmask, double local_dust_to_gas_ratio, grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, - const double* regr, IndexRange idx_range, const double* e_density, - double dom_inv, const double* isrf) { + const double* regr, IndexRange idx_range, double dom_inv) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask[i] != MASK_FALSE) { double cur_regr = regr[logTlininterp_buf.indixe[i] - 1] + From 74073989a6fa92207badc3fab94fb077d21e37bd Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 08:37:10 -0500 Subject: [PATCH 15/30] update_edot_dust_recombination: indexing & docstring --- src/clib/cool1d_multi_g.cpp | 1 - src/clib/dust/gas_heating_cooling.hpp | 24 ++++++++++++++---------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index fa0c4bf3e..e354424f1 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1528,7 +1528,6 @@ void grackle::impl::cool1d_multi_g( my_chemistry, my_rates->gammah, idx_range, dom_inv); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) - if ((my_chemistry->dust_chemistry > 0) || (my_chemistry->dust_recombination_cooling > 0)) { dust::update_edot_dust_recombination( diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index 3ac57c885..0515f08f2 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -93,14 +93,17 @@ inline void update_edot_photoelectric_heat( } } -/// update edot, in place, with contributions from electron recombinations -/// onto dust +/// update edot, in place, with contributions from electron recombination +/// onto positively charged dust grains. +/// +/// Based upon equation 9 of +/// [Wolfire+95](https://ui.adsabs.harvard.edu/abs/1995ApJ...443..152W/abstract). /// /// Each of the 1D arrays (other than @p regr) is only valid for the /// specified @p idx_range /// /// @param [in,out] edot 1D array being used to accumulate the net rate of -/// change in thermal energy +/// change in internal energy of the gas /// @param[in] tgas 1D array of gas temperatures /// @param[in] dust2gas 1D array of ratios between dust & gas densities /// @param[in] rhoH 1D array holding the Hydrogen mass density @@ -113,10 +116,11 @@ inline void update_edot_photoelectric_heat( /// in the local Universe /// @param[in] logTlininterp_buf Hold pre-computed values for each location /// in @p idx_range with values that are used to linearly interpolate -/// tables ov values with respect to the natural log of @p tgas1. Used +/// tables of values with respect to the natural log of @p tgas1. Used /// here to interpolate the @p regr table. /// @param[in] regr 1D table of rate values (precomputed on the shared gas -/// temperature grid) that parameterize the calculation. +/// temperature grid using @ref regr_rate). This is just a cluster of +/// variables taken from the equation. /// @param[in] idx_range Specifies the current index-range /// @param[in] dom_inv /// @@ -129,12 +133,12 @@ inline void update_edot_dust_recombination( const gr_mask_type* itmask, double local_dust_to_gas_ratio, grackle::impl::LogTLinInterpScratchBuf logTlininterp_buf, const double* regr, IndexRange idx_range, double dom_inv) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { if (itmask[i] != MASK_FALSE) { - double cur_regr = regr[logTlininterp_buf.indixe[i] - 1] + - logTlininterp_buf.tdef[i] * - (regr[logTlininterp_buf.indixe[i] + 1 - 1] - - regr[logTlininterp_buf.indixe[i] - 1]); + double cur_regr = + regr[logTlininterp_buf.indixe[i] - 1] + + logTlininterp_buf.tdef[i] * (regr[logTlininterp_buf.indixe[i]] - + regr[logTlininterp_buf.indixe[i] - 1]); double grbeta = 0.74 / std::pow(tgas[i], 0.068); edot[i] = From 1625e7b7664ee1e875d2dd094aa095871d282047 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 09:06:34 -0500 Subject: [PATCH 16/30] update_edot_photoelectric_heat adjust inner loops --- src/clib/dust/gas_heating_cooling.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heating_cooling.hpp index 0515f08f2..2f781967d 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heating_cooling.hpp @@ -61,7 +61,7 @@ inline void update_edot_photoelectric_heat( }; if (my_chemistry->photoelectric_heating == 1) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { if (itmask[i] != MASK_FALSE) { double gammaha_eff = (tgas[i] > 2.e4) ? 0. : gammah; update_edot_fn(i, gammaha_eff); @@ -70,7 +70,7 @@ inline void update_edot_photoelectric_heat( // Use eqn. 1 of Wolfire et al. (1995) } else if (my_chemistry->photoelectric_heating == 2) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { if (itmask[i] != MASK_FALSE) { // Assume constant epsilon = 0.05. double gammaha_eff = (tgas[i] > 2.e4) ? 0. : gammah * 0.05 * isrf[i]; @@ -80,7 +80,7 @@ inline void update_edot_photoelectric_heat( // Full calculation of epsilon (eqn. 2 of Wolfire 1995) } else if (my_chemistry->photoelectric_heating == 3) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { if (itmask[i] != MASK_FALSE) { double pe_X = isrf[i] * dom_inv * std::sqrt(tgas[i]) / e_density[i]; double pe_eps = (4.9e-2 / (1. + std::pow((pe_X / 1925.), 0.73))) + From f086a4442dcc96623610bd817a0e58cc730a4744 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 08:57:02 -0500 Subject: [PATCH 17/30] introduce dust_gas_edot namespace --- src/clib/CMakeLists.txt | 2 +- src/clib/cool1d_multi_g.cpp | 6 +++--- .../{gas_heating_cooling.hpp => gas_heat_cool.hpp} | 14 +++++++++----- 3 files changed, 13 insertions(+), 9 deletions(-) rename src/clib/dust/{gas_heating_cooling.hpp => gas_heat_cool.hpp} (95%) diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 0f7cf58e5..976d52c87 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -109,7 +109,7 @@ add_library(Grackle_Grackle cool1d_multi_g.cpp cool1d_multi_g.hpp cool_multi_time_g.cpp cool_multi_time_g.h dust_props.hpp - dust/gas_heating_cooling.hpp + dust/gas_heat_cool.hpp dust/grain_species_info.cpp dust/grain_species_info.hpp dust/lookup_dust_rates1d.hpp init_misc_species_cool_rates.cpp init_misc_species_cool_rates.hpp diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index e354424f1..dc0af4e01 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -13,7 +13,7 @@ // This file was initially generated automatically during conversion of the // cool1d_multi_g function from FORTRAN to C++ -#include "dust/gas_heating_cooling.hpp" +#include "dust/gas_heat_cool.hpp" #include #include @@ -1523,14 +1523,14 @@ void grackle::impl::cool1d_multi_g( } // Photo-electric heating by UV-irradiated dust - dust::update_edot_photoelectric_heat( + dust_gas_edot::update_edot_photoelectric_heat( edot, tgas, dust2gas, rhoH, cool1dmulti_buf.myde, myisrf.data(), itmask, my_chemistry, my_rates->gammah, idx_range, dom_inv); // Electron recombination onto dust grains (eqn. 9 of Wolfire 1995) if ((my_chemistry->dust_chemistry > 0) || (my_chemistry->dust_recombination_cooling > 0)) { - dust::update_edot_dust_recombination( + dust_gas_edot::update_edot_dust_recombination( edot, tgas, dust2gas, rhoH, cool1dmulti_buf.myde, myisrf.data(), itmask, my_chemistry->local_dust_to_gas_ratio, logTlininterp_buf, my_rates->regr, idx_range, dom_inv); diff --git a/src/clib/dust/gas_heating_cooling.hpp b/src/clib/dust/gas_heat_cool.hpp similarity index 95% rename from src/clib/dust/gas_heating_cooling.hpp rename to src/clib/dust/gas_heat_cool.hpp index 2f781967d..f24c254a0 100644 --- a/src/clib/dust/gas_heating_cooling.hpp +++ b/src/clib/dust/gas_heat_cool.hpp @@ -9,8 +9,8 @@ /// computes contributions to the heating/cooling of gas /// //===----------------------------------------------------------------------===// -#ifndef DUST_GAS_HEATING_COOLING_HPP -#define DUST_GAS_HEATING_COOLING_HPP +#ifndef DUST_GAS_HEAT_COOL_HPP +#define DUST_GAS_HEAT_COOL_HPP #include "fortran_func_decls.h" // gr_mask_type #include "grackle.h" @@ -19,7 +19,11 @@ #include -namespace grackle::impl::dust { +/// this namespace groups logic modeling the various sources of gas +/// heating/cooling from dust-chemistry +/// +/// @todo Come up with a better name? +namespace grackle::impl::dust_gas_edot { /// update edot, in place, with contributions from Photo-electric heating by /// UV-irradiated dust @@ -149,6 +153,6 @@ inline void update_edot_dust_recombination( } } -} // namespace grackle::impl::dust +} // namespace grackle::impl::dust_gas_edot -#endif // DUST_GAS_HEATING_COOLING_HPP \ No newline at end of file +#endif // DUST_GAS_HEAT_COOL_HPP \ No newline at end of file From 4229e6e0b4326600c842f9a523dbd43556a78d59 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 10:14:52 -0500 Subject: [PATCH 18/30] add docstring for gammah_rate --- src/clib/rate_functions.c | 3 +-- src/include/grackle_rate_functions.h | 7 +++++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/clib/rate_functions.c b/src/clib/rate_functions.c index e2bc10597..1b58caa72 100644 --- a/src/clib/rate_functions.c +++ b/src/clib/rate_functions.c @@ -1436,11 +1436,10 @@ double comp_rate(double units, chemistry_data *my_chemistry) return 5.65e-36 / units; } -//Calculation of gammah. double gammah_rate(double units, chemistry_data *my_chemistry) { - //Default is 8.5e-26 for epsilon=0.05, G_0=1.7 (rate in erg s^-1 cm^-3). if (my_chemistry->photoelectric_heating <= 1) { + // Default: 8.5e-26 for epsilon=0.05, G_0=1.7 (rate in erg s^-1 cm^-3) return my_chemistry->photoelectric_heating_rate / units; } else { //photoelectric_heating set to 2 or 3 //User to specify G_0, epsilon set to 0.05 or calculated directly. diff --git a/src/include/grackle_rate_functions.h b/src/include/grackle_rate_functions.h index e97f2d552..dfe17e4f8 100644 --- a/src/include/grackle_rate_functions.h +++ b/src/include/grackle_rate_functions.h @@ -111,6 +111,13 @@ double regr_rate(double T, double units, chemistry_data *my_chemistry); double grain_growth_rate(double T, double units, chemistry_data *my_chemistry); double comp_rate(double units, chemistry_data *my_chemistry); + +/// Calculate quantity closely related to Γ for photo-electric heating from +/// photo-electric heating by dust-grains +/// +/// For added context, we currently use equation 1 (and possibly eqn 2) of +/// [Wolfire+95](https://ui.adsabs.harvard.edu/abs/1995ApJ...443..152W/abstract) +/// to implement photo-electric heating double gammah_rate(double units, chemistry_data *my_chemistry); double gamma_isrf_rate(double units, chemistry_data *my_chemistry); double gamma_isrf2_rate(double units, chemistry_data *my_chemistry); From 0c5e830b77c6d04d780287b10f2991b1dcacabf0 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 9 Jan 2026 10:57:45 -0500 Subject: [PATCH 19/30] improve photoheating docstring --- src/clib/dust/gas_heat_cool.hpp | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/src/clib/dust/gas_heat_cool.hpp b/src/clib/dust/gas_heat_cool.hpp index f24c254a0..0b6e6bda7 100644 --- a/src/clib/dust/gas_heat_cool.hpp +++ b/src/clib/dust/gas_heat_cool.hpp @@ -25,9 +25,25 @@ /// @todo Come up with a better name? namespace grackle::impl::dust_gas_edot { -/// update edot, in place, with contributions from Photo-electric heating by +/// update edot, in place, with contributions from photo-electric heating by /// UV-irradiated dust /// +/// Based on eqn 1 of +/// [Wolfire+95](https://ui.adsabs.harvard.edu/abs/1995ApJ...443..152W/abstract): +/// @code{unparsed} +/// Γeff = ε * Γ +/// Γ = (1e-24 * G₀) erg cm⁻³ s⁻¹ +/// @endcode +/// where G₀ is the interstellar FUV radiation (in Habing units) +/// +/// The value of `photoelectric_heating` config parameter controls the +/// meaning of @p gammah: +/// - when the `photoelectric_heating` is 1, @p gammah directly holds Γ +/// - otherwise, it holds: (Γ/G₀) +/// +/// `photoelectric_heating` also controls modeling of ε +/// +/// @note /// Each of the 1D arrays is only valid for the specified @p idx_range /// /// @param [in,out] edot 1D array being used to accumulate the net rate of @@ -41,7 +57,8 @@ namespace grackle::impl::dust_gas_edot { /// @param[in] itmask Specifies the general iteration-mask of the @p idx_range /// for this calculation. /// @param[in] my_chemistry holds a number of configuration parameters. -/// @param[in] gammah Parameterizes the calculation +/// @param[in] gammah Parameterizes the calculation. Precomputed by +/// @ref gammah_rate. /// @param[in] idx_range Specifies the current index-range /// @param[in] dom_inv /// From 0ea284487c83a2cc9e92e4fa428db10959c86e58 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 15 Jan 2026 23:17:44 -0500 Subject: [PATCH 20/30] incremental commit --- src/clib/cool1d_multi_g.cpp | 157 +++++++++++++++++++----------------- 1 file changed, 85 insertions(+), 72 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index dc0af4e01..db61dfa5d 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -27,6 +27,88 @@ #include "internal_types.hpp" #include "utils-cpp.hpp" +void dust_cooling_rate( + double* edot, const double* tgas, const double* tdust, + const double* dust2gas, const double* rhoH, + const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + IndexRange idx_range, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::View d, const double* gasgr, + double* Ldst, + grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask_metal[i] != MASK_FALSE) { + if (my_chemistry->dust_species == 0) { + Ldst[i] = + -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * rhoH[i]; + + } else { + if (my_chemistry->use_multiple_dust_temperatures == 0) { + Ldst[i] = -gasgr[i] * (tgas[i] - tdust[i]) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } else { + if (my_chemistry->dust_species > 0) { + Ldst[i] = + -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::MgSiO3_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::AC_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } + + if (my_chemistry->dust_species > 1) { + Ldst[i] = + Ldst[i] - + (gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::SiM_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::FeM_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::Mg2SiO4_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::Fe3O4_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::SiO2_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::MgO_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::FeS_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::Al2O3_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } + + if (my_chemistry->dust_species > 2) { + Ldst[i] = + Ldst[i] - + (gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::ref_org_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::vol_org_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::H2O_ice_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } + } + } + + edot[i] = edot[i] + Ldst[i]; + } + } +} void grackle::impl::cool1d_multi_g( int imetal, int iter, double* edot, double* tgas, double* mmw, double* p2d, double* tdust, double* metallicity, double* dust2gas, double* rhoH, @@ -1160,78 +1242,9 @@ void grackle::impl::cool1d_multi_g( // Calculate dust cooling rate if (anydust != MASK_FALSE) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask_metal[i] != MASK_FALSE) { - if (my_chemistry->dust_species == 0) { - Ldst[i] = -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * - rhoH[i]; - - } else { - if (my_chemistry->use_multiple_dust_temperatures == 0) { - Ldst[i] = -gasgr[i] * (tgas[i] - tdust[i]) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } else { - if (my_chemistry->dust_species > 0) { - Ldst[i] = - -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::MgSiO3_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::AC_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } - - if (my_chemistry->dust_species > 1) { - Ldst[i] = - Ldst[i] - - (gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::SiM_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::FeM_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::Mg2SiO4_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::Fe3O4_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::SiO2_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::MgO_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::FeS_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::Al2O3_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } - - if (my_chemistry->dust_species > 2) { - Ldst[i] = - Ldst[i] - - (gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::ref_org_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::vol_org_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::H2O_ice_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } - } - } - - edot[i] = edot[i] + Ldst[i]; - } - } + dust_cooling_rate(edot, tgas, tdust, dust2gas, rhoH, itmask_metal, + my_chemistry, idx_range, grain_temperatures, d, + gasgr.data(), Ldst.data(), gas_grainsp_heatrate); } // Compute continuum opacity From 0e9d5bcc4df1d3ce4b56b593d4afcfec6da891bf Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 15 Jan 2026 23:23:31 -0500 Subject: [PATCH 21/30] rename dust_cooling_rate->dust_gas_edot::update_edot_dust_cooling_rate --- src/clib/cool1d_multi_g.cpp | 89 ++------------------------------- src/clib/dust/gas_heat_cool.hpp | 82 ++++++++++++++++++++++++++++++ 2 files changed, 86 insertions(+), 85 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index db61dfa5d..cc1939e8a 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -27,88 +27,6 @@ #include "internal_types.hpp" #include "utils-cpp.hpp" -void dust_cooling_rate( - double* edot, const double* tgas, const double* tdust, - const double* dust2gas, const double* rhoH, - const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, - IndexRange idx_range, - grackle::impl::GrainSpeciesCollection grain_temperatures, - grackle::impl::View d, const double* gasgr, - double* Ldst, - grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate) { - for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { - if (itmask_metal[i] != MASK_FALSE) { - if (my_chemistry->dust_species == 0) { - Ldst[i] = - -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * rhoH[i]; - - } else { - if (my_chemistry->use_multiple_dust_temperatures == 0) { - Ldst[i] = -gasgr[i] * (tgas[i] - tdust[i]) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } else { - if (my_chemistry->dust_species > 0) { - Ldst[i] = - -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::MgSiO3_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::AC_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } - - if (my_chemistry->dust_species > 1) { - Ldst[i] = - Ldst[i] - - (gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::SiM_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::FeM_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::Mg2SiO4_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::Fe3O4_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::SiO2_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::MgO_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::FeS_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::Al2O3_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } - - if (my_chemistry->dust_species > 2) { - Ldst[i] = - Ldst[i] - - (gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::ref_org_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::vol_org_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::H2O_ice_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } - } - } - - edot[i] = edot[i] + Ldst[i]; - } - } -} void grackle::impl::cool1d_multi_g( int imetal, int iter, double* edot, double* tgas, double* mmw, double* p2d, double* tdust, double* metallicity, double* dust2gas, double* rhoH, @@ -1242,9 +1160,10 @@ void grackle::impl::cool1d_multi_g( // Calculate dust cooling rate if (anydust != MASK_FALSE) { - dust_cooling_rate(edot, tgas, tdust, dust2gas, rhoH, itmask_metal, - my_chemistry, idx_range, grain_temperatures, d, - gasgr.data(), Ldst.data(), gas_grainsp_heatrate); + dust_gas_edot::update_edot_dust_cooling_rate( + edot, tgas, tdust, dust2gas, rhoH, itmask_metal, my_chemistry, + idx_range, grain_temperatures, d, gasgr.data(), Ldst.data(), + gas_grainsp_heatrate); } // Compute continuum opacity diff --git a/src/clib/dust/gas_heat_cool.hpp b/src/clib/dust/gas_heat_cool.hpp index 0b6e6bda7..872b1cb31 100644 --- a/src/clib/dust/gas_heat_cool.hpp +++ b/src/clib/dust/gas_heat_cool.hpp @@ -170,6 +170,88 @@ inline void update_edot_dust_recombination( } } +void update_edot_dust_cooling_rate( + double* edot, const double* tgas, const double* tdust, + const double* dust2gas, const double* rhoH, + const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + IndexRange idx_range, + grackle::impl::GrainSpeciesCollection grain_temperatures, + grackle::impl::View d, const double* gasgr, double* Ldst, + grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + if (itmask_metal[i] != MASK_FALSE) { + if (my_chemistry->dust_species == 0) { + Ldst[i] = + -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * rhoH[i]; + + } else { + if (my_chemistry->use_multiple_dust_temperatures == 0) { + Ldst[i] = -gasgr[i] * (tgas[i] - tdust[i]) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } else { + if (my_chemistry->dust_species > 0) { + Ldst[i] = + -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::MgSiO3_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::AC_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } + + if (my_chemistry->dust_species > 1) { + Ldst[i] = + Ldst[i] - + (gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::SiM_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeM_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::FeM_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Mg2SiO4_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::Mg2SiO4_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Fe3O4_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::Fe3O4_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiO2_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::SiO2_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgO_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::MgO_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::FeS_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::FeS_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::Al2O3_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::Al2O3_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } + + if (my_chemistry->dust_species > 2) { + Ldst[i] = + Ldst[i] - + (gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::ref_org_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::vol_org_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::vol_org_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::H2O_ice_dust][i] * + (tgas[i] - grain_temperatures + .data[OnlyGrainSpLUT::H2O_ice_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; + } + } + } + + edot[i] = edot[i] + Ldst[i]; + } + } +} + } // namespace grackle::impl::dust_gas_edot #endif // DUST_GAS_HEAT_COOL_HPP \ No newline at end of file From 800aa29a534d690970a44aea3e9aee8b35162bee Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 15 Jan 2026 23:28:58 -0500 Subject: [PATCH 22/30] get rid of the Ldst vector --- src/clib/cool1d_multi_g.cpp | 4 +--- src/clib/dust/gas_heat_cool.hpp | 21 +++++++++++---------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index cc1939e8a..a5fb904b9 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -176,7 +176,6 @@ void grackle::impl::cool1d_multi_g( std::vector LCO(my_fields->grid_dimension[0]); std::vector LOH(my_fields->grid_dimension[0]); std::vector LH2O(my_fields->grid_dimension[0]); - std::vector Ldst(my_fields->grid_dimension[0]); std::vector alpha(my_fields->grid_dimension[0]); std::vector alphad(my_fields->grid_dimension[0]); std::vector lshield_con(my_fields->grid_dimension[0]); @@ -1162,8 +1161,7 @@ void grackle::impl::cool1d_multi_g( if (anydust != MASK_FALSE) { dust_gas_edot::update_edot_dust_cooling_rate( edot, tgas, tdust, dust2gas, rhoH, itmask_metal, my_chemistry, - idx_range, grain_temperatures, d, gasgr.data(), Ldst.data(), - gas_grainsp_heatrate); + idx_range, grain_temperatures, d, gasgr.data(), gas_grainsp_heatrate); } // Compute continuum opacity diff --git a/src/clib/dust/gas_heat_cool.hpp b/src/clib/dust/gas_heat_cool.hpp index 872b1cb31..cf7971f28 100644 --- a/src/clib/dust/gas_heat_cool.hpp +++ b/src/clib/dust/gas_heat_cool.hpp @@ -176,21 +176,22 @@ void update_edot_dust_cooling_rate( const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, IndexRange idx_range, grackle::impl::GrainSpeciesCollection grain_temperatures, - grackle::impl::View d, const double* gasgr, double* Ldst, + grackle::impl::View d, const double* gasgr, grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask_metal[i] != MASK_FALSE) { + double Ldst; if (my_chemistry->dust_species == 0) { - Ldst[i] = + Ldst = -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * rhoH[i]; } else { if (my_chemistry->use_multiple_dust_temperatures == 0) { - Ldst[i] = -gasgr[i] * (tgas[i] - tdust[i]) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; + Ldst = -gasgr[i] * (tgas[i] - tdust[i]) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; } else { if (my_chemistry->dust_species > 0) { - Ldst[i] = + Ldst = -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * (tgas[i] - grain_temperatures .data[OnlyGrainSpLUT::MgSiO3_dust][i]) + @@ -201,8 +202,8 @@ void update_edot_dust_cooling_rate( } if (my_chemistry->dust_species > 1) { - Ldst[i] = - Ldst[i] - + Ldst = + Ldst - (gas_grainsp_heatrate.data[OnlyGrainSpLUT::SiM_dust][i] * (tgas[i] - grain_temperatures.data[OnlyGrainSpLUT::SiM_dust][i]) + @@ -231,8 +232,8 @@ void update_edot_dust_cooling_rate( } if (my_chemistry->dust_species > 2) { - Ldst[i] = - Ldst[i] - + Ldst = + Ldst - (gas_grainsp_heatrate.data[OnlyGrainSpLUT::ref_org_dust][i] * (tgas[i] - grain_temperatures .data[OnlyGrainSpLUT::ref_org_dust][i]) + @@ -247,7 +248,7 @@ void update_edot_dust_cooling_rate( } } - edot[i] = edot[i] + Ldst[i]; + edot[i] = edot[i] + Ldst; } } } From c176db508ca2ae4bbc5f460da9c7835f60d1687c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 15:08:52 -0500 Subject: [PATCH 23/30] address compiler/linter warnings by making it clear that Ldst is always initialized --- src/clib/dust/gas_heat_cool.hpp | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/clib/dust/gas_heat_cool.hpp b/src/clib/dust/gas_heat_cool.hpp index cf7971f28..3de5734a2 100644 --- a/src/clib/dust/gas_heat_cool.hpp +++ b/src/clib/dust/gas_heat_cool.hpp @@ -16,6 +16,7 @@ #include "grackle.h" #include "index_helper.h" #include "internal_types.hpp" +#include "status_reporting.h" #include @@ -173,7 +174,7 @@ inline void update_edot_dust_recombination( void update_edot_dust_cooling_rate( double* edot, const double* tgas, const double* tdust, const double* dust2gas, const double* rhoH, - const gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + const gr_mask_type* itmask_metal, const chemistry_data* my_chemistry, IndexRange idx_range, grackle::impl::GrainSpeciesCollection grain_temperatures, grackle::impl::View d, const double* gasgr, @@ -185,21 +186,19 @@ void update_edot_dust_cooling_rate( Ldst = -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * rhoH[i]; - } else { + } else { // my_chemistry->dust_species > 0 if (my_chemistry->use_multiple_dust_temperatures == 0) { Ldst = -gasgr[i] * (tgas[i] - tdust[i]) * d(i, idx_range.j, idx_range.k) * rhoH[i]; } else { - if (my_chemistry->dust_species > 0) { - Ldst = - -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * - (tgas[i] - grain_temperatures - .data[OnlyGrainSpLUT::MgSiO3_dust][i]) + - gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * - (tgas[i] - - grain_temperatures.data[OnlyGrainSpLUT::AC_dust][i])) * - d(i, idx_range.j, idx_range.k) * rhoH[i]; - } + Ldst = + -(gas_grainsp_heatrate.data[OnlyGrainSpLUT::MgSiO3_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::MgSiO3_dust][i]) + + gas_grainsp_heatrate.data[OnlyGrainSpLUT::AC_dust][i] * + (tgas[i] - + grain_temperatures.data[OnlyGrainSpLUT::AC_dust][i])) * + d(i, idx_range.j, idx_range.k) * rhoH[i]; if (my_chemistry->dust_species > 1) { Ldst = From abb1aef298dac9938d918aa7bcce74b6952b4a8b Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 15:46:57 -0500 Subject: [PATCH 24/30] improve docstrings pertaining to update_edot_dust_cooling_rate --- src/clib/cool1d_multi_g.cpp | 4 +-- src/clib/dust/gas_heat_cool.hpp | 48 ++++++++++++++++++++++++---- src/clib/rate_functions.c | 2 -- src/include/grackle_rate_functions.h | 7 ++++ 4 files changed, 51 insertions(+), 10 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index a5fb904b9..7a79d288e 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1160,8 +1160,8 @@ void grackle::impl::cool1d_multi_g( // Calculate dust cooling rate if (anydust != MASK_FALSE) { dust_gas_edot::update_edot_dust_cooling_rate( - edot, tgas, tdust, dust2gas, rhoH, itmask_metal, my_chemistry, - idx_range, grain_temperatures, d, gasgr.data(), gas_grainsp_heatrate); + edot, tgas, tdust, grain_temperatures, dust2gas, rhoH, itmask_metal, + my_chemistry, idx_range, d, gasgr.data(), gas_grainsp_heatrate); } // Compute continuum opacity diff --git a/src/clib/dust/gas_heat_cool.hpp b/src/clib/dust/gas_heat_cool.hpp index 3de5734a2..9db6d6c18 100644 --- a/src/clib/dust/gas_heat_cool.hpp +++ b/src/clib/dust/gas_heat_cool.hpp @@ -171,14 +171,50 @@ inline void update_edot_dust_recombination( } } +/// update edot, in place, to account for energy transfer between gas and dust +/// grains. +/// +/// Each of the 1D arrays is only valid for the specified @p idx_range. +/// +/// The rates of energy transfer are computed, from interpolation tables, at +/// the same time that @p tdust (or @p grain_temperatures) is computed. Based +/// on the configuration, the original interpolation table was computed by +/// @ref gasGrain_rate or @ref gasGrain2_rate. +/// +/// @todo +/// Consider refactoring so that the updates to @p edot are directly computed +/// at the same time as the dust temperatures. If we opt not to do this, then +/// we should probably refactor to ensure that the definition of the +/// gas-to-grain rate is more consistent between configurations (and so that we +/// don't need to pass d to this function). +/// +/// @param [in,out] edot 1D array being used to accumulate the net rate of +/// change in internal energy of the gas +/// @param[in] tgas 1D array of gas temperatures +/// @param[in] tdust, grain_temperatures Specifies dust temperatures. The +/// former is used in some configurations to provide a 1D array of dust +/// temperatures for all dust. The latter is used in other configurations +/// to provide separate 1d dust temperature arrays for each modeled grain +/// species. +/// @param[in] dust2gas 1D array of ratios between dust & gas densities +/// @param[in] rhoH 1D array holding the Hydrogen mass density +/// @param[in] itmask_metal Specifies the iteration-mask of the @p idx_range +/// for this calculation. +/// @param[in] my_chemistry Holds various configuration parameters +/// @param[in] idx_range Specifies the current index-range +/// @param[in] d The 3D mass density array +/// @param[in] gasgr, gas_grainsp_heatrate Specifies gas-to-grain heat transfer +/// rates (the precise definition of this "rate" varies a lot!). The +/// former is used in some configurations to provide a 1D array of rates +/// for all dust. The latter is used in other configurations to provide +/// separate 1d dust temperature arrays for each modeled grain species. void update_edot_dust_cooling_rate( double* edot, const double* tgas, const double* tdust, - const double* dust2gas, const double* rhoH, - const gr_mask_type* itmask_metal, const chemistry_data* my_chemistry, - IndexRange idx_range, - grackle::impl::GrainSpeciesCollection grain_temperatures, - grackle::impl::View d, const double* gasgr, - grackle::impl::GrainSpeciesCollection gas_grainsp_heatrate) { + const GrainSpeciesCollection& grain_temperatures, const double* dust2gas, + const double* rhoH, const gr_mask_type* itmask_metal, + const chemistry_data* my_chemistry, IndexRange idx_range, + View d, const double* gasgr, + const GrainSpeciesCollection& gas_grainsp_heatrate) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { if (itmask_metal[i] != MASK_FALSE) { double Ldst; diff --git a/src/clib/rate_functions.c b/src/clib/rate_functions.c index 1b58caa72..2fc656fdd 100644 --- a/src/clib/rate_functions.c +++ b/src/clib/rate_functions.c @@ -1375,7 +1375,6 @@ double cieco_rate(double T, double units, chemistry_data *my_chemistry) return cierate * (mh/2.0) / units; } -//Calculation of gas_grain. double gasGrain_rate(double T, double units, chemistry_data *my_chemistry) { double grain_coef; @@ -1409,7 +1408,6 @@ double gasGrain_rate(double T, double units, chemistry_data *my_chemistry) } } -//Calculation of gas_grain2. double gasGrain2_rate(double T, double units, chemistry_data *my_chemistry) { //Variables. diff --git a/src/include/grackle_rate_functions.h b/src/include/grackle_rate_functions.h index dfe17e4f8..ca07bd5f8 100644 --- a/src/include/grackle_rate_functions.h +++ b/src/include/grackle_rate_functions.h @@ -99,7 +99,14 @@ double HDlow_rate(double T, double units, chemistry_data *my_chemistry); double cie_thin_cooling_rate(double T); double cieco_rate(double T, double units, chemistry_data *my_chemistry); +/// Calculate gas_grain, the Gas/grain energy transfer rate double gasGrain_rate(double T, double units, chemistry_data *my_chemistry); + +/// Calculate gas_grain2 +/// +/// The resulting value is similar to the value returned by @ref gasGrain_rate. +/// It is used to compute the Gas/grain energy transfer rate for arbitrary size +/// distributions. double gasGrain2_rate(double T, double units, chemistry_data *my_chemistry); /// Calculate regr (as in GRain REcombination cooling) From cb3a4104217c0a6af789ff31ec15ddaf906f1bd0 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 17:21:43 -0500 Subject: [PATCH 25/30] Factor out dust_related_props --- src/clib/CMakeLists.txt | 1 + src/clib/cool1d_multi_g.cpp | 64 +++------------------- src/clib/dust/misc.hpp | 103 ++++++++++++++++++++++++++++++++++++ 3 files changed, 111 insertions(+), 57 deletions(-) create mode 100644 src/clib/dust/misc.hpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 976d52c87..0074437c6 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -112,6 +112,7 @@ add_library(Grackle_Grackle dust/gas_heat_cool.hpp dust/grain_species_info.cpp dust/grain_species_info.hpp dust/lookup_dust_rates1d.hpp + dust/misc.hpp init_misc_species_cool_rates.cpp init_misc_species_cool_rates.hpp initialize_chemistry_data.cpp initialize_dust_yields.cpp initialize_dust_yields.hpp diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 7a79d288e..9dc2911e2 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -13,6 +13,7 @@ // This file was initially generated automatically during conversion of the // cool1d_multi_g function from FORTRAN to C++ +#include "dust/misc.hpp" #include "dust/gas_heat_cool.hpp" #include @@ -1099,63 +1100,12 @@ void grackle::impl::cool1d_multi_g( } } - // Compute grain size increment - if ((my_chemistry->use_dust_density_field > 0) && - (my_chemistry->dust_species > 0)) { - grackle::impl::fortran_wrapper::calc_grain_size_increment_1d( - dom, idx_range, itmask_metal, my_chemistry, my_rates, my_fields, - internal_dust_prop_buf); - } - - // Calculate dust to gas ratio AND interstellar radiation field - // -> an earlier version of this logic would store values @ indices - // where `itmask_metal(i) .ne. MASK_FALSE` - // -> this was undesirable, b/c these quantities are required for - // photo-electric heating, which can occur when - // `itmask_metal(i) .eq. MASK_FALSE` (we can revisit this choice - // later). Moreover, in most cases, these calculations will be - // faster when there is no branching - - if ((anydust != MASK_FALSE) || (my_chemistry->photoelectric_heating > 0)) { - if (my_chemistry->use_dust_density_field > 0) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - // REMINDER: use of `itmask` over `itmask_metal` is - // currently required by Photo-electric heating - if (itmask[i] != MASK_FALSE) { - // it may be faster to remove this branching - dust2gas[i] = dust(i, idx_range.j, idx_range.k) / - d(i, idx_range.j, idx_range.k); - } - } - } else { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - dust2gas[i] = my_chemistry->local_dust_to_gas_ratio * metallicity[i]; - } - } - } - - if ((anydust != MASK_FALSE) || (my_chemistry->photoelectric_heating > 1)) { - if (my_chemistry->use_isrf_field > 0) { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - myisrf[i] = isrf_habing(i, idx_range.j, idx_range.k); - } - } else { - for (i = idx_range.i_start; i <= idx_range.i_end; i++) { - myisrf[i] = my_chemistry->interstellar_radiation_field; - } - } - } - - // compute dust temperature and cooling due to dust - if (anydust != MASK_FALSE) { - // TODO: trad -> comp2 - grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( - comp2, tgas, tdust, metallicity, dust2gas, cool1dmulti_buf.mynh, - cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr.data(), - myisrf.data(), kappa_tot.data(), my_chemistry, my_rates, my_fields, - idx_range, grain_temperatures, gas_grainsp_heatrate, grain_kappa, - logTlininterp_buf, internal_dust_prop_buf); - } + dust_related_props(tgas, metallicity, itmask, itmask_metal, my_chemistry, + my_rates, my_fields, idx_range, logTlininterp_buf, + cool1dmulti_buf, d, dust, isrf_habing, dom, coolunit, + comp2, dust2gas, myisrf, tdust, grain_temperatures, gasgr, + gas_grainsp_heatrate, internal_dust_prop_buf, grain_kappa, + kappa_tot, anydust); // Calculate dust cooling rate if (anydust != MASK_FALSE) { diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp new file mode 100644 index 000000000..55cd2ae5d --- /dev/null +++ b/src/clib/dust/misc.hpp @@ -0,0 +1,103 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Miscellaneous functionality +/// +//===----------------------------------------------------------------------===// +#ifndef GRACKLE_DUST_MISC_HPP +#define GRACKLE_DUST_MISC_HPP + +#include +#include "grackle.h" +#include "fortran_func_decls.h" // gr_mask_type +#include "fortran_func_wrappers.hpp" +#include "index_helper.h" // IndexHelper +#include "utils-cpp.hpp" // View +#include "internal_types.hpp" +#include "dust_props.hpp" + +namespace grackle::impl { + +/// Calculate various dust related properties +/// +/// @param[in] tgas gas temperature +inline void dust_related_props( + double* tgas, double* metallicity, const gr_mask_type* itmask, + gr_mask_type* itmask_metal, chemistry_data* my_chemistry, + chemistry_data_storage* my_rates, grackle_field_data* my_fields, + IndexRange idx_range, LogTLinInterpScratchBuf logTlininterp_buf, + Cool1DMultiScratchBuf cool1dmulti_buf, View d, + View dust, View isrf_habing, double dom, + double coolunit, double comp2, double* dust2gas, + std::vector& myisrf, double* tdust, + GrainSpeciesCollection grain_temperatures, std::vector& gasgr, + GrainSpeciesCollection gas_grainsp_heatrate, + InternalDustPropBuf internal_dust_prop_buf, + GrainSpeciesCollection grain_kappa, std::vector& kappa_tot, + gr_mask_type anydust) { + // Compute grain size increment + if ((my_chemistry->use_dust_density_field > 0) && + (my_chemistry->dust_species > 0)) { + grackle::impl::fortran_wrapper::calc_grain_size_increment_1d( + dom, idx_range, itmask_metal, my_chemistry, my_rates, my_fields, + internal_dust_prop_buf); + } + + // Calculate dust to gas ratio AND interstellar radiation field + // -> an earlier version of this logic would store values @ indices + // where `itmask_metal(i) .ne. MASK_FALSE` + // -> this was undesirable, b/c these quantities are required for + // photo-electric heating, which can occur when + // `itmask_metal(i) .eq. MASK_FALSE` (we can revisit this choice + // later). Moreover, in most cases, these calculations will be + // faster when there is no branching + + if ((anydust != MASK_FALSE) || (my_chemistry->photoelectric_heating > 0)) { + if (my_chemistry->use_dust_density_field > 0) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + // REMINDER: use of `itmask` over `itmask_metal` is + // currently required by Photo-electric heating + if (itmask[i] != MASK_FALSE) { + // it may be faster to remove this branching + dust2gas[i] = dust(i, idx_range.j, idx_range.k) / + d(i, idx_range.j, idx_range.k); + } + } + } else { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + dust2gas[i] = my_chemistry->local_dust_to_gas_ratio * metallicity[i]; + } + } + } + + if ((anydust != MASK_FALSE) || (my_chemistry->photoelectric_heating > 1)) { + if (my_chemistry->use_isrf_field > 0) { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + myisrf[i] = isrf_habing(i, idx_range.j, idx_range.k); + } + } else { + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { + myisrf[i] = my_chemistry->interstellar_radiation_field; + } + } + } + + // compute dust temperature and cooling due to dust + if (anydust != MASK_FALSE) { + // TODO: trad -> comp2 + grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( + comp2, tgas, tdust, metallicity, dust2gas, cool1dmulti_buf.mynh, + cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr.data(), + myisrf.data(), kappa_tot.data(), my_chemistry, my_rates, my_fields, + idx_range, grain_temperatures, gas_grainsp_heatrate, grain_kappa, + logTlininterp_buf, internal_dust_prop_buf); + } +} +} // namespace grackle::impl + +#endif // GRACKLE_DUST_MISC_HPP From a69782f490ae322f061afb39fe1d0b688a138596 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 17:39:10 -0500 Subject: [PATCH 26/30] adjust a few more arguments --- src/clib/cool1d_multi_g.cpp | 17 +++++------------ src/clib/dust/misc.hpp | 30 ++++++++++++++++++++++-------- 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 9dc2911e2..51ec2c8d8 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -78,9 +78,6 @@ void grackle::impl::cool1d_multi_g( grackle::impl::View metal( my_fields->metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); - grackle::impl::View dust( - my_fields->dust_density, my_fields->grid_dimension[0], - my_fields->grid_dimension[1], my_fields->grid_dimension[2]); grackle::impl::View Vheat( my_fields->volumetric_heating_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); @@ -93,9 +90,6 @@ void grackle::impl::cool1d_multi_g( grackle::impl::View photogamma( my_fields->RT_heating_rate, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); - grackle::impl::View isrf_habing( - my_fields->isrf_habing, my_fields->grid_dimension[0], - my_fields->grid_dimension[1], my_fields->grid_dimension[2]); grackle::impl::View CI( my_fields->CI_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); @@ -1100,12 +1094,11 @@ void grackle::impl::cool1d_multi_g( } } - dust_related_props(tgas, metallicity, itmask, itmask_metal, my_chemistry, - my_rates, my_fields, idx_range, logTlininterp_buf, - cool1dmulti_buf, d, dust, isrf_habing, dom, coolunit, - comp2, dust2gas, myisrf, tdust, grain_temperatures, gasgr, - gas_grainsp_heatrate, internal_dust_prop_buf, grain_kappa, - kappa_tot, anydust); + dust_related_props( + anydust, tgas, metallicity, itmask, itmask_metal, my_chemistry, my_rates, + my_fields, internalu, idx_range, logTlininterp_buf, cool1dmulti_buf, + comp2, dust2gas, myisrf, tdust, grain_temperatures, gasgr, + gas_grainsp_heatrate, grain_kappa, kappa_tot, internal_dust_prop_buf); // Calculate dust cooling rate if (anydust != MASK_FALSE) { diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp index 55cd2ae5d..51b2cc6fc 100644 --- a/src/clib/dust/misc.hpp +++ b/src/clib/dust/misc.hpp @@ -25,21 +25,24 @@ namespace grackle::impl { /// Calculate various dust related properties /// +/// @param[in] anydust Whether dust chemistry is enabled /// @param[in] tgas gas temperature inline void dust_related_props( - double* tgas, double* metallicity, const gr_mask_type* itmask, - gr_mask_type* itmask_metal, chemistry_data* my_chemistry, - chemistry_data_storage* my_rates, grackle_field_data* my_fields, + gr_mask_type anydust, double* tgas, double* metallicity, + const gr_mask_type* itmask, gr_mask_type* itmask_metal, + chemistry_data* my_chemistry, chemistry_data_storage* my_rates, + grackle_field_data* my_fields, InternalGrUnits internalu, IndexRange idx_range, LogTLinInterpScratchBuf logTlininterp_buf, - Cool1DMultiScratchBuf cool1dmulti_buf, View d, - View dust, View isrf_habing, double dom, - double coolunit, double comp2, double* dust2gas, + Cool1DMultiScratchBuf cool1dmulti_buf, double comp2, double* dust2gas, std::vector& myisrf, double* tdust, GrainSpeciesCollection grain_temperatures, std::vector& gasgr, GrainSpeciesCollection gas_grainsp_heatrate, - InternalDustPropBuf internal_dust_prop_buf, GrainSpeciesCollection grain_kappa, std::vector& kappa_tot, - gr_mask_type anydust) { + InternalDustPropBuf internal_dust_prop_buf) { + // get relevant unit values + double dom = internalu_calc_dom_(internalu); + double coolunit = internalu.coolunit; + // Compute grain size increment if ((my_chemistry->use_dust_density_field > 0) && (my_chemistry->dust_species > 0)) { @@ -58,6 +61,13 @@ inline void dust_related_props( // faster when there is no branching if ((anydust != MASK_FALSE) || (my_chemistry->photoelectric_heating > 0)) { + grackle::impl::View d( + my_fields->density, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust( + my_fields->dust_density, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + if (my_chemistry->use_dust_density_field > 0) { for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { // REMINDER: use of `itmask` over `itmask_metal` is @@ -77,6 +87,10 @@ inline void dust_related_props( if ((anydust != MASK_FALSE) || (my_chemistry->photoelectric_heating > 1)) { if (my_chemistry->use_isrf_field > 0) { + grackle::impl::View isrf_habing( + my_fields->isrf_habing, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int i = idx_range.i_start; i <= idx_range.i_end; i++) { myisrf[i] = isrf_habing(i, idx_range.j, idx_range.k); } From 3343a3d546fa68a61d0369cd8a6303e73475cb7d Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 17:53:02 -0500 Subject: [PATCH 27/30] adjust the final few arguments --- src/clib/cool1d_multi_g.cpp | 11 ++++++----- src/clib/dust/misc.hpp | 14 ++++++-------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 51ec2c8d8..43227ce4c 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1094,11 +1094,12 @@ void grackle::impl::cool1d_multi_g( } } - dust_related_props( - anydust, tgas, metallicity, itmask, itmask_metal, my_chemistry, my_rates, - my_fields, internalu, idx_range, logTlininterp_buf, cool1dmulti_buf, - comp2, dust2gas, myisrf, tdust, grain_temperatures, gasgr, - gas_grainsp_heatrate, grain_kappa, kappa_tot, internal_dust_prop_buf); + dust_related_props(anydust, tgas, metallicity, itmask, itmask_metal, + my_chemistry, my_rates, my_fields, internalu, idx_range, + logTlininterp_buf, cool1dmulti_buf, comp2, dust2gas, + myisrf.data(), tdust, grain_temperatures, gasgr.data(), + gas_grainsp_heatrate, grain_kappa, kappa_tot.data(), + internal_dust_prop_buf); // Calculate dust cooling rate if (anydust != MASK_FALSE) { diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp index 51b2cc6fc..77ff4d1ef 100644 --- a/src/clib/dust/misc.hpp +++ b/src/clib/dust/misc.hpp @@ -12,7 +12,6 @@ #ifndef GRACKLE_DUST_MISC_HPP #define GRACKLE_DUST_MISC_HPP -#include #include "grackle.h" #include "fortran_func_decls.h" // gr_mask_type #include "fortran_func_wrappers.hpp" @@ -34,10 +33,9 @@ inline void dust_related_props( grackle_field_data* my_fields, InternalGrUnits internalu, IndexRange idx_range, LogTLinInterpScratchBuf logTlininterp_buf, Cool1DMultiScratchBuf cool1dmulti_buf, double comp2, double* dust2gas, - std::vector& myisrf, double* tdust, - GrainSpeciesCollection grain_temperatures, std::vector& gasgr, - GrainSpeciesCollection gas_grainsp_heatrate, - GrainSpeciesCollection grain_kappa, std::vector& kappa_tot, + double* myisrf, double* tdust, GrainSpeciesCollection grain_temperatures, + double* gasgr, GrainSpeciesCollection gas_grainsp_heatrate, + GrainSpeciesCollection grain_kappa, double* kappa_tot, InternalDustPropBuf internal_dust_prop_buf) { // get relevant unit values double dom = internalu_calc_dom_(internalu); @@ -106,9 +104,9 @@ inline void dust_related_props( // TODO: trad -> comp2 grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( comp2, tgas, tdust, metallicity, dust2gas, cool1dmulti_buf.mynh, - cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr.data(), - myisrf.data(), kappa_tot.data(), my_chemistry, my_rates, my_fields, - idx_range, grain_temperatures, gas_grainsp_heatrate, grain_kappa, + cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr, myisrf, + kappa_tot, my_chemistry, my_rates, my_fields, idx_range, + grain_temperatures, gas_grainsp_heatrate, grain_kappa, logTlininterp_buf, internal_dust_prop_buf); } } From 37732297ba3dbe422cd81b48c5fb4394dd931798 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 18:25:18 -0500 Subject: [PATCH 28/30] add a docstring for dust_related_props --- src/clib/dust/misc.hpp | 40 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp index 77ff4d1ef..2587e6a6e 100644 --- a/src/clib/dust/misc.hpp +++ b/src/clib/dust/misc.hpp @@ -24,8 +24,46 @@ namespace grackle::impl { /// Calculate various dust related properties /// +/// @todo +/// - Stop passing @p cool1dmulti_buf and instead pass the `mynh` and +/// `gasgr_tdust` members directly. +/// - The former is unmodified while the latter is used as a scratch buffer +/// - to minimize conflicts wait for the `calc_all_tdust_gasgr_1d_g` +/// transcription PR +/// - Make all non-modified `double*` arguments `const double` (this must wait +/// until after the transcription PRs for the wrapped fortran routines) +/// /// @param[in] anydust Whether dust chemistry is enabled -/// @param[in] tgas gas temperature +/// @param[in] tgas 1d array of gas temperature +/// @param[in] itmask Specifies the general iteration-mask of the @p idx_range +/// for this calculation. +/// @param[in] itmask_metal Specifies the metal/dust-specific iteration-mask of +/// the @p idx_range for this calculation. +/// @param[in] my_chemistry holds a number of configuration parameters. +/// @param[in] my_rates Holds assorted rate data and other internal +/// configuration info. +/// @param[in] my_fields Specifies the field data. +/// @param[in] internalu Specifies Grackle's internal unit-system +/// @param[in] idx_range Specifies the current index-range +/// @param[in] logTlininterp_buf hold values for each location in @p idx_range +/// that are used to linearly interpolate tables with respect to the +/// natural log of @p tgas. +/// @param[in,out] cool1dmulti_buf The mynh member holds a 1D array of number +/// densities while the gasgr_tdust member is just a scratch buffer +/// @param[in] comp2 Holds the CMB temperature at the current redshift +/// @param[out] dust2gas Holds the computed dust-to-gas ratio at each +/// location in the index range. In other words, this holds the dust mass +/// per unit gas mass (only used in certain configuration) +/// @param[in,out] myisrf a scratch buffer that may be used to temporarily +/// record the interstellar radiation field +/// @param[out] tdust, grain_temperatures dust temperatures may be written +/// to one of these variables, based on configuration +/// @param[out] gasgr, gas_grainsp_heatrate Grain/gas energy transfer rates may +/// be written to one of these variables, based on configuration +/// @param[out] grain_kappa, kappa_tot Opacity-related information may be +/// written to one of these variables, based on configuration +/// @param[in,out] internal_dust_prop_buf Holds scratch-space for holding +/// grain-specific information inline void dust_related_props( gr_mask_type anydust, double* tgas, double* metallicity, const gr_mask_type* itmask, gr_mask_type* itmask_metal, From e537cfff3e85da474a195796b219176474fc534c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 16 Jan 2026 18:34:05 -0500 Subject: [PATCH 29/30] Shuffle a few last things --- src/clib/cool1d_multi_g.cpp | 8 ++++---- src/clib/dust/misc.hpp | 33 ++++++++++++++++----------------- 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 43227ce4c..404eca659 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1096,10 +1096,10 @@ void grackle::impl::cool1d_multi_g( dust_related_props(anydust, tgas, metallicity, itmask, itmask_metal, my_chemistry, my_rates, my_fields, internalu, idx_range, - logTlininterp_buf, cool1dmulti_buf, comp2, dust2gas, - myisrf.data(), tdust, grain_temperatures, gasgr.data(), - gas_grainsp_heatrate, grain_kappa, kappa_tot.data(), - internal_dust_prop_buf); + logTlininterp_buf, comp2, dust2gas, tdust, + grain_temperatures, gasgr.data(), gas_grainsp_heatrate, + kappa_tot.data(), grain_kappa, cool1dmulti_buf, + myisrf.data(), internal_dust_prop_buf); // Calculate dust cooling rate if (anydust != MASK_FALSE) { diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp index 2587e6a6e..86175567c 100644 --- a/src/clib/dust/misc.hpp +++ b/src/clib/dust/misc.hpp @@ -25,13 +25,13 @@ namespace grackle::impl { /// Calculate various dust related properties /// /// @todo +/// All the following items should wait until after the transcription PRs for +/// the wrapped Fortran routines are merged: /// - Stop passing @p cool1dmulti_buf and instead pass the `mynh` and -/// `gasgr_tdust` members directly. -/// - The former is unmodified while the latter is used as a scratch buffer -/// - to minimize conflicts wait for the `calc_all_tdust_gasgr_1d_g` -/// transcription PR -/// - Make all non-modified `double*` arguments `const double` (this must wait -/// until after the transcription PRs for the wrapped fortran routines) +/// `gasgr_tdust` members directly. The former is unmodified while the latter +/// is used as a scratch buffer +/// - Make all non-modified `double*` arguments `const double*` +/// - rename the @p comp2 argument so that it's known as Tcmb /// /// @param[in] anydust Whether dust chemistry is enabled /// @param[in] tgas 1d array of gas temperature @@ -48,20 +48,20 @@ namespace grackle::impl { /// @param[in] logTlininterp_buf hold values for each location in @p idx_range /// that are used to linearly interpolate tables with respect to the /// natural log of @p tgas. -/// @param[in,out] cool1dmulti_buf The mynh member holds a 1D array of number -/// densities while the gasgr_tdust member is just a scratch buffer /// @param[in] comp2 Holds the CMB temperature at the current redshift /// @param[out] dust2gas Holds the computed dust-to-gas ratio at each /// location in the index range. In other words, this holds the dust mass /// per unit gas mass (only used in certain configuration) -/// @param[in,out] myisrf a scratch buffer that may be used to temporarily -/// record the interstellar radiation field /// @param[out] tdust, grain_temperatures dust temperatures may be written /// to one of these variables, based on configuration /// @param[out] gasgr, gas_grainsp_heatrate Grain/gas energy transfer rates may /// be written to one of these variables, based on configuration -/// @param[out] grain_kappa, kappa_tot Opacity-related information may be +/// @param[out] kappa_tot, grain_kappa Opacity-related information may be /// written to one of these variables, based on configuration +/// @param[in,out] cool1dmulti_buf The mynh member holds a 1D array of number +/// densities while the gasgr_tdust member is just a scratch buffer +/// @param[in,out] myisrf a scratch buffer that may be used to temporarily +/// record the interstellar radiation field /// @param[in,out] internal_dust_prop_buf Holds scratch-space for holding /// grain-specific information inline void dust_related_props( @@ -70,11 +70,11 @@ inline void dust_related_props( chemistry_data* my_chemistry, chemistry_data_storage* my_rates, grackle_field_data* my_fields, InternalGrUnits internalu, IndexRange idx_range, LogTLinInterpScratchBuf logTlininterp_buf, - Cool1DMultiScratchBuf cool1dmulti_buf, double comp2, double* dust2gas, - double* myisrf, double* tdust, GrainSpeciesCollection grain_temperatures, - double* gasgr, GrainSpeciesCollection gas_grainsp_heatrate, - GrainSpeciesCollection grain_kappa, double* kappa_tot, - InternalDustPropBuf internal_dust_prop_buf) { + double comp2, double* dust2gas, double* tdust, + GrainSpeciesCollection grain_temperatures, double* gasgr, + GrainSpeciesCollection gas_grainsp_heatrate, double* kappa_tot, + GrainSpeciesCollection grain_kappa, Cool1DMultiScratchBuf cool1dmulti_buf, + double* myisrf, InternalDustPropBuf internal_dust_prop_buf) { // get relevant unit values double dom = internalu_calc_dom_(internalu); double coolunit = internalu.coolunit; @@ -139,7 +139,6 @@ inline void dust_related_props( // compute dust temperature and cooling due to dust if (anydust != MASK_FALSE) { - // TODO: trad -> comp2 grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( comp2, tgas, tdust, metallicity, dust2gas, cool1dmulti_buf.mynh, cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr, myisrf, From 1daecfd9f32573900dbb1e774df5e5ee9ee79b8f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sat, 17 Jan 2026 11:23:27 -0500 Subject: [PATCH 30/30] After a littlen thought, I wanted to make a few remaining changes to the arguments --- src/clib/cool1d_multi_g.cpp | 8 ++++---- src/clib/dust/misc.hpp | 32 ++++++++++++++++---------------- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 404eca659..a4d7583d6 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -1094,11 +1094,11 @@ void grackle::impl::cool1d_multi_g( } } - dust_related_props(anydust, tgas, metallicity, itmask, itmask_metal, - my_chemistry, my_rates, my_fields, internalu, idx_range, - logTlininterp_buf, comp2, dust2gas, tdust, + dust_related_props(anydust, tgas, cool1dmulti_buf.mynh, metallicity, itmask, + itmask_metal, my_chemistry, my_rates, my_fields, internalu, + idx_range, logTlininterp_buf, comp2, dust2gas, tdust, grain_temperatures, gasgr.data(), gas_grainsp_heatrate, - kappa_tot.data(), grain_kappa, cool1dmulti_buf, + kappa_tot.data(), grain_kappa, cool1dmulti_buf.gasgr_tdust, myisrf.data(), internal_dust_prop_buf); // Calculate dust cooling rate diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp index 86175567c..bffb8752a 100644 --- a/src/clib/dust/misc.hpp +++ b/src/clib/dust/misc.hpp @@ -27,14 +27,15 @@ namespace grackle::impl { /// @todo /// All the following items should wait until after the transcription PRs for /// the wrapped Fortran routines are merged: -/// - Stop passing @p cool1dmulti_buf and instead pass the `mynh` and -/// `gasgr_tdust` members directly. The former is unmodified while the latter -/// is used as a scratch buffer -/// - Make all non-modified `double*` arguments `const double*` -/// - rename the @p comp2 argument so that it's known as Tcmb +/// - Make all non-modified `double*` arguments `const double*`. +/// - consider accepting `rhoH` instead of @p nH and computing the value of `nH` +/// just when we need it (to be clear, there's a chance that this plan may +/// not be entirely viable). /// /// @param[in] anydust Whether dust chemistry is enabled /// @param[in] tgas 1d array of gas temperature +/// @param[in] nH 1d array of Hydrogen number densities +/// @param[in] metallicity 1d array of metallicities /// @param[in] itmask Specifies the general iteration-mask of the @p idx_range /// for this calculation. /// @param[in] itmask_metal Specifies the metal/dust-specific iteration-mask of @@ -48,7 +49,7 @@ namespace grackle::impl { /// @param[in] logTlininterp_buf hold values for each location in @p idx_range /// that are used to linearly interpolate tables with respect to the /// natural log of @p tgas. -/// @param[in] comp2 Holds the CMB temperature at the current redshift +/// @param[in] trad Holds the CMB temperature at the current redshift /// @param[out] dust2gas Holds the computed dust-to-gas ratio at each /// location in the index range. In other words, this holds the dust mass /// per unit gas mass (only used in certain configuration) @@ -58,23 +59,23 @@ namespace grackle::impl { /// be written to one of these variables, based on configuration /// @param[out] kappa_tot, grain_kappa Opacity-related information may be /// written to one of these variables, based on configuration -/// @param[in,out] cool1dmulti_buf The mynh member holds a 1D array of number -/// densities while the gasgr_tdust member is just a scratch buffer +/// @param[in,out] gasgr_tdust A 1D array of that acts as a scratch buffer +/// (with some refactoring, this can probably be removed) /// @param[in,out] myisrf a scratch buffer that may be used to temporarily /// record the interstellar radiation field /// @param[in,out] internal_dust_prop_buf Holds scratch-space for holding /// grain-specific information inline void dust_related_props( - gr_mask_type anydust, double* tgas, double* metallicity, + gr_mask_type anydust, double* tgas, double* nH, double* metallicity, const gr_mask_type* itmask, gr_mask_type* itmask_metal, chemistry_data* my_chemistry, chemistry_data_storage* my_rates, grackle_field_data* my_fields, InternalGrUnits internalu, IndexRange idx_range, LogTLinInterpScratchBuf logTlininterp_buf, - double comp2, double* dust2gas, double* tdust, + double trad, double* dust2gas, double* tdust, GrainSpeciesCollection grain_temperatures, double* gasgr, GrainSpeciesCollection gas_grainsp_heatrate, double* kappa_tot, - GrainSpeciesCollection grain_kappa, Cool1DMultiScratchBuf cool1dmulti_buf, - double* myisrf, InternalDustPropBuf internal_dust_prop_buf) { + GrainSpeciesCollection grain_kappa, double* gasgr_tdust, double* myisrf, + InternalDustPropBuf internal_dust_prop_buf) { // get relevant unit values double dom = internalu_calc_dom_(internalu); double coolunit = internalu.coolunit; @@ -140,10 +141,9 @@ inline void dust_related_props( // compute dust temperature and cooling due to dust if (anydust != MASK_FALSE) { grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( - comp2, tgas, tdust, metallicity, dust2gas, cool1dmulti_buf.mynh, - cool1dmulti_buf.gasgr_tdust, itmask_metal, coolunit, gasgr, myisrf, - kappa_tot, my_chemistry, my_rates, my_fields, idx_range, - grain_temperatures, gas_grainsp_heatrate, grain_kappa, + trad, tgas, tdust, metallicity, dust2gas, nH, gasgr_tdust, itmask_metal, + coolunit, gasgr, myisrf, kappa_tot, my_chemistry, my_rates, my_fields, + idx_range, grain_temperatures, gas_grainsp_heatrate, grain_kappa, logTlininterp_buf, internal_dust_prop_buf); } }