diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 7036b28e8..0074437c6 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -109,8 +109,10 @@ 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_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 02793cc90..a4d7583d6 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -13,6 +13,9 @@ // 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 #include #include @@ -75,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]); @@ -90,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]); @@ -120,8 +117,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, ih2cox, + min_metallicity; double comp1, comp2; // Performing heap allocations for all of the subsequent buffers within this @@ -174,7 +171,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]); @@ -1098,138 +1094,18 @@ 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(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.gasgr_tdust, + myisrf.data(), internal_dust_prop_buf); // 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_gas_edot::update_edot_dust_cooling_rate( + edot, tgas, tdust, grain_temperatures, dust2gas, rhoH, itmask_metal, + my_chemistry, idx_range, d, gasgr.data(), gas_grainsp_heatrate); } // Compute continuum opacity @@ -1521,80 +1397,17 @@ 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; - } - } - } + 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)) { - 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_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); } // Compton cooling or heating and X-ray compton heating @@ -1969,4 +1782,4 @@ void grackle::impl::cool1d_multi_g( grackle::impl::drop_GrainSpeciesCollection(&gas_grainsp_heatrate); return; -} \ No newline at end of file +} diff --git a/src/clib/dust/gas_heat_cool.hpp b/src/clib/dust/gas_heat_cool.hpp new file mode 100644 index 000000000..9db6d6c18 --- /dev/null +++ b/src/clib/dust/gas_heat_cool.hpp @@ -0,0 +1,293 @@ +//===----------------------------------------------------------------------===// +// +// 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_HEAT_COOL_HPP +#define DUST_GAS_HEAT_COOL_HPP + +#include "fortran_func_decls.h" // gr_mask_type +#include "grackle.h" +#include "index_helper.h" +#include "internal_types.hpp" +#include "status_reporting.h" + +#include + +/// 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 +/// +/// 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 +/// 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] 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. Precomputed by +/// @ref gammah_rate. +/// @param[in] idx_range Specifies the current index-range +/// @param[in] dom_inv +/// +/// @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 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 + // 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; + }; + + if (my_chemistry->photoelectric_heating == 1) { + 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); + } + } + + // 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_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]; + update_edot_fn(i, gammaha_eff); + } + } + + // 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_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))) + + ((3.7e-2 * std::pow((tgas[i] / 1.e4), 0.7)) / + (1. + (pe_X / 5000.))); + double gammaha_eff = gammah * pe_eps * isrf[i]; + update_edot_fn(i, gammaha_eff); + } + } + } +} + +/// 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 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 +/// @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 +/// 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 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 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 +/// +/// @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 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, double dom_inv) { + 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]] - + regr[logTlininterp_buf.indixe[i] - 1]); + + double grbeta = 0.74 / std::pow(tgas[i], 0.068); + edot[i] = + edot[i] - + cur_regr * std::pow((isrf[i] * dom_inv / e_density[i]), grbeta) * + e_density[i] * rhoH[i] * dust2gas[i] / local_dust_to_gas_ratio; + } + } +} + +/// 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 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; + if (my_chemistry->dust_species == 0) { + Ldst = + -gasgr[i] * (tgas[i] - tdust[i]) * dust2gas[i] * rhoH[i] * rhoH[i]; + + } 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 { + 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 = + Ldst - + (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 = + Ldst - + (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; + } + } +} + +} // namespace grackle::impl::dust_gas_edot + +#endif // DUST_GAS_HEAT_COOL_HPP \ No newline at end of file diff --git a/src/clib/dust/misc.hpp b/src/clib/dust/misc.hpp new file mode 100644 index 000000000..bffb8752a --- /dev/null +++ b/src/clib/dust/misc.hpp @@ -0,0 +1,152 @@ +//===----------------------------------------------------------------------===// +// +// 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 "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 +/// +/// @todo +/// All the following items should wait until after the transcription PRs for +/// the wrapped Fortran routines are merged: +/// - 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 +/// 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] 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) +/// @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] kappa_tot, grain_kappa Opacity-related information may be +/// written to one of these variables, based on configuration +/// @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* 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 trad, double* dust2gas, double* tdust, + GrainSpeciesCollection grain_temperatures, double* gasgr, + GrainSpeciesCollection gas_grainsp_heatrate, double* kappa_tot, + 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; + + // 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)) { + 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 + // 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) { + 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); + } + } 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) { + grackle::impl::fortran_wrapper::calc_all_tdust_gasgr_1d_g( + 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); + } +} +} // namespace grackle::impl + +#endif // GRACKLE_DUST_MISC_HPP diff --git a/src/clib/internal_types.hpp b/src/clib/internal_types.hpp index 7aaa30fa7..618ab8346 100644 --- a/src/clib/internal_types.hpp +++ b/src/clib/internal_types.hpp @@ -125,9 +125,7 @@ struct Cool1DMultiScratchBuf { double* tgasold = nullptr; double* mynh = nullptr; double* myde = nullptr; - double* gammaha_eff = nullptr; double* gasgr_tdust = nullptr; - double* regr = nullptr; }; /// used to help implement the visitor design pattern @@ -143,9 +141,7 @@ 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); } diff --git a/src/clib/rate_functions.c b/src/clib/rate_functions.c index 154b82d31..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. @@ -1421,7 +1419,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) @@ -1437,11 +1434,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 9d3647481..ca07bd5f8 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 @@ -88,12 +99,32 @@ 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) +/// +/// 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); + +/// 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);