diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 39e8d0a60..3da0abf95 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -106,6 +106,7 @@ add_library(Grackle_Grackle dust/calc_kappa_grain.cpp dust/calc_kappa_grain.hpp dust/calc_tdust_1d_g.cpp dust/calc_tdust_1d_g.hpp dust/calc_tdust_3d.cpp dust/calc_tdust_3d.hpp + dust/dust_growth_and_destruction.cpp dust/dust_growth_and_destruction.hpp dust/gas_heat_cool.hpp dust/grain_species_info.cpp dust/grain_species_info.hpp dust/lookup_dust_rates1d.hpp diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 926770c60..5a3c9eb6c 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -31,6 +31,7 @@ OBJS_CONFIG_LIB = \ dust/multi_grain_species/calc_grain_size_increment_species_1d.lo \ dust/calc_kappa_grain.lo \ dust/calc_tdust_1d_g.lo \ + dust/dust_growth_and_destruction.lo \ dust/grain_species_info.lo \ dynamic_api.lo \ api/units.lo \ @@ -49,7 +50,6 @@ OBJS_CONFIG_LIB = \ solve_rate_cool.lo \ support/status_reporting.lo \ update_UVbackground_rates.lo \ - calc_tdust_1d_g.lo \ dust/calc_tdust_3d.lo \ calc_grain_size_increment_1d.lo \ rate_functions.lo \ diff --git a/src/clib/ceiling_species.hpp b/src/clib/ceiling_species.hpp index 8cc78a79b..c0ae42eff 100644 --- a/src/clib/ceiling_species.hpp +++ b/src/clib/ceiling_species.hpp @@ -69,7 +69,7 @@ inline void ceiling_species(int imetal, chemistry_data* my_chemistry, GRIMPL_NS::View metal( my_fields->metal_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); - GRIMPL_NS::View dust( + grackle::impl::View dust( my_fields->dust_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); GRIMPL_NS::View DM( diff --git a/src/clib/cool1d_multi_g.cpp b/src/clib/cool1d_multi_g.cpp index 06e7e7975..41eccd334 100644 --- a/src/clib/cool1d_multi_g.cpp +++ b/src/clib/cool1d_multi_g.cpp @@ -883,7 +883,6 @@ void grackle::impl::cool1d_multi_g( itmask_metal[i] = MASK_FALSE; } } - 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, diff --git a/src/clib/dust/calc_tdust_3d.cpp b/src/clib/dust/calc_tdust_3d.cpp index 7cd50c376..6c821b705 100644 --- a/src/clib/dust/calc_tdust_3d.cpp +++ b/src/clib/dust/calc_tdust_3d.cpp @@ -154,17 +154,24 @@ void calc_tdust_3d( } } - // Compute grain size increment - - if ( (my_chemistry->use_dust_density_field > 0) && (my_chemistry->dust_species > 0) ) { + // // Set itmask to false for dust-poor cells + // if (my_chemistry->use_dust_density_field > 0) { + // for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + // if (dust(i,j,k) < 1.e-9 * d(i,j,k)) { + // itmask_metal[i] = MASK_FALSE; + // } + // } + // } - calc_grain_size_increment_1d ( + // Compute grain size increment + if ( (my_chemistry->use_dust_density_field > 0) && (my_chemistry->dust_species > 0) + && (my_chemistry->dust_model == 0) ) { + grackle::impl::calc_grain_size_increment_1d ( dom, idx_range, itmask_metal.data(), my_chemistry, my_rates->opaque_storage->grain_species_info, my_rates->opaque_storage->inject_pathway_props, my_fields, internal_dust_prop_buf ); - } for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { @@ -202,7 +209,7 @@ void calc_tdust_3d( // endif if (my_chemistry->use_dust_density_field > 0) { - dust2gas[i] = dust(i,j,k) / d(i,j,k); + dust2gas[i] = dust(i,j,k) / (d(i,j,k)); } else { dust2gas[i] = my_chemistry->local_dust_to_gas_ratio * metallicity[i]; } diff --git a/src/clib/dust/dust_growth_and_destruction.cpp b/src/clib/dust/dust_growth_and_destruction.cpp new file mode 100644 index 000000000..adafea51c --- /dev/null +++ b/src/clib/dust/dust_growth_and_destruction.cpp @@ -0,0 +1,912 @@ +#include +#include +#include +#include "dust_growth_and_destruction.hpp" +#include "internal_types.hpp" +#include "utils-cpp.hpp" + +namespace { +const double sec_per_year = 3.155e7; +const double sec_per_Myr_local = 1.0e6 * sec_per_year; + +const double tiny_value = 1.0e-20; +const double huge_value = 1.0e+20; + +// Reference temperature for the bulk (dust_model=1, non-species) accretion +// timescale. +const double bulk_growth_t_ref = 20.0; + +// COLIBRE accretion / sputtering reference values +// [REF: Trayford+2026 MNRAS 545, staf2040, section 3.3]. +const double colibre_growth_t_ref = 10.0; // K +const double colibre_growth_nH_ref = 10.0; // cm^-3 +const double colibre_sputter_t_ref = 2.0e6; // K + +const double atomic_C = 12.01; +const double atomic_O = 16.00; +const double atomic_Mg = 24.31; +const double atomic_Si = 28.09; +const double atomic_Fe = 55.85; + +double limited_expm1(double x) { + if (!std::isfinite(x) || x <= 0.0) return 0.0; + return std::expm1(std::min(x, 50.0)); +} + +// Equivalent per-time rates whose product with dt reproduces one exponential +// e-folding update: rate * dt = rho_dust * (exp(+-dt/tau) - 1). +double e_fold_growth_rate(double rho_dust, double dt, double tau) { + if (rho_dust <= 0.0 || dt <= 0.0 || tau <= 0.0) return 0.0; + return rho_dust * limited_expm1(dt / tau) / dt; +} + +double e_fold_loss_rate(double rho_dust, double dt, double inv_tau) { + if (rho_dust <= 0.0 || dt <= 0.0 || inv_tau <= 0.0) return 0.0; + double x = std::min(dt * inv_tau, 50.0); + return rho_dust * (-std::expm1(-x)) / dt; +} + +// Solar metal mass fractions for the tracked dust-forming elements, relative +// to total solar metals. These match gracklepy.utilities.convenience, which +// seeds the dust_species_track gas reservoirs from metal_density. +const double solar_frac_C = 0.15925782394660776; +const double solar_frac_O = 0.4242932765702842; +const double solar_frac_Mg = 0.045644817372018066; +const double solar_frac_Si = 0.052744600629574714; +const double solar_frac_Fe = 0.08523143041944482; + +// Gate thresholds for dust evolution: skip cells where dust or metals are a +// negligible fraction of the baryon density. The 1e-9 ratio matches the +// metal-poor threshold used in make_consistent; below it, dust evolution +// would just integrate numerical noise. +const double dust_gate_threshold = 1.0e-9; +const double metal_gate_threshold = 1.0e-9; + +// Solar-to-local abundance ratio entering the COLIBRE accretion timescale. +// Degenerate inputs yield an effectively infinite timescale (no growth). +double abundance_ratio(double solar_eps, double local_eps) { + if (solar_eps <= 0.0 || local_eps <= 0.0) return huge_value; + return solar_eps / local_eps; +} + +// Hydrogen-nuclei mass density census. Counts H from whichever species fields +// are evolved so that the n_H entering the dust rates stays consistent with +// the active network (primordial_chemistry / metal_chemistry); falls back to +// HydrogenFractionByMass * (rho_gas - rho_metal) when no primordial species +// are evolved. Returns 0 for pathological cells (metals >= total density). +struct HNucleiCensus { + bool use_H, use_H2, use_HD, use_HeH, use_metal_H; + double hfrac; + grackle::impl::View d, metal; + grackle::impl::View HI, HII, HM, H2I, H2II, HDI, HDII, HeHII; + grackle::impl::View OH, H2O, CH, CH2, OHII, H2OII, H3OII; + + HNucleiCensus(const chemistry_data* my_chemistry, + const grackle_field_data* my_fields) + : use_H(my_chemistry->primordial_chemistry > 0), + use_H2(my_chemistry->primordial_chemistry > 1), + use_HD(my_chemistry->primordial_chemistry > 2), + use_HeH(my_chemistry->primordial_chemistry > 3), + use_metal_H(my_chemistry->metal_chemistry == 1), + hfrac(my_chemistry->HydrogenFractionByMass), + d(view_(my_fields, my_fields->density, true)), + metal(view_(my_fields, my_fields->metal_density, true)), + HI(view_(my_fields, my_fields->HI_density, use_H)), + HII(view_(my_fields, my_fields->HII_density, use_H)), + HM(view_(my_fields, my_fields->HM_density, use_H2)), + H2I(view_(my_fields, my_fields->H2I_density, use_H2)), + H2II(view_(my_fields, my_fields->H2II_density, use_H2)), + HDI(view_(my_fields, my_fields->HDI_density, use_HD)), + HDII(view_(my_fields, my_fields->HDII_density, use_HeH)), + HeHII(view_(my_fields, my_fields->HeHII_density, use_HeH)), + OH(view_(my_fields, my_fields->OH_density, use_metal_H)), + H2O(view_(my_fields, my_fields->H2O_density, use_metal_H)), + CH(view_(my_fields, my_fields->CH_density, use_metal_H)), + CH2(view_(my_fields, my_fields->CH2_density, use_metal_H)), + OHII(view_(my_fields, my_fields->OHII_density, use_metal_H)), + H2OII(view_(my_fields, my_fields->H2OII_density, use_metal_H)), + H3OII(view_(my_fields, my_fields->H3OII_density, use_metal_H)) {} + + double rho_H(int i, int j, int k) const { + double rho_nonmetal = d(i, j, k) - metal(i, j, k); + if (rho_nonmetal <= 0.0) return 0.0; + if (!use_H) return hfrac * rho_nonmetal; + + double rho = HI(i, j, k) + HII(i, j, k); + if (use_H2) { + rho += HM(i, j, k) + H2I(i, j, k) + H2II(i, j, k); + } + if (use_HD) { + rho += HDI(i, j, k) / 3.0; + } + if (use_HeH) { + rho += HDII(i, j, k) / 3.0 + HeHII(i, j, k) / 5.0; + } + if (use_metal_H) { + rho += OH(i, j, k) / 17.0 + H2O(i, j, k) * (2.0 / 18.0) + + CH(i, j, k) / 13.0 + CH2(i, j, k) * (2.0 / 14.0) + + OHII(i, j, k) / 17.0 + H2OII(i, j, k) * (2.0 / 18.0) + + H3OII(i, j, k) * (3.0 / 19.0); + } + return rho; + } + +private: + static grackle::impl::View view_( + const grackle_field_data* my_fields, gr_float* ptr, bool active) { + return grackle::impl::View( + active ? ptr : my_fields->density, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + } +}; + +} // namespace + +// ========================================== +// DUST GROWTH (ACCRETION) +// ========================================== +void grackle::impl::dust_growth(chemistry_data* my_chemistry, + grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, + const gr_mask_type* itmask, + const double* dt_value, const double* t_gas, + double* growth_dM) { + 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]); + grackle::impl::View metal( + my_fields->metal_density, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + + double dens_proper = internalu.urho * std::pow(internalu.a_value, 3); + double tau_ref = + my_chemistry->dust_growth_tauref * 1e9 * sec_per_year / internalu.tbase1; + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + growth_dM[i] = 0.0; + + if (itmask[i] != MASK_FALSE) { + double rho_dust = dust(i, idx_range.j, idx_range.k); + // metal_density holds total gas-phase metals (C, O, and all others) + double rho_metal = metal(i, idx_range.j, idx_range.k); + double rho_gas = d(i, idx_range.j, idx_range.k); + + // No metals to accrete onto grains — skip growth + if (rho_metal < metal_gate_threshold * rho_gas) { + continue; + } + + double temp = t_gas[i]; + double dt = dt_value[i]; + + // tau_accr = tau_ref * (densref / rho_metal,phys) * sqrt(T_ref / T) + // * Z_sun: the density and metallicity scalings combine into + // the single rho_metal term since rho_metal = Z * rho_gas. + double tau_accr0 = tau_ref * + (my_chemistry->dust_growth_densref / dens_proper) * + std::pow(bulk_growth_t_ref / temp, 0.5); + double rho_metal_eff = std::max(rho_metal, tiny_value); + double tau_accr = + tau_accr0 * (my_chemistry->SolarMetalFractionByMass / rho_metal_eff); + tau_accr = std::clamp(tau_accr, tiny_value, huge_value); + + // Fraction of the combined dust+metal reservoir still in the gas + // phase; throttles growth as metals deplete. + double frac_metal_available = + std::clamp(rho_metal / (rho_dust + rho_metal), 0.0, 1.0); + double growth_rate = frac_metal_available * (rho_dust / tau_accr); + growth_dM[i] = std::min(growth_rate, rho_metal / dt); + } + } +} + +// ========================================== +// DUST GROWTH (SPECIES-SPECIFIC: Mg-silicate + Fe-silicate + carbonaceous) +// ========================================== +// COLIBRE-style accretion [REF: Trayford+2026 MNRAS 545, staf2040, section +// 3.3]. Carbonaceous growth is bottlenecked by gas-phase C. Silicate growth +// is computed once for the total silicate reservoir using the maximum +// depletion factor over O, Si, and Mg+Fe, then split between Mg2SiO4 and +// Fe2SiO4 endmembers according to the local gas-phase Mg/Fe number abundance. +// The dense-gas rate uses a clumped hydrogen density n_H' = C n_H; thermal +// sputtering in dust_destruction_species deliberately uses the unclumped n_H. +// Growth is applied as an exponential e-folding update, returned as an +// equivalent per-time rate for dust_update_species(). +void grackle::impl::dust_growth_species( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, const double* t_gas, double* growth_dM_mg_silicate, + double* growth_dM_fe_silicate, double* growth_dM_carbon) { + 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_mg_sil( + my_fields->dust_density_mg_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_fe_sil( + my_fields->dust_density_fe_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_carb( + my_fields->dust_density_carbonaceous, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mC( + my_fields->metal_density_carbon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mO( + my_fields->metal_density_oxygen, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mMg( + my_fields->metal_density_magnesium, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mSi( + my_fields->metal_density_silicon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mFe( + my_fields->metal_density_iron, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + HNucleiCensus h_census(my_chemistry, my_fields); + + double dens_proper = internalu.urho * std::pow(internalu.a_value, 3); + // tau_ref values are stored in Myr (COLIBRE Table 2). + double tau_ref_sil = my_chemistry->dust_growth_tauref_silicate * + sec_per_Myr_local / internalu.tbase1; + double tau_ref_carb = my_chemistry->dust_growth_tauref_carbon * + sec_per_Myr_local / internalu.tbase1; + // COLIBRE reference values S_ref = 0.3, a_ref = 0.1 micron; non-positive + // parameters disable growth via an effectively infinite timescale. + double sticking_factor = (my_chemistry->dust_growth_sticking_coeff > 0.0) + ? 0.3 / my_chemistry->dust_growth_sticking_coeff + : huge_value; + double grain_size_factor = (my_chemistry->dust_grainsize > 0.0) + ? my_chemistry->dust_grainsize / 0.1 + : huge_value; + + // Solar abundances (per H mass, then per H nucleus) of the tracked + // elements, used as the depletion reference in the tau scalings. + double solar_nonmetal = + std::max(1.0 - my_chemistry->SolarMetalFractionByMass, tiny_value); + double solar_H = std::max( + my_chemistry->HydrogenFractionByMass * solar_nonmetal, tiny_value); + double solar_mass_C = + my_chemistry->SolarMetalFractionByMass * solar_frac_C / solar_H; + double solar_mass_O = + my_chemistry->SolarMetalFractionByMass * solar_frac_O / solar_H; + double solar_mass_Mg = + my_chemistry->SolarMetalFractionByMass * solar_frac_Mg / solar_H; + double solar_mass_Si = + my_chemistry->SolarMetalFractionByMass * solar_frac_Si / solar_H; + double solar_mass_Fe = + my_chemistry->SolarMetalFractionByMass * solar_frac_Fe / solar_H; + double solar_eps_C = solar_mass_C / atomic_C; + double solar_eps_O = solar_mass_O / atomic_O; + double solar_eps_Mg = solar_mass_Mg / atomic_Mg; + double solar_eps_Si = solar_mass_Si / atomic_Si; + double solar_eps_Fe = solar_mass_Fe / atomic_Fe; + double solar_eps_MgFe = solar_eps_Mg + solar_eps_Fe; + + // Sub-grid clumping boost: C rises log-linearly in n_H from 1 at nH_min to + // factor_max at nH_max. + auto clumping_factor = [&](double nH) -> double { + double cmax = std::max(my_chemistry->dust_growth_clumping_factor_max, 1.0); + double nmin = + std::max(my_chemistry->dust_growth_clumping_nH_min, tiny_value); + double nmax = std::max(my_chemistry->dust_growth_clumping_nH_max, + nmin * (1.0 + 1.0e-12)); + if (nH <= nmin || cmax <= 1.0) return 1.0; + if (nH >= nmax) return cmax; + double x = (std::log10(nH) - std::log10(nmin)) / + (std::log10(nmax) - std::log10(nmin)); + return std::pow(cmax, std::clamp(x, 0.0, 1.0)); + }; + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + growth_dM_mg_silicate[i] = 0.0; + growth_dM_fe_silicate[i] = 0.0; + growth_dM_carbon[i] = 0.0; + + if (itmask[i] != MASK_FALSE) { + double rho_gas = d(i, idx_range.j, idx_range.k); + double rho_dust_mg_sil_i = dust_mg_sil(i, idx_range.j, idx_range.k); + double rho_dust_fe_sil_i = dust_fe_sil(i, idx_range.j, idx_range.k); + double rho_dust_carb_i = dust_carb(i, idx_range.j, idx_range.k); + + double rho_C = mC(i, idx_range.j, idx_range.k); + double rho_O = mO(i, idx_range.j, idx_range.k); + double rho_Mg = mMg(i, idx_range.j, idx_range.k); + double rho_Si = mSi(i, idx_range.j, idx_range.k); + double rho_Fe = mFe(i, idx_range.j, idx_range.k); + + double T = std::max(t_gas[i], tiny_value); + double dt = dt_value[i]; + + double rho_H_nuclei = h_census.rho_H(i, idx_range.j, idx_range.k); + if (rho_gas <= 0.0 || rho_H_nuclei <= 0.0) { + continue; + } + double nH = rho_H_nuclei * dens_proper / mh; + if (nH <= 0.0) { + continue; + } + + double nH_eff = nH * clumping_factor(nH); + double accr_struct = (colibre_growth_nH_ref / nH_eff) * + std::pow(colibre_growth_t_ref / T, 0.5) * + grain_size_factor * sticking_factor; + double inv_rho_H = 1.0 / rho_H_nuclei; + double eps_C = rho_C * inv_rho_H / atomic_C; + double eps_O = rho_O * inv_rho_H / atomic_O; + double eps_Mg = rho_Mg * inv_rho_H / atomic_Mg; + double eps_Si = rho_Si * inv_rho_H / atomic_Si; + double eps_Fe = rho_Fe * inv_rho_H / atomic_Fe; + double eps_MgFe = eps_Mg + eps_Fe; + + // ---------- Carbonaceous: rate-limited by gas-phase C ---------- + if (rho_C >= metal_gate_threshold * rho_gas && rho_dust_carb_i > 0.0 && + dt > 0.0) { + double tau_accr_carb = + tau_ref_carb * accr_struct * abundance_ratio(solar_eps_C, eps_C); + tau_accr_carb = std::clamp(tau_accr_carb, tiny_value, huge_value); + double rate = e_fold_growth_rate(rho_dust_carb_i, dt, tau_accr_carb); + growth_dM_carbon[i] = std::min(rate, rho_C / dt); + } + + // ---------- Silicate: COLIBRE Mg+Fe composite bottleneck ---------- + double rho_dust_sil_i = rho_dust_mg_sil_i + rho_dust_fe_sil_i; + if (rho_dust_sil_i > 0.0 && dt > 0.0 && eps_O > 0.0 && eps_Si > 0.0 && + eps_MgFe > 0.0) { + double sil_ratio = + std::max({abundance_ratio(solar_eps_O, eps_O), + abundance_ratio(solar_eps_Si, eps_Si), + abundance_ratio(solar_eps_MgFe, eps_MgFe)}); + double tau_accr = tau_ref_sil * accr_struct * sil_ratio; + tau_accr = std::clamp(tau_accr, tiny_value, huge_value); + double sil_rate = e_fold_growth_rate(rho_dust_sil_i, dt, tau_accr); + + // Split total silicate accretion by endmember molecule abundance. + // Mass fractions include the different Mg2SiO4 / Fe2SiO4 molecule + // weights, so equal Mg and Fe number abundance reproduces the + // COLIBRE equal-molecule seed split. + double mg_weight = + eps_Mg * (2.0 * atomic_Mg + atomic_Si + 4.0 * atomic_O); + double fe_weight = + eps_Fe * (2.0 * atomic_Fe + atomic_Si + 4.0 * atomic_O); + double endmember_weight = mg_weight + fe_weight; + if (endmember_weight > 0.0) { + double mg_frac = std::clamp(mg_weight / endmember_weight, 0.0, 1.0); + growth_dM_mg_silicate[i] = sil_rate * mg_frac; + growth_dM_fe_silicate[i] = sil_rate * (1.0 - mg_frac); + } + } + } + } +} + +// ========================================== +// DUST DESTRUCTION (SNe + SPUTTERING) +// ========================================== +void grackle::impl::dust_destruction( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, double dt_full, const double* t_gas, + double* destruction_dM) { + 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]); + bool use_sne = (my_chemistry->use_sne_field > 0); + grackle::impl::View sne( + use_sne ? my_fields->sne_rate : my_fields->density, + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + bool use_tau_dest = (my_chemistry->use_tau_dest_field > 0); + grackle::impl::View tau_dest_field( + use_tau_dest ? my_fields->tau_dest : my_fields->density, + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + + double dens_proper = internalu.urho * std::pow(internalu.a_value, 3); + + // Gas mass shocked to >= 100 km/s per SN, M_s(100) = 6800 Msun + // (v_s/100 km/s)^-2 [REF: McKee 1989; used as in Li, Narayanan & Dave + // 2019], converted to code density * code volume. + double Ms100 = 6800.0 * my_chemistry->sne_coeff * + (100.0 / my_chemistry->sne_shockspeed) * + (100.0 / my_chemistry->sne_shockspeed) * SolarMass / + (internalu.urho * std::pow(internalu.uxyz, 3)); + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + destruction_dM[i] = 0.0; + + if (itmask[i] != MASK_FALSE) { + double rho_gas = d(i, idx_range.j, idx_range.k); + double rho_dust = dust(i, idx_range.j, idx_range.k); + + // No dust to destroy — skip destruction + if (rho_dust < dust_gate_threshold * rho_gas) { + continue; + } + + double sne_this = use_sne ? sne(i, idx_range.j, idx_range.k) : 0.0; + double temp = t_gas[i]; + double dt = dt_value[i]; + double dM_shock = 0.0; + + if (use_tau_dest) { + // user-provided destruction timescale + double tau_dest = tau_dest_field(i, idx_range.j, idx_range.k); + if (tau_dest <= 0) { + tau_dest = huge_value; + } + dM_shock = std::min(rho_dust / tau_dest, rho_dust / dt); + } else if (use_sne && sne_this > 0) { + // SN shock destruction: sne_rate holds the number of SNe per code + // volume delivered over the external timestep dt_full, so + // Ms100 * sne_this / dt_full is the shocked-gas mass rate and + // tau_dest the time to shock all gas in the cell. + double tau_dest = + rho_gas / (Ms100 * sne_this * my_chemistry->dust_destruction_eff) * + dt_full; + dM_shock = std::min(rho_dust / tau_dest, rho_dust / dt); + } + + // Thermal sputtering in hot gas [REF: Tsai & Mathews 1995 ApJ 448, + // 84]: tau_sp = 0.17 Gyr (a/0.1 um) (rho/1e-27 g cm^-3)^-1 + // [(2e6 K / T)^2.5 + 1]. The factor 3 converts the grain-radius + // e-folding time to a mass-loss time (m ~ a^3). + if (temp >= 1.0e5 && dM_shock < rho_dust / dt) { + double tau_sput = 1.7e8 * sec_per_year / internalu.tbase1 * + (my_chemistry->dust_grainsize / 0.1) * + (1.0e-27 / (dens_proper * rho_gas)) * + (std::pow((2.0e6 / temp), 2.5) + 1.0); + dM_shock = + std::min(dM_shock + 3.0 * rho_dust / tau_sput, rho_dust / dt); + } + + if (std::isnan(dM_shock)) { + std::cout << "dust_destruction: dM calculated as NaN" << std::endl; + dM_shock = 0.0; + } + + destruction_dM[i] = -dM_shock; + } + } +} + +// ========================================== +// DUST DESTRUCTION (SPECIES-SPECIFIC: Mg-silicate + Fe-silicate + carbonaceous) +// ========================================== +// SN-shock destruction + thermal sputtering applied independently to each +// dust species, gated by dust_species_track == 1. Carbonaceous (graphite) is +// the shock-vulnerability baseline; both silicate endmembers use the same +// silicate shock coefficient. Slavin+2015 Table 2 gives gas-cleared masses of +// 990 Msun for silicates and 600 Msun for carbonaceous grains in their +// standard SNR model, so silicates are destroyed about 1.65x faster +// [REF: Slavin, Dwek, Jones 2015 ApJ 803, 7; Jones+1996 ApJ 469, 740]. +// Thermal sputtering follows the species-independent COLIBRE/Tsai-Mathews +// timescale tau_sp = tau_ref (a/0.1) (n_H/cm^-3)^-1 [1 + (T/2e6 K)^-2.5]. +// Destruction is applied as exponential decay, returned as an equivalent +// per-time rate for dust_update_species(). +void grackle::impl::dust_destruction_species( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, double dt_full, const double* t_gas, + double* destruction_dM_mg_silicate, double* destruction_dM_fe_silicate, + double* destruction_dM_carbon) { + 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_mg_sil( + my_fields->dust_density_mg_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_fe_sil( + my_fields->dust_density_fe_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_carb( + my_fields->dust_density_carbonaceous, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + bool use_sne = (my_chemistry->use_sne_field > 0); + grackle::impl::View sne( + use_sne ? my_fields->sne_rate : my_fields->density, + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + bool use_tau_dest = (my_chemistry->use_tau_dest_field > 0); + grackle::impl::View tau_dest_field( + use_tau_dest ? my_fields->tau_dest : my_fields->density, + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + HNucleiCensus h_census(my_chemistry, my_fields); + + double dens_proper = internalu.urho * std::pow(internalu.a_value, 3); + + // Gas mass shocked per SN, M_s(100) = 6800 Msun (v_s/100 km/s)^-2 + // [REF: McKee 1989; used as in Li, Narayanan & Dave 2019]. + double Ms100 = 6800.0 * my_chemistry->sne_coeff * + (100.0 / my_chemistry->sne_shockspeed) * + (100.0 / my_chemistry->sne_shockspeed) * SolarMass / + (internalu.urho * std::pow(internalu.uxyz, 3)); + + // Species-specific shock-vulnerability multipliers. Graphite is the + // baseline (1.0); silicate follows the Slavin+2015 SNR gas-cleared mass + // ratio, 990/600 = 1.65. + const double shock_factor_carbon = 1.0; + const double shock_factor_silicate = 1.65; + + // COLIBRE/Tsai-Mathews thermal sputtering tau_ref (stored in Myr). + double tau_sput_ref = + my_chemistry->dust_sputter_tauref * sec_per_Myr_local / internalu.tbase1; + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + destruction_dM_mg_silicate[i] = 0.0; + destruction_dM_fe_silicate[i] = 0.0; + destruction_dM_carbon[i] = 0.0; + + if (itmask[i] != MASK_FALSE) { + double rho_gas = d(i, idx_range.j, idx_range.k); + double rho_dust_mg_sil_i = dust_mg_sil(i, idx_range.j, idx_range.k); + double rho_dust_fe_sil_i = dust_fe_sil(i, idx_range.j, idx_range.k); + double rho_dust_carb_i = dust_carb(i, idx_range.j, idx_range.k); + + bool mg_sil_active = (rho_dust_mg_sil_i >= dust_gate_threshold * rho_gas); + bool fe_sil_active = (rho_dust_fe_sil_i >= dust_gate_threshold * rho_gas); + bool carb_active = (rho_dust_carb_i >= dust_gate_threshold * rho_gas); + if (!mg_sil_active && !fe_sil_active && !carb_active) continue; + if (rho_gas <= 0.0) continue; + + double sne_this = use_sne ? sne(i, idx_range.j, idx_range.k) : 0.0; + double temp = std::max(t_gas[i], tiny_value); + double dt = dt_value[i]; + if (dt <= 0.0) continue; + + double rho_H_nuclei = h_census.rho_H(i, idx_range.j, idx_range.k); + if (rho_H_nuclei <= 0.0) continue; + double nH = rho_H_nuclei * dens_proper / mh; + if (nH <= 0.0) continue; + + // Common (species-independent) sputtering structural factor. + double sput_struct = (my_chemistry->dust_grainsize / 0.1) * (1.0 / nH) * + (1.0 + std::pow(temp / colibre_sputter_t_ref, -2.5)); + + // Equivalent per-time destruction rate for one species; its product + // with dt is an exponential e-folding mass loss. + auto compute_dM = [&](double rho_dust, double shock_factor) -> double { + double inv_tau_loss = 0.0; + + if (use_tau_dest) { + double tau_dest = tau_dest_field(i, idx_range.j, idx_range.k); + if (tau_dest <= 0) tau_dest = huge_value; + // Apply species multiplier on top of user-supplied tau_dest: + // higher shock_factor -> shorter tau -> faster destruction. + tau_dest = tau_dest / shock_factor; + if (tau_dest > 0.0 && std::isfinite(tau_dest)) { + inv_tau_loss += 1.0 / tau_dest; + } + } else if (use_sne && sne_this > 0.0) { + // sne_rate holds the number of SNe per code volume delivered over + // the external timestep dt_full; Ms100 * sne_this is then the + // shocked-gas mass per volume in dt_full, and dividing by + // (rho_gas * dt_full) gives the inverse shock-destruction + // timescale. The exponential update supplies the dt dependence. + double inv_tau_shock = Ms100 * shock_factor * sne_this * + my_chemistry->dust_destruction_eff / + (rho_gas * dt_full); + if (inv_tau_shock > 0.0 && std::isfinite(inv_tau_shock)) { + inv_tau_loss += inv_tau_shock; + } + } + + double tau_sput = tau_sput_ref * sput_struct; + if (tau_sput > 0.0 && std::isfinite(tau_sput)) { + inv_tau_loss += 1.0 / tau_sput; + } + + double dM = -e_fold_loss_rate(rho_dust, dt, inv_tau_loss); + if (std::isnan(dM)) { + std::cout << "dM (species) calculated as NaN, " << dM << std::endl; + dM = 0.0; + } + return dM; + }; + + if (carb_active) { + destruction_dM_carbon[i] = + compute_dM(rho_dust_carb_i, shock_factor_carbon); + } + if (mg_sil_active) { + destruction_dM_mg_silicate[i] = + compute_dM(rho_dust_mg_sil_i, shock_factor_silicate); + } + if (fe_sil_active) { + destruction_dM_fe_silicate[i] = + compute_dM(rho_dust_fe_sil_i, shock_factor_silicate); + } + } + } +} + +// ========================================== +// DUST UPDATE (SPECIES-SPECIFIC: Mg-silicate + Fe-silicate + carbonaceous) +// ========================================== +// Per-channel mass exchange between dust species and their corresponding +// gas-phase reactant pools. Active when dust_species_track == 1. +// +// carbon channel: rho_dust_carbonaceous <-> rho_metal_carbon +// Mg-sil channel: rho_dust_mg_silicate <-> {Mg, Si, O} as Mg2SiO4 +// Fe-sil channel: rho_dust_fe_silicate <-> {Fe, Si, O} as Fe2SiO4 +// +// Growth is capped by the limiting gas-phase reactant, destruction by the +// available dust mass. No SN injection on this path — host code seeds the +// species directly. +void grackle::impl::dust_update_species( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, const double* growth_dM_mg_silicate, + const double* growth_dM_fe_silicate, const double* growth_dM_carbon, + const double* destruction_dM_mg_silicate, + const double* destruction_dM_fe_silicate, + const double* destruction_dM_carbon, bool dryrun) { + 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]); + grackle::impl::View dust_sil( + my_fields->dust_density_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_mg_sil( + my_fields->dust_density_mg_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_fe_sil( + my_fields->dust_density_fe_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_carb( + my_fields->dust_density_carbonaceous, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + 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 mC( + my_fields->metal_density_carbon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mO( + my_fields->metal_density_oxygen, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mMg( + my_fields->metal_density_magnesium, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mSi( + my_fields->metal_density_silicon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View mFe( + my_fields->metal_density_iron, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + + const double f_mg_sil_Mg = my_chemistry->dust_mg_silicate_f_Mg; + const double f_mg_sil_Fe = my_chemistry->dust_mg_silicate_f_Fe; + const double f_mg_sil_Si = my_chemistry->dust_mg_silicate_f_Si; + const double f_mg_sil_O = my_chemistry->dust_mg_silicate_f_O; + const double f_fe_sil_Mg = my_chemistry->dust_fe_silicate_f_Mg; + const double f_fe_sil_Fe = my_chemistry->dust_fe_silicate_f_Fe; + const double f_fe_sil_Si = my_chemistry->dust_fe_silicate_f_Si; + const double f_fe_sil_O = my_chemistry->dust_fe_silicate_f_O; + + for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { + if (itmask[i] == MASK_FALSE) continue; + + double rho_gas = d(i, idx_range.j, idx_range.k); + double rho_dust_mg_sil_i = dust_mg_sil(i, idx_range.j, idx_range.k); + double rho_dust_fe_sil_i = dust_fe_sil(i, idx_range.j, idx_range.k); + double rho_dust_carb_i = dust_carb(i, idx_range.j, idx_range.k); + double rho_metal_total = metal(i, idx_range.j, idx_range.k); + double rho_C = mC(i, idx_range.j, idx_range.k); + double rho_O = mO(i, idx_range.j, idx_range.k); + double rho_Mg = mMg(i, idx_range.j, idx_range.k); + double rho_Si = mSi(i, idx_range.j, idx_range.k); + double rho_Fe = mFe(i, idx_range.j, idx_range.k); + double dt = dt_value[i]; + if (dt <= 0.0) continue; + + double rho_tracked_before = rho_C + rho_O + rho_Mg + rho_Si + rho_Fe; + // Untracked gas-phase metals (Al, S, Ne, Na, ...). If host code has seeded + // per-element fields independently and they exceed rho_metal_total, + // rho_metal_other would go negative and propagate into rho_metal_new on + // the rebuild below — clamp at zero to prevent silent mass leakage. + double rho_metal_other = + std::max(rho_metal_total - rho_tracked_before, 0.0); + + // dM > 0 means net growth (gas reactants -> dust); + // dM < 0 means net destruction (dust -> gas reactants). + double dM_carb = (growth_dM_carbon[i] + destruction_dM_carbon[i]) * dt; + double dM_mg_sil = + (growth_dM_mg_silicate[i] + destruction_dM_mg_silicate[i]) * dt; + double dM_fe_sil = + (growth_dM_fe_silicate[i] + destruction_dM_fe_silicate[i]) * dt; + + // ---------- Per-channel caps ---------- + // Carbon channel: bounded by rho_C (growth) or rho_dust_carb + // (destruction). Leave a 10% margin on growth so the gas-C reservoir is + // not driven to exactly zero in a single subcycle — make_consistent's + // gas-phase C correction (correctCg = Cg/totalCg) would otherwise zero + // CI/CII/CO/... and destabilize the metal-chemistry implicit solve that + // follows. + if (dM_carb > 0.0) { + dM_carb = std::min(dM_carb, 0.9 * rho_C); + } else { + dM_carb = std::max(dM_carb, -rho_dust_carb_i); + } + + // Mg- and Fe-silicate channels share Si/O, but Mg and Fe are independent. + // Cap destruction by each species' own dust mass. + if (dM_mg_sil < 0.0) { + dM_mg_sil = std::max(dM_mg_sil, -rho_dust_mg_sil_i); + } + if (dM_fe_sil < 0.0) { + dM_fe_sil = std::max(dM_fe_sil, -rho_dust_fe_sil_i); + } + + // Cap positive growth jointly against the gas reservoirs: each pass + // rescales only the channels that consume the most-limiting element, so + // Mg-silicate growth survives when Fe alone limits Fe-silicate (and vice + // versa) while shared Si/O is never double-used. The 0.9 safety margin + // (same rationale as the carbon channel) leaves headroom in each gas + // reservoir so make_consistent's per-element corrections don't divide by + // ~zero in the next subcycle. + const double sil_grow_safety = 0.9; + double dM_mg_sil_grow = std::max(dM_mg_sil, 0.0); + double dM_fe_sil_grow = std::max(dM_fe_sil, 0.0); + for (int cap_iter = 0; cap_iter < 4; cap_iter++) { + double scale = 1.0; + int limiter = -1; // 0=Mg, 1=Fe, 2=Si, 3=O + auto consider_element = [&](double rho_X, double c_mg_sil, + double c_fe_sil, int element_id) { + double need = dM_mg_sil_grow * c_mg_sil + dM_fe_sil_grow * c_fe_sil; + double budget = sil_grow_safety * rho_X; + if (need > budget && need > 0.0) { + double trial = std::max(0.0, budget / need); + if (trial < scale) { + scale = trial; + limiter = element_id; + } + } + }; + consider_element(rho_Mg, f_mg_sil_Mg, f_fe_sil_Mg, 0); + consider_element(rho_Fe, f_mg_sil_Fe, f_fe_sil_Fe, 1); + consider_element(rho_Si, f_mg_sil_Si, f_fe_sil_Si, 2); + consider_element(rho_O, f_mg_sil_O, f_fe_sil_O, 3); + if (limiter < 0 || scale >= 1.0) break; + + if (limiter == 0) { + if (f_mg_sil_Mg > 0.0) dM_mg_sil_grow *= scale; + if (f_fe_sil_Mg > 0.0) dM_fe_sil_grow *= scale; + } else if (limiter == 1) { + if (f_mg_sil_Fe > 0.0) dM_mg_sil_grow *= scale; + if (f_fe_sil_Fe > 0.0) dM_fe_sil_grow *= scale; + } else if (limiter == 2) { + if (f_mg_sil_Si > 0.0) dM_mg_sil_grow *= scale; + if (f_fe_sil_Si > 0.0) dM_fe_sil_grow *= scale; + } else { + if (f_mg_sil_O > 0.0) dM_mg_sil_grow *= scale; + if (f_fe_sil_O > 0.0) dM_fe_sil_grow *= scale; + } + } + if (dM_mg_sil > 0.0) { + dM_mg_sil = dM_mg_sil_grow; + } + if (dM_fe_sil > 0.0) { + dM_fe_sil = dM_fe_sil_grow; + } + + // ---------- Apply ---------- + rho_dust_carb_i += dM_carb; + rho_C -= dM_carb; + + rho_dust_mg_sil_i += dM_mg_sil; + rho_dust_fe_sil_i += dM_fe_sil; + rho_Mg -= dM_mg_sil * f_mg_sil_Mg + dM_fe_sil * f_fe_sil_Mg; + rho_Fe -= dM_mg_sil * f_mg_sil_Fe + dM_fe_sil * f_fe_sil_Fe; + rho_Si -= dM_mg_sil * f_mg_sil_Si + dM_fe_sil * f_fe_sil_Si; + rho_O -= dM_mg_sil * f_mg_sil_O + dM_fe_sil * f_fe_sil_O; + + // Floors / NaN guard + rho_dust_carb_i = std::max(0.0, rho_dust_carb_i); + rho_dust_mg_sil_i = std::max(0.0, rho_dust_mg_sil_i); + rho_dust_fe_sil_i = std::max(0.0, rho_dust_fe_sil_i); + rho_C = std::max(0.0, rho_C); + rho_O = std::max(0.0, rho_O); + rho_Mg = std::max(0.0, rho_Mg); + rho_Si = std::max(0.0, rho_Si); + rho_Fe = std::max(0.0, rho_Fe); + + // Bulk dust = Mg-silicate + Fe-silicate + carbonaceous. The compatibility + // silicate field is the sum of both silicate endmembers. + double rho_dust_sil_i = rho_dust_mg_sil_i + rho_dust_fe_sil_i; + double rho_dust_new = rho_dust_sil_i + rho_dust_carb_i; + + // Rebuild total metals from the updated tracked fields plus the unchanged + // untracked remainder; this stays consistent even if user-edited silicate + // fractions do not sum to exactly one. + double rho_tracked_after = rho_C + rho_O + rho_Mg + rho_Si + rho_Fe; + double rho_metal_new = std::max(rho_metal_other + rho_tracked_after, 0.0); + double delta_metal_total = rho_metal_new - rho_metal_total; + + // Gas density tracks metals as a subset (dust is not part of `density`). + rho_gas += delta_metal_total; + + if (std::isnan(rho_dust_new) || std::isnan(rho_metal_new) || + std::isnan(rho_gas)) { + std::cout << "dust_update_species: NaN at cell " << i + << " dM_carb=" << dM_carb << " dM_mg_sil=" << dM_mg_sil + << " dM_fe_sil=" << dM_fe_sil << std::endl; + continue; + } + + if (!dryrun) { + dust(i, idx_range.j, idx_range.k) = (gr_float)rho_dust_new; + dust_sil(i, idx_range.j, idx_range.k) = (gr_float)rho_dust_sil_i; + dust_mg_sil(i, idx_range.j, idx_range.k) = (gr_float)rho_dust_mg_sil_i; + dust_fe_sil(i, idx_range.j, idx_range.k) = (gr_float)rho_dust_fe_sil_i; + dust_carb(i, idx_range.j, idx_range.k) = (gr_float)rho_dust_carb_i; + metal(i, idx_range.j, idx_range.k) = (gr_float)rho_metal_new; + mC(i, idx_range.j, idx_range.k) = (gr_float)rho_C; + mO(i, idx_range.j, idx_range.k) = (gr_float)rho_O; + mMg(i, idx_range.j, idx_range.k) = (gr_float)rho_Mg; + mSi(i, idx_range.j, idx_range.k) = (gr_float)rho_Si; + mFe(i, idx_range.j, idx_range.k) = (gr_float)rho_Fe; + d(i, idx_range.j, idx_range.k) = (gr_float)rho_gas; + } + } +} + +// ========================================== +// DUST UPDATE (BULK) +// ========================================== +void grackle::impl::dust_update(chemistry_data* my_chemistry, + grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, + const gr_mask_type* itmask, + const double* dt_value, const double* growth_dM, + const double* destruction_dM, bool dryrun) { + 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]); + grackle::impl::View metal( + my_fields->metal_density, 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_stop; i++) { + if (itmask[i] == MASK_FALSE) continue; + + double rho_gas = d(i, idx_range.j, idx_range.k); + double rho_dust = dust(i, idx_range.j, idx_range.k); + double rho_metal_total = metal(i, idx_range.j, idx_range.k); + double dt = dt_value[i]; + + // dM_exchange > 0: net growth (metal -> dust); < 0: destruction + // (dust -> metal). Cap destruction by available dust and growth by 90% + // of gas-phase metals. + double dM_exchange = (growth_dM[i] + destruction_dM[i]) * dt; + dM_exchange = std::max(-rho_dust, dM_exchange); + dM_exchange = std::min(0.9 * rho_metal_total, dM_exchange); + + rho_dust += dM_exchange; + rho_metal_total -= dM_exchange; + rho_gas -= dM_exchange; // gas tracks metals as a subset + + rho_dust = std::max(0.0, rho_dust); + rho_metal_total = std::max(0.0, rho_metal_total); + + if (!dryrun) { + dust(i, idx_range.j, idx_range.k) = (gr_float)rho_dust; + metal(i, idx_range.j, idx_range.k) = (gr_float)rho_metal_total; + d(i, idx_range.j, idx_range.k) = (gr_float)rho_gas; + } + } +} diff --git a/src/clib/dust/dust_growth_and_destruction.hpp b/src/clib/dust/dust_growth_and_destruction.hpp new file mode 100644 index 000000000..777c77365 --- /dev/null +++ b/src/clib/dust/dust_growth_and_destruction.hpp @@ -0,0 +1,108 @@ +#ifndef DUST_GROWTH_AND_DESTRUCTION_HPP +#define DUST_GROWTH_AND_DESTRUCTION_HPP + +#include "grackle_types.h" +#include "grackle_chemistry_data.h" +#include "support/index_helper.hpp" +#include "internal_units.hpp" +#include "phys_constants.h" +#include "fortran_func_decls.h" + +namespace grackle::impl { + +// Calculates dust growth rates (accretion) onto grain surfaces. +// Stores the mass change rate for each cell in growth_dM. +void dust_growth(chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, + const gr_mask_type* itmask, const double* dt_value, + const double* t_gas, + double* growth_dM // output: mass change rate for each cell +); + +// Species-specific accretion onto three pre-existing dust populations +// (Mg-silicate + Fe-silicate + carbonaceous). Active when +// dust_species_track == 1. +// - carbonaceous: bottlenecked by gas-phase carbon +// - total silicate: bottlenecked by O, Si, and Mg+Fe, then split between +// Mg2SiO4 and Fe2SiO4 endmembers according to local Mg/Fe abundance +// Uses the COLIBRE dense-gas accretion normalization with clumped +// n_H' = C n_H, T_ref = 10 K, n_H_ref = 10 cm^-3, S_ref = 0.3, a_ref = 0.1 +// micron, and tau_ref values stored in Myr +// [REF: Trayford+2026 MNRAS 545, staf2040]. +// The returned rates are equivalent rates whose product with dt gives the +// exponential e-folding mass update. +void dust_growth_species( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, const double* t_gas, + double* growth_dM_mg_silicate, // output: Mg-silicate accretion rate + double* growth_dM_fe_silicate, // output: Fe-silicate accretion rate + double* growth_dM_carbon // output: carbonaceous accretion rate +); + +// Calculates dust destruction rates from SNe shocks and thermal sputtering. +// dt_full is the full external timestep (sne_rate is normalized per external +// step, not per subcycle). Stores the mass change rate for each cell in +// destruction_dM. +void dust_destruction( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, double dt_full, const double* t_gas, + double* destruction_dM // output: mass change rate for each cell +); + +// Species-specific destruction (SN shocks + thermal sputtering) of the three +// dust populations. Active when dust_species_track == 1. +// - shock yield: graphite is baseline (factor 1.0); both silicate endmembers +// follow the Slavin+2015 standard SNR gas-cleared mass ratio +// 990/600 = 1.65 +// [REF: Slavin, Dwek, Jones 2015 ApJ 803, 7; Jones+1996 ApJ 469, 740] +// - thermal sputtering: common COLIBRE/Tsai-Mathews tau_ref for all species +// with (a/0.1) (n_H/cm^-3)^-1 [1 + (T/2e6 K)^-2.5] scaling +// [REF: Tsai & Mathews 1995 ApJ 448, 84] +// Destruction uses exponential decay; the returned rates are equivalent rates +// whose product with dt gives the e-folding mass loss. +void dust_destruction_species( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, double dt_full, const double* t_gas, + double* destruction_dM_mg_silicate, // output: Mg-silicate destruction rate + double* destruction_dM_fe_silicate, // output: Fe-silicate destruction rate + double* destruction_dM_carbon // output: carbonaceous destruction rate +); + +// Update the density fields using calculated mass changes. +void dust_update( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, + const double* growth_dM, // input: mass change from growth + const double* destruction_dM, // input: mass change from destruction + bool dryrun); + +// Species-specific field update for the split-silicate path +// (dust_species_track == 1). Per-channel mass exchange: +// - carbon channel: rho_dust_carbonaceous <-> metal_density_carbon +// - Mg-sil channel: rho_dust_mg_silicate <-> {Mg, Si, O} as Mg2SiO4 +// - Fe-sil channel: rho_dust_fe_silicate <-> {Fe, Si, O} as Fe2SiO4 +// Growth is capped by the limiting gas-phase reactant, destruction by the +// available dust mass. No SN injection here — host code seeds dust species +// directly (or via the inject_pathway machinery in make_consistent). +void dust_update_species( + chemistry_data* my_chemistry, grackle_field_data* my_fields, + InternalGrUnits internalu, IndexRange idx_range, const gr_mask_type* itmask, + const double* dt_value, + const double* growth_dM_mg_silicate, // input: Mg-silicate accretion rate + const double* growth_dM_fe_silicate, // input: Fe-silicate accretion rate + const double* growth_dM_carbon, // input: carbonaceous accretion rate + const double* destruction_dM_mg_silicate, // input: Mg-silicate destruction + // rate (<=0) + const double* destruction_dM_fe_silicate, // input: Fe-silicate destruction + // rate (<=0) + const double* + destruction_dM_carbon, // input: carbonaceous destruction rate (<=0) + bool dryrun); + +} // namespace grackle::impl + +#endif // DUST_GROWTH_AND_DESTRUCTION_HPP diff --git a/src/clib/dust/lookup_dust_rates1d.hpp b/src/clib/dust/lookup_dust_rates1d.hpp index 5e6be9a40..7fec65092 100644 --- a/src/clib/dust/lookup_dust_rates1d.hpp +++ b/src/clib/dust/lookup_dust_rates1d.hpp @@ -196,10 +196,13 @@ inline void lookup_dust_rates1d( GRIMPL_REQUIRE(my_chemistry->metal_chemistry == 1, "sanity-check!"); // Compute grain size increment - calc_grain_size_increment_1d(dom, idx_range, itmask_metal, my_chemistry, - my_rates->opaque_storage->grain_species_info, - my_rates->opaque_storage->inject_pathway_props, - my_fields, internal_dust_prop_scratch_buf); + if (my_chemistry->dust_model == 0) { + calc_grain_size_increment_1d( + dom, idx_range, itmask_metal, my_chemistry, + my_rates->opaque_storage->grain_species_info, + my_rates->opaque_storage->inject_pathway_props, my_fields, + internal_dust_prop_scratch_buf); + } grackle::impl::View d( my_fields->density, my_fields->grid_dimension[0], diff --git a/src/clib/dust_growth_and_destruction.F b/src/clib/dust_growth_and_destruction.F deleted file mode 100644 index bd74989aa..000000000 --- a/src/clib/dust_growth_and_destruction.F +++ /dev/null @@ -1,361 +0,0 @@ -#include "phys_const.def" -#include "dust_evol.def" - -!//////////////////////////////////////////////////////////// -! calculate characteristic timescales used for subcycling -!//////////////////////////////////////////////////////////// - subroutine dust_tau_for_timestep(in, jn, kn, i, j, k, - & Density, Tgas, Tdust, metal_Density, - & gas_metal1, gas_metal2, gas_metal3, - & gas_metal4, gas_metal5, gas_metal6, - & gas_metal7, gas_metal8, gas_metal9, - & gas_metal10, SNe_ThisTimeStep, - & dust_destruction_eff, - & sne_coeff, sne_shockspeed, - & dust_grainsize, dust_growth_densref, - & dust_growth_tauref, - & tau_accr, tau_dest, tau_sput, - & z_solar, SolarAbundances, - & time_units, dens_units, len_units, aye, dtime) - implicit NONE -#include "grackle_fortran_types.def" - - character(80) :: userinp - integer in, jn, kn, i, j, k - real*8 tau_accr(in), tau_dest(in), tau_sput(in), - & z_solar, SolarAbundances(NUM_METAL_SPECIES_GRACKLE), - & dust_destruction_eff, sne_coeff, sne_shockspeed, - & dust_grainsize, dust_growth_densref, - & dust_growth_tauref, - & time_units, dens_units, len_units, aye, dtime - R_PREC Density(in,jn,kn), Tgas(in,jn,kn), Tdust(in,jn,kn), metal_Density(in,jn,kn), - & gas_metal1(in,jn,kn), gas_metal2(in,jn,kn), - & gas_metal3(in,jn,kn), gas_metal4(in,jn,kn), - & gas_metal5(in,jn,kn), gas_metal6(in,jn,kn), - & gas_metal7(in,jn,kn), gas_metal8(in,jn,kn), - & gas_metal9(in,jn,kn), gas_metal10(in,jn,kn), - & SNe_ThisTimeStep(in,jn,kn) - -! local variables - real*8 tau_ref, Ms100, tau_accr0, dens_proper - - dens_proper = dens_units * aye**3 - - Ms100 = 6800.0*sne_coeff - & *(100.0/sne_shockspeed)*(100.0/sne_shockspeed) - & * SolarMass / (dens_units * len_units**3) ! gas mass shocked per SNe (Sedov-Taylor) -! check: - tau_ref = dust_growth_tauref * 1e9 * SEC_PER_YEAR_GRACKLE / time_units - -! growth - tau_accr0 = tau_ref * (dust_growth_densref / dens_proper) - & * (Tdust(i,j,k)/Tgas(i,j,k))**0.5 - tau_accr(i) = 1.d+20 - tau_accr(i) = min(tau_accr0 * (z_solar - & / metal_Density(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(1) - & / gas_metal1(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(2) - & / gas_metal2(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(3) - & / gas_metal3(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(4) - & / gas_metal4(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(5) - & / gas_metal5(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(6) - & / gas_metal6(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(7) - & / gas_metal7(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(8) - & / gas_metal8(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(9) - & / gas_metal9(i,j,k)), tau_accr(i)) - tau_accr(i) = min(tau_accr0 * (SolarAbundances(10) - & / gas_metal10(i,j,k)), tau_accr(i)) - - if (tau_accr(i) .le. 0) then - write(*,*) "-- tau_accr neagative! --", tau_accr(i) - write(*,*) "tau_accr0=", tau_accr0 - write(*,*) "metal_density=", metal_density(i,j,k) - write(*,*) "gas_metal1=", gas_metal1(i,j,k) - write(*,*) "gas_metal2=", gas_metal2(i,j,k) - write(*,*) "gas_metal3=", gas_metal3(i,j,k) - write(*,*) "gas_metal4=", gas_metal4(i,j,k) - write(*,*) "gas_metal5=", gas_metal5(i,j,k) - write(*,*) "gas_metal6=", gas_metal6(i,j,k) - write(*,*) "gas_metal7=", gas_metal7(i,j,k) - write(*,*) "gas_metal8=", gas_metal8(i,j,k) - write(*,*) "gas_metal9=", gas_metal9(i,j,k) - write(*,*) "gas_metal10=", gas_metal10(i,j,k) - write(*,*) "SolarAbundances1=", SolarAbundances(1) - write(*,*) "SolarAbundances2=", SolarAbundances(2) - write(*,*) "SolarAbundances3=", SolarAbundances(3) - write(*,*) "SolarAbundances4=", SolarAbundances(4) - write(*,*) "SolarAbundances5=", SolarAbundances(5) - write(*,*) "SolarAbundances6=", SolarAbundances(6) - write(*,*) "SolarAbundances7=", SolarAbundances(7) - write(*,*) "SolarAbundances8=", SolarAbundances(8) - write(*,*) "SolarAbundances9=", SolarAbundances(9) - write(*,*) "SolarAbundances10=", SolarAbundances(10) - call exit(74) - endif - -! destruction by SN shocks - if (SNe_ThisTimeStep(i,j,k).le.0.0) then - tau_dest(i) = 1.d+20 - else - tau_dest(i) = Density(i,j,k)/(Ms100*SNe_ThisTimeStep(i,j,k) - & *dust_destruction_eff) * dtime - endif - -! destruction by thermal sputtering - tau_sput(i) = 1.7e8 * SEC_PER_YEAR_GRACKLE / time_units - & * (dust_grainsize/0.1) * (1e-27/(dens_proper*Density(i,j,k))) - & * ((2e6/Tgas(i,j,k))**2.5+1.0) ! Tsai & Mathews (1995) - - return - end ! end of dust_tau_for_timestep - - -!//////////////////////////////////////////////////////////// -! main routine for dust growth and destruction -!//////////////////////////////////////////////////////////// - subroutine dust_growth_and_destruction(in, jn, kn, i, j, k, - & gas_Density, dust_Density, metal_Density, Tgas, Tdust, - & gas_metal1, gas_metal2, gas_metal3, - & gas_metal4, gas_metal5, gas_metal6, - & gas_metal7, gas_metal8, gas_metal9, - & gas_metal10, - & dust_metal1, dust_metal2, dust_metal3, - & dust_metal4, dust_metal5, dust_metal6, - & dust_metal7, dust_metal8, dust_metal9, - & dust_metal10, - & SNe_ThisTimeStep, dust_destruction_eff, - & sne_coeff, sne_shockspeed, - & dust_grainsize, dust_growth_densref, - & dust_growth_tauref, - & tau_dest, tau_sput, - & z_solar, SolarAbundances, - & time_units, dens_units, len_units, aye, dtime, metpt, dmetpt) - implicit NONE -#include "grackle_fortran_types.def" -! arguments - integer in, jn, kn, i, j, k - real*8 tau_dest(in), tau_sput(in), - & z_solar, SolarAbundances(NUM_METAL_SPECIES_GRACKLE), - & dust_destruction_eff, sne_coeff, sne_shockspeed, - & dust_grainsize, dust_growth_densref, - & dust_growth_tauref, - & time_units, dens_units, len_units, aye, dtime, dM_conserv - R_PREC Tgas(in,jn,kn), Tdust(in,jn,kn), - & gas_metal1(in,jn,kn), gas_metal2(in,jn,kn), - & gas_metal3(in,jn,kn), gas_metal4(in,jn,kn), - & gas_metal5(in,jn,kn), gas_metal6(in,jn,kn), - & gas_metal7(in,jn,kn), gas_metal8(in,jn,kn), - & gas_metal9(in,jn,kn), gas_metal10(in,jn,kn), - & dust_metal1(in,jn,kn), dust_metal2(in,jn,kn), - & dust_metal3(in,jn,kn), dust_metal4(in,jn,kn), - & dust_metal5(in,jn,kn), dust_metal6(in,jn,kn), - & dust_metal7(in,jn,kn), dust_metal8(in,jn,kn), - & dust_metal9(in,jn,kn), dust_metal10(in,jn,kn), - & metal_Density(in,jn,kn), dust_Density(in,jn,kn), - & gas_Density(in,jn,kn), SNe_ThisTimeStep(in,jn,kn), - & metpt(NUM_METAL_SPECIES_GRACKLE+1), dmetpt(NUM_METAL_SPECIES_GRACKLE+1) - -! local variables - integer n - character(80) :: userinp - real*8 tau_ref, Ms100, tau_accr0, tau_accr - R_PREC dM(NUM_METAL_SPECIES_GRACKLE+1), dMs, - & Mmet(NUM_METAL_SPECIES_GRACKLE+1), - & dens_proper - -! debug variables - real*8 dM_tau_accr - R_PREC total_density_init, total_density_final - - dens_proper = dens_units * aye**3 - -! copy metallicites into a single array - metpt(1) = metal_Density(i,j,k) - metpt(2) = gas_metal1(i,j,k) - metpt(3) = gas_metal2(i,j,k) - metpt(4) = gas_metal3(i,j,k) - metpt(5) = gas_metal4(i,j,k) - metpt(6) = gas_metal5(i,j,k) - metpt(7) = gas_metal6(i,j,k) - metpt(8) = gas_metal7(i,j,k) - metpt(9) = gas_metal8(i,j,k) - metpt(10) = gas_metal9(i,j,k) - metpt(11) = gas_metal10(i,j,k) - dmetpt(1) = dust_Density(i,j,k) - dmetpt(2) = dust_metal1(i,j,k) - dmetpt(3) = dust_metal2(i,j,k) - dmetpt(4) = dust_metal3(i,j,k) - dmetpt(5) = dust_metal4(i,j,k) - dmetpt(6) = dust_metal5(i,j,k) - dmetpt(7) = dust_metal6(i,j,k) - dmetpt(8) = dust_metal7(i,j,k) - dmetpt(9) = dust_metal8(i,j,k) - dmetpt(10) = dust_metal9(i,j,k) - dmetpt(11) = dust_metal10(i,j,k) - - total_density_init = metpt(1) + dmetpt(1) - - Ms100 = 6800.0*sne_coeff - & *(100.0/sne_shockspeed)*(100.0/sne_shockspeed) - & * SolarMass / (dens_units * len_units**3) ! gas mass shocked per SNe (Sedov-Taylor) - tau_ref = dust_growth_tauref * 1e9 * SEC_PER_YEAR_GRACKLE / time_units - - do n=0+1,NUM_METAL_SPECIES_GRACKLE+1 - dM(n) = 0. - enddo - -! metals accreted from gas to dust - tau_accr0 = tau_ref * (dust_growth_densref/dens_proper) - & * (Tdust(i,j,k)/Tgas(i,j,k))**0.5 - do n=1+1,NUM_METAL_SPECIES_GRACKLE+1 - Mmet(n) = metpt(n) + dmetpt(n) - tau_accr = tau_accr0 * SolarAbundances(n-1) / metpt(n) - - if (dmetpt(n) .ne. dmetpt(n)) then - write(*,*) "dmetpt at index", n, - & "calculated as NaN in line 174 of dust model", dmetpt(n) - endif - - if (Mmet(n) .ne. Mmet(n)) then - write(*,*) "Mmet at index", n, - & "calculated as NaN in line 174 of dust model", Mmet(n) - endif - - if (tau_accr .ne. tau_accr) then - write(*,*) "tau_accr", - & "calculated as NaN in line 174 of dust model", tau_accr - endif - - if (dtime .ne. dtime) then - write(*,*) "dtime", - & "calculated as NaN in line 174 of dust model", dtime - endif - - dM(n) = dM(n) + min((1 - dmetpt(n) / Mmet(n)) - & * (dmetpt(n)/tau_accr) * dtime, - & Mmet(n) - dmetpt(n)) ! growth capped by gas in a single cell - - if (dM(n) .ne. dM(n)) then - write(*,*) "dM at index", n, - & "calculated as NaN in line 174 of dust model", dM(n) - endif - - dM(1) = dM(1) + dM(n) - enddo - - dM_tau_accr = dM(1) - -! SNe shock and thermal sputtering - if (SNe_ThisTimeStep(i,j,k) .le. 0.0) then - dMs = 0.0 - else - dMs = min(dust_Density(i,j,k)/tau_dest(i)*dtime, - & dust_Density(i,j,k)) - endif - if (dMs .ge. dust_Density(i,j,k)) then - if (dMs .gt. dust_Density(i,j,k)) then - write (*,*) "WARNING: dMs > Mdust SNe shock destruction:", - & SNe_ThisTimeStep(i,j,k), tau_dest(i) - endif - else - dMs = dMs + dust_Density(i,j,k) / tau_sput(i) *3.0* dtime - dMs = min(dMs,dust_Density(i,j,k)) - endif - do n=0+1,NUM_METAL_SPECIES_GRACKLE+1 - dM(n) = dM(n) - dmetpt(n) * dMs - - if (dM(n) .ne. dM(n)) then - write(*,*) "dM at index", n, - & "calculated as NaN in line 202 of dust model", dM(n) - endif - enddo - -! recalculate metallicity - do n=0+1,NUM_METAL_SPECIES_GRACKLE+1 -! 1. Place conditions on dM such that mass is conserved but dust and gas -! mass cannot be negative -! 2. Make sure that only 90% of the gaseous metals can become dust - dM(n) = max(-1*dmetpt(n), dM(n)) - dM(n) = min(0.9 * metpt(n), dM(n)) - -! Transfer metals between gaseous and dust phase only if dust present - dM_conserv = 0.0 - if (dust_Density(i,j,k) .gt. 0.0) then - dmetpt(n) = dmetpt(n) + dM(n) - metpt(n) = metpt(n) - dM(n) - else -! If dust density is <= 0 then makes sure metals are conserved accordingly - dM_conserv = dmetpt(n) - dmetpt(n) = dmetpt(n) - dM_conserv - metpt(n) = metpt(n) + dM_conserv - endif -! surely this should be only happen if you add to the dust metallicity, -! else you are removing gas metallicity and breaking conservation - ! metpt(n) = metpt(n) - dM(n) - enddo - -! update gas density to account for changes in metal density - gas_Density(i,j,k) = gas_Density(i,j,k) + (metpt(1) - metal_Density(i,j,k)) - - if (metpt(1) .lt. 0) then - write(*,*) "metpt(1)=", metpt(1), - & "dM(1)=", dM(1) - endif - - if (dmetpt(1) .lt. 0) then - write(*,*) "final metal density = ", metpt(1) - write(*,*) "inital metal density = ", metal_Density(i,j,k) - write(*,*) "DELTA metal density = ", -1*dM(1) - write(*,*) "final dust density = ", dmetpt(1) - write(*,*) "initial dust density = ", dust_Density(i,j,k) - write(*,*) "DELTA dust density = ", dM(1) - write(*,*) "------------" - write(*,*) "dM_tau_accr = ", dM_tau_accr - call exit(21) - endif - - total_density_final = metpt(1) + dmetpt(1) - - if (abs(total_density_final - total_density_init) .gt. 1d-8) then - write(*,*) "Total dust + metal density not conserved!" - write(*,*) "Initial value --> ", total_density_init - write(*,*) "Final value --> ", total_density_final - write(*,*) "dmetpt = ", dmetpt(1), ", metpt = ", metpt(1), ", dM = ", dM(1) - write(*,*) "Ending run!!" - call exit(21) - endif - -! copy metallicites back into fields arrays - metal_Density(i,j,k) = metpt(1) - gas_metal1(i,j,k) = metpt(2) - gas_metal2(i,j,k) = metpt(3) - gas_metal3(i,j,k) = metpt(4) - gas_metal4(i,j,k) = metpt(5) - gas_metal5(i,j,k) = metpt(6) - gas_metal6(i,j,k) = metpt(7) - gas_metal7(i,j,k) = metpt(8) - gas_metal8(i,j,k) = metpt(9) - gas_metal9(i,j,k) = metpt(10) - gas_metal10(i,j,k) = metpt(11) - dust_Density(i,j,k) = dmetpt(1) - dust_metal1(i,j,k) = dmetpt(2) - dust_metal2(i,j,k) = dmetpt(3) - dust_metal3(i,j,k) = dmetpt(4) - dust_metal4(i,j,k) = dmetpt(5) - dust_metal5(i,j,k) = dmetpt(6) - dust_metal6(i,j,k) = dmetpt(7) - dust_metal7(i,j,k) = dmetpt(8) - dust_metal8(i,j,k) = dmetpt(9) - dust_metal9(i,j,k) = dmetpt(10) - dust_metal10(i,j,k) = dmetpt(11) - - return - end ! end of dust growth and destruction module diff --git a/src/clib/field_data_misc_fdatamembers.def b/src/clib/field_data_misc_fdatamembers.def index 3832322b2..057097d0a 100644 --- a/src/clib/field_data_misc_fdatamembers.def +++ b/src/clib/field_data_misc_fdatamembers.def @@ -28,9 +28,24 @@ ENTRY(z_velocity) // metal_cooling = 1 ENTRY(metal_density) + // dust_species_track = 1 +ENTRY(metal_density_carbon) +ENTRY(metal_density_oxygen) + + // dust_species_track = 1 (5-element gas tracking) +ENTRY(metal_density_magnesium) +ENTRY(metal_density_silicon) +ENTRY(metal_density_iron) + // use_dust_density_field = 1 ENTRY(dust_density) + // dust_species_track = 1 (Mg-silicate + Fe-silicate + carbonaceous dust) +ENTRY(dust_density_silicate) +ENTRY(dust_density_mg_silicate) +ENTRY(dust_density_fe_silicate) +ENTRY(dust_density_carbonaceous) + // use_volumetric_heating_rate = 1 ENTRY(volumetric_heating_rate) // use_specific_heating_rate = 1 @@ -79,3 +94,6 @@ ENTRY(Al2O3_dust_temperature) ENTRY(ref_org_dust_temperature) ENTRY(vol_org_dust_temperature) ENTRY(H2O_ice_dust_temperature) + + // dust model parameter + ENTRY(sne_rate) diff --git a/src/clib/grackle_chemistry_data_fields.def b/src/clib/grackle_chemistry_data_fields.def index 7c2887df9..245bb9d92 100644 --- a/src/clib/grackle_chemistry_data_fields.def +++ b/src/clib/grackle_chemistry_data_fields.def @@ -333,3 +333,71 @@ ENTRY(omp_nthreads, INT, 0) * 3: Newton-Raphson-only */ ENTRY(solver_method, INT, 1) + +/* dust model selection +* 0: default (runs calc_grain_size_increment_1d as usual) +* 1: harrison dust model +*/ +ENTRY(dust_model, INT, 0) + +/* Flag to use snetimestep */ +ENTRY(use_sne_field, INT, 0) + +/* Flag to use user-provided dust destruction timescale */ +ENTRY(use_tau_dest_field, INT, 0) + +/* Li et al 2019 dust model parameters*/ +ENTRY(dust_destruction_eff, DOUBLE, 3.0e-1) +ENTRY(sne_coeff, DOUBLE, 1.0) +ENTRY(sne_shockspeed, DOUBLE, 1.0e2) +ENTRY(dust_grainsize, DOUBLE, 1.0e-1) +ENTRY(dust_growth_densref, DOUBLE, 2.3e-22) +ENTRY(dust_growth_tauref, DOUBLE, 0.004) +ENTRY(dust_growth_sticking_coeff, DOUBLE, 3.0e-1) + +/* Dust creation by stellar feedback parameters */ +ENTRY(dust_condensation_eff, DOUBLE, 1.5e-1) +ENTRY(sne_metal_yield, DOUBLE, 3.0) + +/* Species-resolved dust tracking. + Requires dust_model=1. + 0) off (default — bulk dust_density) + 1) on — evolves dust_density_mg_silicate, dust_density_fe_silicate, + dust_density_carbonaceous and 5-element gas tracking + (C, O, Mg, Si, Fe). dust_density_silicate is maintained as + Mg-silicate + Fe-silicate for compatibility. + REF: Trayford+2026 MNRAS 545, staf2040 */ +ENTRY(dust_species_track, INT, 0) + +/* COLIBRE-style growth reference timescales [Myr]. + Normalized at n_H' = 10 cm^-3, T = 10 K, S = 0.3, a = 0.1 micron, + and solar bottleneck abundance. REF: Trayford+2026 Table 2, eq. 4. */ +ENTRY(dust_growth_tauref_silicate, DOUBLE, 99.3) +ENTRY(dust_growth_tauref_carbon, DOUBLE, 180.0) + +/* COLIBRE-style unresolved-density boost for grain growth. + C rises from 1 to 100 between n_H = 0.1 and 100 cm^-3, clipped outside. */ +ENTRY(dust_growth_clumping_factor_max, DOUBLE, 100.0) +ENTRY(dust_growth_clumping_nH_min, DOUBLE, 1.0e-1) +ENTRY(dust_growth_clumping_nH_max, DOUBLE, 1.0e2) + +/* COLIBRE / Tsai-Mathews common thermal sputtering reference time [Myr]. + tau_sp = tau_ref * (a/0.1) * (n_H/1 cm^-3)^-1 + * [1 + (T/2e6 K)^-2.5]. */ +ENTRY(dust_sputter_tauref, DOUBLE, 0.85) + +/* Initial/fallback mass fraction of silicate dust in Mg-silicate. COLIBRE + seeds equal numbers of Mg2SiO4 and Fe2SiO4 molecules, not equal masses. */ +ENTRY(dust_silicate_mg_fraction, DOUBLE, 0.408428) + +/* Mg-rich silicate endmember stoichiometric mass fractions for Mg2SiO4. */ +ENTRY(dust_mg_silicate_f_Mg, DOUBLE, 0.345521) +ENTRY(dust_mg_silicate_f_Fe, DOUBLE, 0.0) +ENTRY(dust_mg_silicate_f_Si, DOUBLE, 0.199630) +ENTRY(dust_mg_silicate_f_O, DOUBLE, 0.454849) + +/* Fe-rich silicate endmember stoichiometric mass fractions for Fe2SiO4. */ +ENTRY(dust_fe_silicate_f_Mg, DOUBLE, 0.0) +ENTRY(dust_fe_silicate_f_Fe, DOUBLE, 0.548111) +ENTRY(dust_fe_silicate_f_Si, DOUBLE, 0.137825) +ENTRY(dust_fe_silicate_f_O, DOUBLE, 0.314064) diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index bfa80a354..030ad43a1 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -238,6 +238,34 @@ static int local_initialize_chemistry_data_( return GR_FAIL; } + // dust_species_track=1 is the new species path (Mg-silicate + Fe-silicate + // + carbonaceous) and is mutually exclusive with the legacy grain_growth / + // dust_sublimation correction loop in make_consistent — that block scales + // legacy dust fields (MgSiO3, AC, Mg2SiO4, Fe3O4, SiO2D, MgO, FeS, Al2O3, + // SiM, FeM) against per-element ratios that the Phase F hijack has just + // overwritten with values derived from the new species fields. Running + // both paths simultaneously silently corrupts both reservoirs. + if (my_chemistry->dust_species_track == 1 && + my_chemistry->grain_growth == 1) { + fprintf(stderr, + "ERROR: dust_species_track = 1 is mutually exclusive with " + "grain_growth = 1 (set grain_growth = 0).\n"); + return GR_FAIL; + } + if (my_chemistry->dust_species_track == 1 && + my_chemistry->dust_sublimation == 1) { + fprintf(stderr, + "ERROR: dust_species_track = 1 is mutually exclusive with " + "dust_sublimation = 1 (set dust_sublimation = 0).\n"); + return GR_FAIL; + } + if (my_chemistry->dust_species_track == 1 && + my_chemistry->dust_model != 1) { + fprintf(stderr, + "ERROR: dust_species_track = 1 requires dust_model = 1.\n"); + return GR_FAIL; + } + // Default photo-electric heating to off if unset. if (my_chemistry->photoelectric_heating < 0) { my_chemistry->photoelectric_heating = 0; diff --git a/src/clib/make_consistent.cpp b/src/clib/make_consistent.cpp index 77b253893..c5f2fc3cb 100644 --- a/src/clib/make_consistent.cpp +++ b/src/clib/make_consistent.cpp @@ -13,7 +13,9 @@ // This file was initially generated automatically during conversion of the // make_consistent_g function from FORTRAN to C++ +#include #include +#include #include #include "grackle.h" @@ -59,6 +61,7 @@ void make_consistent( const_cast(my_fields->metal_density), my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View HM( my_fields->HM_density, my_fields->grid_dimension[0], my_fields->grid_dimension[1], my_fields->grid_dimension[2]); @@ -234,6 +237,52 @@ void make_consistent( } } + // Phase F hijack: when dust_species_track == 1, the per-element + // gas/dust/total nuclide arrays (Cg/Cd/Ct etc.) are populated directly + // from our tracked metal_density_X + dust species fields, overriding the + // SN-yield-derived computation below. Views are constructed once here. + grackle::impl::View mC_view, mO_view, mMg_view, mSi_view, + mFe_view, dust_sil_view, dust_mg_sil_view, dust_fe_sil_view, + dust_carb_view; + if (my_chemistry->dust_species_track == 1) { + mC_view = grackle::impl::View( + const_cast(my_fields->metal_density_carbon), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + mO_view = grackle::impl::View( + const_cast(my_fields->metal_density_oxygen), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + mMg_view = grackle::impl::View( + const_cast(my_fields->metal_density_magnesium), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + mSi_view = grackle::impl::View( + const_cast(my_fields->metal_density_silicon), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + mFe_view = grackle::impl::View( + const_cast(my_fields->metal_density_iron), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + dust_sil_view = grackle::impl::View( + const_cast(my_fields->dust_density_silicate), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + dust_mg_sil_view = grackle::impl::View( + const_cast(my_fields->dust_density_mg_silicate), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + dust_fe_sil_view = grackle::impl::View( + const_cast(my_fields->dust_density_fe_silicate), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + dust_carb_view = grackle::impl::View( + const_cast(my_fields->dust_density_carbonaceous), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + } + std::vector Ct(my_fields->grid_dimension[0]); std::vector Ot(my_fields->grid_dimension[0]); std::vector Mgt(my_fields->grid_dimension[0]); @@ -379,6 +428,55 @@ void make_consistent( Fed[i] = Fet[i] - Feg[i]; } + // Phase F hijack: replace SN-yield-derived per-element arrays with + // values pulled directly from our tracked fields. Carbonaceous dust + // is pure C; Mg- and Fe-silicates contribute Mg/Fe/Si/O at their own + // stoichiometric mass fractions. Al and S are not part of this + // architecture, so their per-element arrays are zeroed. + if (my_chemistry->dust_species_track == 1) { + const double f_mg_sil = + std::clamp(my_chemistry->dust_silicate_mg_fraction, 0.0, 1.0); + const double f_mg_sil_Mg = my_chemistry->dust_mg_silicate_f_Mg; + const double f_mg_sil_Fe = my_chemistry->dust_mg_silicate_f_Fe; + const double f_mg_sil_Si = my_chemistry->dust_mg_silicate_f_Si; + const double f_mg_sil_O = my_chemistry->dust_mg_silicate_f_O; + const double f_fe_sil_Mg = my_chemistry->dust_fe_silicate_f_Mg; + const double f_fe_sil_Fe = my_chemistry->dust_fe_silicate_f_Fe; + const double f_fe_sil_Si = my_chemistry->dust_fe_silicate_f_Si; + const double f_fe_sil_O = my_chemistry->dust_fe_silicate_f_O; + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { + double mg_sil = dust_mg_sil_view(i, j, k); + double fe_sil = dust_fe_sil_view(i, j, k); + const double sil_compat = dust_sil_view(i, j, k); + if (mg_sil + fe_sil <= 0.0 && sil_compat > 0.0) { + mg_sil = f_mg_sil * sil_compat; + fe_sil = (1.0 - f_mg_sil) * sil_compat; + } + const double carb = dust_carb_view(i, j, k); + Cg[i] = mC_view(i, j, k); + Og[i] = mO_view(i, j, k); + Mgg[i] = mMg_view(i, j, k); + Sig[i] = mSi_view(i, j, k); + Feg[i] = mFe_view(i, j, k); + Cd[i] = carb; + Od[i] = mg_sil * f_mg_sil_O + fe_sil * f_fe_sil_O; + Mgd[i] = mg_sil * f_mg_sil_Mg + fe_sil * f_fe_sil_Mg; + Sid[i] = mg_sil * f_mg_sil_Si + fe_sil * f_fe_sil_Si; + Fed[i] = mg_sil * f_mg_sil_Fe + fe_sil * f_fe_sil_Fe; + Ct[i] = Cg[i] + Cd[i]; + Ot[i] = Og[i] + Od[i]; + Mgt[i] = Mgg[i] + Mgd[i]; + Sit[i] = Sig[i] + Sid[i]; + Fet[i] = Feg[i] + Fed[i]; + Alt[i] = 0.; + Alg[i] = 0.; + Ald[i] = 0.; + St[i] = 0.; + Sg[i] = 0.; + Sd[i] = 0.; + } + } + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { // if (itmask_metal(i)) then OH(i, j, k) = std::fabs(OH(i, j, k)); @@ -522,7 +620,8 @@ void make_consistent( // ! if (d(i,j,k)*dom .lt. // ! & min(1.e6_DKIND/(metal(i,j,k)/d(i,j,k)/0.02d-4)**2 // ! & ,1.e6_DKIND)) then - if (((imetal == 0) && (d(i, j, k) * dom < 1.e8)) || + if (my_chemistry->dust_species_track || + ((imetal == 0) && (d(i, j, k) * dom < 1.e8)) || ((imetal == 1) && (((metal(i, j, k) <= 1.e-9 * d(i, j, k)) && (d(i, j, k) * dom < 1.e8)) || ((metal(i, j, k) > 1.e-9 * d(i, j, k)) && @@ -979,6 +1078,51 @@ void make_consistent( } } + // Phase E invariant: bulk dust_density = Mg-silicate + Fe-silicate + // + carbonaceous, with dust_density_silicate maintained as the silicate sum. + // dust_update_species() maintains this per-cell, but external mutations + // (host injection, inject_pathway writes) can break it before + // make_consistent. Re-derive here so downstream consumers (calc_tdust_3d.cpp, + // cooling tables) see a consistent bulk field. + if (my_chemistry->dust_species_track == 1) { + const double f_mg_sil = + std::clamp(my_chemistry->dust_silicate_mg_fraction, 0.0, 1.0); + grackle::impl::View dust_bulk( + my_fields->dust_density, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_sil( + my_fields->dust_density_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_mg_sil( + my_fields->dust_density_mg_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_fe_sil( + my_fields->dust_density_fe_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_carb( + const_cast(my_fields->dust_density_carbonaceous), + my_fields->grid_dimension[0], my_fields->grid_dimension[1], + my_fields->grid_dimension[2]); + for (k = my_fields->grid_start[2]; k <= my_fields->grid_end[2]; k++) { + for (j = my_fields->grid_start[1]; j <= my_fields->grid_end[1]; j++) { + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { + double mg_sil = dust_mg_sil(i, j, k); + double fe_sil = dust_fe_sil(i, j, k); + const double sil_compat = dust_sil(i, j, k); + if (mg_sil + fe_sil <= 0.0 && sil_compat > 0.0) { + mg_sil = f_mg_sil * sil_compat; + fe_sil = (1.0 - f_mg_sil) * sil_compat; + dust_mg_sil(i, j, k) = (gr_float)mg_sil; + dust_fe_sil(i, j, k) = (gr_float)fe_sil; + } + const double sil = mg_sil + fe_sil; + dust_sil(i, j, k) = (gr_float)sil; + dust_bulk(i, j, k) = (gr_float)(sil + dust_carb(i, j, k)); + } + } + } + } + return; } diff --git a/src/clib/scale_fields.cpp b/src/clib/scale_fields.cpp index bd5f85b72..a0cb1c120 100644 --- a/src/clib/scale_fields.cpp +++ b/src/clib/scale_fields.cpp @@ -109,6 +109,33 @@ void scale_fields(int imetal, gr_float factor, chemistry_data* my_chemistry, 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 metal_C( + my_fields->metal_density_carbon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View metal_O( + my_fields->metal_density_oxygen, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View metal_Mg( + my_fields->metal_density_magnesium, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View metal_Si( + my_fields->metal_density_silicon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View metal_Fe( + my_fields->metal_density_iron, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_sil( + my_fields->dust_density_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_mg_sil( + my_fields->dust_density_mg_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_fe_sil( + my_fields->dust_density_fe_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + grackle::impl::View dust_carb( + my_fields->dust_density_carbonaceous, 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]); @@ -268,11 +295,21 @@ void scale_fields(int imetal, gr_float factor, chemistry_data* my_chemistry, } } - if (imetal == 1) { + if ((imetal == 1) || (my_chemistry->dust_species_track == 1)) { for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { metal(i, j, k) = metal(i, j, k) * factor; } - + } + if (my_chemistry->dust_species_track == 1) { + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { + metal_C(i, j, k) = metal_C(i, j, k) * factor; + metal_O(i, j, k) = metal_O(i, j, k) * factor; + metal_Mg(i, j, k) = metal_Mg(i, j, k) * factor; + metal_Si(i, j, k) = metal_Si(i, j, k) * factor; + metal_Fe(i, j, k) = metal_Fe(i, j, k) * factor; + } + } + if (imetal == 1) { if (my_chemistry->metal_chemistry == 1) { for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { CI(i, j, k) = CI(i, j, k) * factor; @@ -314,41 +351,50 @@ void scale_fields(int imetal, gr_float factor, chemistry_data* my_chemistry, } } } + } - if (my_chemistry->use_dust_density_field == 1) { + if ((my_chemistry->use_dust_density_field == 1) || + (my_chemistry->dust_species_track == 1)) { + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { + dust(i, j, k) = dust(i, j, k) * factor; + } + if (my_chemistry->dust_species_track == 1) { for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; i++) { - dust(i, j, k) = dust(i, j, k) * factor; + dust_sil(i, j, k) = dust_sil(i, j, k) * factor; + dust_mg_sil(i, j, k) = dust_mg_sil(i, j, k) * factor; + dust_fe_sil(i, j, k) = dust_fe_sil(i, j, k) * factor; + dust_carb(i, j, k) = dust_carb(i, j, k) * factor; } + } - if ((my_chemistry->grain_growth == 1) || - (my_chemistry->dust_sublimation == 1)) { - if (my_chemistry->dust_species > 0) { - for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; - i++) { - MgSiO3(i, j, k) = MgSiO3(i, j, k) * factor; - AC(i, j, k) = AC(i, j, k) * factor; - } + if ((my_chemistry->grain_growth == 1) || + (my_chemistry->dust_sublimation == 1)) { + if (my_chemistry->dust_species > 0) { + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; + i++) { + MgSiO3(i, j, k) = MgSiO3(i, j, k) * factor; + AC(i, j, k) = AC(i, j, k) * factor; } - if (my_chemistry->dust_species > 1) { - for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; - i++) { - SiM(i, j, k) = SiM(i, j, k) * factor; - FeM(i, j, k) = FeM(i, j, k) * factor; - Mg2SiO4(i, j, k) = Mg2SiO4(i, j, k) * factor; - Fe3O4(i, j, k) = Fe3O4(i, j, k) * factor; - SiO2D(i, j, k) = SiO2D(i, j, k) * factor; - MgO(i, j, k) = MgO(i, j, k) * factor; - FeS(i, j, k) = FeS(i, j, k) * factor; - Al2O3(i, j, k) = Al2O3(i, j, k) * factor; - } + } + if (my_chemistry->dust_species > 1) { + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; + i++) { + SiM(i, j, k) = SiM(i, j, k) * factor; + FeM(i, j, k) = FeM(i, j, k) * factor; + Mg2SiO4(i, j, k) = Mg2SiO4(i, j, k) * factor; + Fe3O4(i, j, k) = Fe3O4(i, j, k) * factor; + SiO2D(i, j, k) = SiO2D(i, j, k) * factor; + MgO(i, j, k) = MgO(i, j, k) * factor; + FeS(i, j, k) = FeS(i, j, k) * factor; + Al2O3(i, j, k) = Al2O3(i, j, k) * factor; } - if (my_chemistry->dust_species > 2) { - for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; - i++) { - reforg(i, j, k) = reforg(i, j, k) * factor; - volorg(i, j, k) = volorg(i, j, k) * factor; - H2Oice(i, j, k) = H2Oice(i, j, k) * factor; - } + } + if (my_chemistry->dust_species > 2) { + for (i = my_fields->grid_start[0]; i <= my_fields->grid_end[0]; + i++) { + reforg(i, j, k) = reforg(i, j, k) * factor; + volorg(i, j, k) = volorg(i, j, k) * factor; + H2Oice(i, j, k) = H2Oice(i, j, k) * factor; } } } diff --git a/src/clib/scale_fields.hpp b/src/clib/scale_fields.hpp index d7499d887..c8ea32e0a 100644 --- a/src/clib/scale_fields.hpp +++ b/src/clib/scale_fields.hpp @@ -53,6 +53,114 @@ inline void scale_fields_table(grackle_field_data* my_fields, double factor) { } } } + if (my_fields->metal_density_carbon != nullptr) { + grackle::impl::View metal_C( + my_fields->metal_density_carbon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + metal_C(i, j, k) = metal_C(i, j, k) * factor; + } + } + } + } + if (my_fields->metal_density_oxygen != nullptr) { + grackle::impl::View metal_O( + my_fields->metal_density_oxygen, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + metal_O(i, j, k) = metal_O(i, j, k) * factor; + } + } + } + } + if (my_fields->metal_density_magnesium != nullptr) { + grackle::impl::View metal_Mg( + my_fields->metal_density_magnesium, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + metal_Mg(i, j, k) = metal_Mg(i, j, k) * factor; + } + } + } + } + if (my_fields->metal_density_silicon != nullptr) { + grackle::impl::View metal_Si( + my_fields->metal_density_silicon, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + metal_Si(i, j, k) = metal_Si(i, j, k) * factor; + } + } + } + } + if (my_fields->metal_density_iron != nullptr) { + grackle::impl::View metal_Fe( + my_fields->metal_density_iron, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + metal_Fe(i, j, k) = metal_Fe(i, j, k) * factor; + } + } + } + } + if (my_fields->dust_density_silicate != nullptr) { + grackle::impl::View dust_sil( + my_fields->dust_density_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + dust_sil(i, j, k) = dust_sil(i, j, k) * factor; + } + } + } + } + if (my_fields->dust_density_mg_silicate != nullptr) { + grackle::impl::View dust_mg_sil( + my_fields->dust_density_mg_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + dust_mg_sil(i, j, k) = dust_mg_sil(i, j, k) * factor; + } + } + } + } + if (my_fields->dust_density_fe_silicate != nullptr) { + grackle::impl::View dust_fe_sil( + my_fields->dust_density_fe_silicate, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + dust_fe_sil(i, j, k) = dust_fe_sil(i, j, k) * factor; + } + } + } + } + if (my_fields->dust_density_carbonaceous != nullptr) { + grackle::impl::View dust_carb( + my_fields->dust_density_carbonaceous, my_fields->grid_dimension[0], + my_fields->grid_dimension[1], my_fields->grid_dimension[2]); + for (int k = grid_start[2]; k <= grid_end[2]; k++) { + for (int j = grid_start[1]; j <= grid_end[1]; j++) { + for (int i = grid_start[0]; i <= grid_end[0]; i++) { + dust_carb(i, j, k) = dust_carb(i, j, k) * factor; + } + } + } + } } /// A helper function for scaling the injection pathway metal density fields @@ -65,13 +173,6 @@ void scale_inject_path_metal_densities_(grackle_field_data* my_fields, gr_float factor, int n_inj_path_ptrs); /// Scales fields related to computing dust temperature -/// -/// @param[in] my_chemistry holds a number of configuration parameters -/// @param[inout] my_fields holds the fields that will be updated in-place -/// @param[in] imetal Specifies whether the metal_density was specified -/// @param[in] factor The factor that is multiplied by the fields -/// @param[in] n_inj_path_ptrs The number of pointers tracked by -/// `my_fields->inject_pathway_metal_density` inline void scale_fields_dust(chemistry_data* my_chemistry, grackle_field_data* my_fields, int imetal, gr_float factor, int n_inj_path_ptrs) { @@ -133,6 +234,20 @@ inline void scale_fields_dust(chemistry_data* my_chemistry, if (imetal == 1) { for (int i = idx_range.i_start; i < idx_range.i_stop; i++) { metal(i, j, k) = metal(i, j, k) * factor; + // if (my_chemistry->multi_metals > 0) { + // metal_loc(i, j, k) = metal_loc(i, j, k) * factor; + // metal_C13(i, j, k) = metal_C13(i, j, k) * factor; + // metal_C20(i, j, k) = metal_C20(i, j, k) * factor; + // metal_C25(i, j, k) = metal_C25(i, j, k) * factor; + // metal_C30(i, j, k) = metal_C30(i, j, k) * factor; + // metal_F13(i, j, k) = metal_F13(i, j, k) * factor; + // metal_F15(i, j, k) = metal_F15(i, j, k) * factor; + // metal_F50(i, j, k) = metal_F50(i, j, k) * factor; + // metal_F80(i, j, k) = metal_F80(i, j, k) * factor; + // metal_P170(i, j, k) = metal_P170(i, j, k) * factor; + // metal_P200(i, j, k) = metal_P200(i, j, k) * factor; + // metal_Y19(i, j, k) = metal_Y19(i, j, k) * factor; + // } } } if (my_chemistry->use_dust_density_field == 1) { diff --git a/src/clib/solve_rate_cool.cpp b/src/clib/solve_rate_cool.cpp index 3abafb52c..1bef6f186 100644 --- a/src/clib/solve_rate_cool.cpp +++ b/src/clib/solve_rate_cool.cpp @@ -40,6 +40,7 @@ #include "cool1d_multi_g.hpp" #include "scale_fields.hpp" #include "solve_rate_cool.hpp" +#include "dust/dust_growth_and_destruction.hpp" namespace GRIMPL_NAMESPACE_DECL { @@ -790,6 +791,17 @@ int solve_rate_cool( std::vector nelec_times_mH(my_fields->grid_dimension[0]); std::vector edot(my_fields->grid_dimension[0]); + // Arrays to store dust growth and destruction mass changes + std::vector growth_dM(my_fields->grid_dimension[0]); + std::vector destruction_dM(my_fields->grid_dimension[0]); + // species-specific growth & destruction outputs (dust_species_track == 1) + std::vector growth_dM_mg_silicate(my_fields->grid_dimension[0]); + std::vector growth_dM_fe_silicate(my_fields->grid_dimension[0]); + std::vector growth_dM_carbon(my_fields->grid_dimension[0]); + std::vector destruction_dM_mg_silicate(my_fields->grid_dimension[0]); + std::vector destruction_dM_fe_silicate(my_fields->grid_dimension[0]); + std::vector destruction_dM_carbon(my_fields->grid_dimension[0]); + // iteration masks std::vector itmask(my_fields->grid_dimension[0]); std::vector itmask_metal(my_fields->grid_dimension[0]); @@ -998,6 +1010,46 @@ int solve_rate_cool( ); } + // TEMPORARY: dust growth/destruction is currently invoked here as its + // own block. Eventually, the growth and destruction rates should be + // computed alongside the other dust rates (stored together in the + // newly-created FullRxnRateBuf), and the dust density updates should + // happen alongside the other density updates rather than as a + // separate pass. The placement below is a short-term stopgap and + // will be restructured once that integration lands. + if (my_chemistry->dust_model == 1){ + if (my_chemistry->dust_species_track == 1) { + // Compute and apply Mg-silicate + Fe-silicate + carbonaceous rates + grackle::impl::dust_growth_species( + my_chemistry, my_fields, internalu, idx_range, itmask.data(), + dtit.data(), tgas.data(), + growth_dM_mg_silicate.data(), growth_dM_fe_silicate.data(), + growth_dM_carbon.data()); + grackle::impl::dust_destruction_species( + my_chemistry, my_fields, internalu, idx_range, itmask.data(), + dtit.data(), dt, tgas.data(), + destruction_dM_mg_silicate.data(), destruction_dM_fe_silicate.data(), + destruction_dM_carbon.data()); + grackle::impl::dust_update_species( + my_chemistry, my_fields, internalu, idx_range, itmask.data(), + dtit.data(), + growth_dM_mg_silicate.data(), growth_dM_fe_silicate.data(), + growth_dM_carbon.data(), + destruction_dM_mg_silicate.data(), destruction_dM_fe_silicate.data(), + destruction_dM_carbon.data(), + false); + } else { + grackle::impl::dust_growth( + my_chemistry, my_fields, internalu, idx_range, itmask.data(), + dtit.data(), tgas.data(), growth_dM.data()); + grackle::impl::dust_destruction( + my_chemistry, my_fields, internalu, idx_range, itmask.data(), + dtit.data(), dt, tgas.data(), destruction_dM.data()); + grackle::impl::dust_update( + my_chemistry, my_fields, internalu, idx_range, itmask.data(), dtit.data(), + growth_dM.data(), destruction_dM.data(), false); + } + } // Add the timestep to the elapsed time for each cell and find // minimum elapsed time step in this row @@ -1013,6 +1065,7 @@ int solve_rate_cool( // If all cells are done (in idx_range), break out of subcycle loop if (std::fabs(dt-ttmin) < tolerance*dt) { break; } + } // subcycle iteration loop (for current idx_range) // review number of iterations that were spent in the subcycle loop diff --git a/src/clib/tabulated/calc_temp1d_cloudy.cpp b/src/clib/tabulated/calc_temp1d_cloudy.cpp index 74d0f5786..f471cb36b 100644 --- a/src/clib/tabulated/calc_temp1d_cloudy.cpp +++ b/src/clib/tabulated/calc_temp1d_cloudy.cpp @@ -139,11 +139,10 @@ void calc_temp1d_cloudy(const double* rhoH, double* tgas, double* mmw, // Add metal species to mean molecular weight if (imetal == 1) { + double total_metal_1d = metal(i, idx_range.j, idx_range.k); munew = d(i, idx_range.j, idx_range.k) / - ((d(i, idx_range.j, idx_range.k) - - metal(i, idx_range.j, idx_range.k)) / - munew + - metal(i, idx_range.j, idx_range.k) / mu_metal); + ((d(i, idx_range.j, idx_range.k) - total_metal_1d) / munew + + total_metal_1d / mu_metal); tgas[i] = tgas[i] * munew / muold; } diff --git a/src/include/grackle_chemistry_data.h b/src/include/grackle_chemistry_data.h index 499de6095..13c4b6d0a 100644 --- a/src/include/grackle_chemistry_data.h +++ b/src/include/grackle_chemistry_data.h @@ -312,6 +312,79 @@ typedef struct */ int solver_method; + /* dust model selection + * 0: default (runs calc_grain_size_increment_1d as usual) + * 1: harrison dust model + */ + int dust_model; + + /* Flag to use snetimestep */ + int use_sne_field; + + /* Flag to use user-provided dust destruction timescale field */ + int use_tau_dest_field; + + /* flag and parameters for Li+ 2019 dust growth and destruction */ + double dust_destruction_eff; + double sne_coeff; + double sne_shockspeed; + double dust_grainsize; + double dust_growth_densref; + double dust_growth_tauref; + /* Sticking coefficient S for species-specific dust growth. + dust_growth_species rescales tau_ref by 0.3 / S. */ + double dust_growth_sticking_coeff; + + /* parameters for dust creation*/ + double dust_condensation_eff; + double sne_metal_yield; + + /* Species-resolved dust tracking. + Requires dust_model=1. + 0) off — bulk dust_density (default) + 1) on — evolves dust_density_mg_silicate, dust_density_fe_silicate, + dust_density_carbonaceous and 5-element gas tracking + (C, O, Mg, Si, Fe). dust_density_silicate is maintained as + Mg-silicate + Fe-silicate for compatibility. + REF: Trayford+2026 MNRAS 545, staf2040 */ + int dust_species_track; + + /* COLIBRE-style species growth reference timescales [Myr]. + Normalized at n_H' = 10 cm^-3, T = 10 K, S = 0.3, a = 0.1 micron, + and solar bottleneck abundance. Rescaled by dust_grainsize and + dust_growth_sticking_coeff. REF: Trayford+2026 MNRAS 545, staf2040. */ + double dust_growth_tauref_silicate; + double dust_growth_tauref_carbon; + + /* COLIBRE-style unresolved-density boost for accretion: C rises from 1 to + dust_growth_clumping_factor_max between dust_growth_clumping_nH_min and + dust_growth_clumping_nH_max, using log-linear interpolation in n_H. */ + double dust_growth_clumping_factor_max; + double dust_growth_clumping_nH_min; + double dust_growth_clumping_nH_max; + + /* COLIBRE / Tsai-Mathews thermal sputtering reference time [Myr]. + tau_sp = tau_ref * (a/0.1) * (n_H/1 cm^-3)^-1 + * [1 + (T/2e6 K)^-2.5]. */ + double dust_sputter_tauref; + + /* Initial/fallback mass fraction of silicate dust in Mg-rich silicate. + COLIBRE seeds equal numbers of Mg2SiO4 and Fe2SiO4 molecules, which + corresponds to this mass fraction rather than an equal-mass split. */ + double dust_silicate_mg_fraction; + + /* Mg-rich silicate endmember stoichiometric mass fractions for Mg2SiO4. */ + double dust_mg_silicate_f_Mg; + double dust_mg_silicate_f_Fe; + double dust_mg_silicate_f_Si; + double dust_mg_silicate_f_O; + + /* Fe-rich silicate endmember stoichiometric mass fractions for Fe2SiO4. */ + double dust_fe_silicate_f_Mg; + double dust_fe_silicate_f_Fe; + double dust_fe_silicate_f_Si; + double dust_fe_silicate_f_O; + } chemistry_data; /***************************** diff --git a/src/include/grackle_fortran_interface.def b/src/include/grackle_fortran_interface.def index 7b541d1c3..391b3c7b3 100644 --- a/src/include/grackle_fortran_interface.def +++ b/src/include/grackle_fortran_interface.def @@ -49,7 +49,16 @@ c This is the fortran definition of grackle_field_data TYPE(C_PTR) :: y_velocity TYPE(C_PTR) :: z_velocity TYPE(C_PTR) :: metal_density + TYPE(C_PTR) :: metal_density_carbon + TYPE(C_PTR) :: metal_density_oxygen + TYPE(C_PTR) :: metal_density_magnesium + TYPE(C_PTR) :: metal_density_silicon + TYPE(C_PTR) :: metal_density_iron TYPE(C_PTR) :: dust_density + TYPE(C_PTR) :: dust_density_silicate + TYPE(C_PTR) :: dust_density_mg_silicate + TYPE(C_PTR) :: dust_density_fe_silicate + TYPE(C_PTR) :: dust_density_carbonaceous TYPE(C_PTR) :: e_density TYPE(C_PTR) :: HI_density TYPE(C_PTR) :: HII_density @@ -133,6 +142,8 @@ c This is the fortran definition of grackle_field_data TYPE(C_PTR) :: ref_org_dust_temperature TYPE(C_PTR) :: vol_org_dust_temperature TYPE(C_PTR) :: H2O_ice_dust_temperature + TYPE(C_PTR) :: sne_rate + TYPE(C_PTR) :: tau_dest END TYPE c This is the fortran definition of grackle_chemistry_data @@ -216,7 +227,36 @@ c This is the fortran definition of grackle_chemistry_data INTEGER(C_INT) :: uniform_grain_isrf_heating_rate INTEGER(C_INT) :: max_iterations INTEGER(C_INT) :: exit_after_iterations_exceeded -cc INTEGER(C_INT) :: omp_nthreads // not supported in fortran + INTEGER(C_INT) :: omp_nthreads + INTEGER(C_INT) :: solver_method + INTEGER(C_INT) :: dust_model + INTEGER(C_INT) :: use_sne_field + INTEGER(C_INT) :: use_tau_dest_field + REAL(C_DOUBLE) :: dust_destruction_eff + REAL(C_DOUBLE) :: sne_coeff + REAL(C_DOUBLE) :: sne_shockspeed + REAL(C_DOUBLE) :: dust_grainsize + REAL(C_DOUBLE) :: dust_growth_densref + REAL(C_DOUBLE) :: dust_growth_tauref + REAL(C_DOUBLE) :: dust_growth_sticking_coeff + REAL(C_DOUBLE) :: dust_condensation_eff + REAL(C_DOUBLE) :: sne_metal_yield + INTEGER(C_INT) :: dust_species_track + REAL(C_DOUBLE) :: dust_growth_tauref_silicate + REAL(C_DOUBLE) :: dust_growth_tauref_carbon + REAL(C_DOUBLE) :: dust_growth_clumping_factor_max + REAL(C_DOUBLE) :: dust_growth_clumping_nH_min + REAL(C_DOUBLE) :: dust_growth_clumping_nH_max + REAL(C_DOUBLE) :: dust_sputter_tauref + REAL(C_DOUBLE) :: dust_silicate_mg_fraction + REAL(C_DOUBLE) :: dust_mg_silicate_f_Mg + REAL(C_DOUBLE) :: dust_mg_silicate_f_Fe + REAL(C_DOUBLE) :: dust_mg_silicate_f_Si + REAL(C_DOUBLE) :: dust_mg_silicate_f_O + REAL(C_DOUBLE) :: dust_fe_silicate_f_Mg + REAL(C_DOUBLE) :: dust_fe_silicate_f_Fe + REAL(C_DOUBLE) :: dust_fe_silicate_f_Si + REAL(C_DOUBLE) :: dust_fe_silicate_f_O END TYPE c The following define the fortran interfaces to the C routines diff --git a/src/include/grackle_types.h b/src/include/grackle_types.h index e65e9f81c..fc95a488f 100644 --- a/src/include/grackle_types.h +++ b/src/include/grackle_types.h @@ -72,9 +72,31 @@ typedef struct // metal_cooling = 1 gr_float *metal_density; + // dust_species_track = 1 + // These are *subsets* of metal_density (not separate fields). + gr_float *metal_density_carbon; // gas-phase C (subset of metal_density) + gr_float *metal_density_oxygen; // gas-phase O (subset of metal_density) + + // dust_species_track = 1 + // 5-element gas tracking adds Mg, Si, Fe alongside C, O above. + // Subsets of metal_density. + gr_float *metal_density_magnesium; + gr_float *metal_density_silicon; + gr_float *metal_density_iron; + // use_dust_density_field = 1 gr_float *dust_density; + // dust_species_track = 1 + // Three-species dust: bulk dust_density = Mg-silicate + Fe-silicate + + // carbonaceous. dust_density_silicate is kept as the derived + // Mg-silicate + Fe-silicate sum for compatibility. + // REF: Trayford+2026 MNRAS 545, staf2040. + gr_float *dust_density_silicate; + gr_float *dust_density_mg_silicate; + gr_float *dust_density_fe_silicate; + gr_float *dust_density_carbonaceous; + // primordial_chemistry = 1 gr_float *e_density; gr_float *HI_density; @@ -203,6 +225,11 @@ typedef struct gr_float *vol_org_dust_temperature; gr_float *H2O_ice_dust_temperature; + // use_snetimestep = 1 + gr_float *sne_rate; + + // use_tau_dest_field = 1 + gr_float *tau_dest; } grackle_field_data; diff --git a/src/python/examples/cooling_cell.py b/src/python/examples/cooling_cell.py index a6b9dd045..ba59f0c51 100644 --- a/src/python/examples/cooling_cell.py +++ b/src/python/examples/cooling_cell.py @@ -77,6 +77,8 @@ def main(args=None): my_chemistry.primordial_chemistry = 0 my_chemistry.metal_cooling = 1 my_chemistry.UVbackground = 1 + my_chemistry.dust_model = 1 + my_chemistry.dust_species_track = 0 my_chemistry.grackle_data_file = \ os.path.join(grackle_data_dir, "CloudyData_UVB=HM2012.h5") @@ -121,4 +123,4 @@ def main(args=None): if __name__ == '__main__': - sys.exit(main()) + sys.exit(main()) \ No newline at end of file diff --git a/src/python/examples/freefall.py b/src/python/examples/freefall.py index aabbff225..01507821d 100644 --- a/src/python/examples/freefall.py +++ b/src/python/examples/freefall.py @@ -139,4 +139,4 @@ def main(args=None): if __name__ == "__main__": - sys.exit(main()) + sys.exit(main()) \ No newline at end of file diff --git a/src/python/gracklepy/fluid_container.py b/src/python/gracklepy/fluid_container.py index c1df6f845..6d5d1afac 100644 --- a/src/python/gracklepy/fluid_container.py +++ b/src/python/gracklepy/fluid_container.py @@ -284,6 +284,28 @@ def _required_density_fields(my_chemistry): my_fields.append("metal_density") if my_chemistry.dust_chemistry == 1: my_fields.append("dust_density") + if my_chemistry.dust_species_track == 1: + # Species-resolved dust (Mg-silicate + Fe-silicate + carbonaceous) + # with 5-element gas tracking. dust_density_silicate is a compatibility + # sum of Mg-silicate + Fe-silicate. + # REF: Trayford+2026 MNRAS 545, staf2040. + # Bulk metal_density and dust_density are required invariants for this + # path, even when cooling itself is disabled. + if "metal_density" not in my_fields: + my_fields.append("metal_density") + if "dust_density" not in my_fields: + my_fields.append("dust_density") + my_fields.extend([ + "metal_density_carbon", + "metal_density_oxygen", + "metal_density_magnesium", + "metal_density_silicon", + "metal_density_iron", + "dust_density_silicate", + "dust_density_mg_silicate", + "dust_density_fe_silicate", + "dust_density_carbonaceous", + ]) if my_chemistry.metal_chemistry > 0: my_fields.extend(_dust_metal_densities[my_chemistry.dust_species]) my_fields.extend(_dust_densities[my_chemistry.dust_species]) @@ -319,6 +341,10 @@ def _required_extra_fields(my_chemistry): my_fields.append("H2_custom_shielding_factor") if my_chemistry.use_isrf_field == 1: my_fields.append("isrf_habing") + if my_chemistry.use_sne_field == 1: + my_fields.append("sne_rate") + if my_chemistry.use_tau_dest_field == 1: + my_fields.append("tau_dest") return my_fields def _required_calculated_fields(my_chemistry): diff --git a/src/python/gracklepy/grackle_defs.pxd b/src/python/gracklepy/grackle_defs.pxd index 9b9215366..5fa65edde 100644 --- a/src/python/gracklepy/grackle_defs.pxd +++ b/src/python/gracklepy/grackle_defs.pxd @@ -109,7 +109,16 @@ cdef extern from "grackle.h": gr_float *y_velocity; gr_float *z_velocity; gr_float *metal_density; + gr_float *metal_density_carbon; + gr_float *metal_density_oxygen; + gr_float *metal_density_magnesium; + gr_float *metal_density_silicon; + gr_float *metal_density_iron; gr_float *dust_density; + gr_float *dust_density_silicate; + gr_float *dust_density_mg_silicate; + gr_float *dust_density_fe_silicate; + gr_float *dust_density_carbonaceous; gr_float *e_density; gr_float *HI_density; gr_float *HII_density; @@ -200,6 +209,8 @@ cdef extern from "grackle.h": gr_float *ref_org_dust_temperature; gr_float *vol_org_dust_temperature; gr_float *H2O_ice_dust_temperature; + gr_float *sne_rate; + gr_float *tau_dest; ctypedef struct c_grackle_version "grackle_version": const char* version; diff --git a/src/python/gracklepy/grackle_wrapper.pyx b/src/python/gracklepy/grackle_wrapper.pyx index a1f6a7a78..64e858b94 100644 --- a/src/python/gracklepy/grackle_wrapper.pyx +++ b/src/python/gracklepy/grackle_wrapper.pyx @@ -655,7 +655,16 @@ cdef c_field_data setup_field_data(object fc, int[::1] buf, my_fields.z_velocity = get_field(fc, "z_velocity") my_fields.metal_density = get_field(fc, "metal_density") + my_fields.metal_density_carbon = get_field(fc, "metal_density_carbon") + my_fields.metal_density_oxygen = get_field(fc, "metal_density_oxygen") + my_fields.metal_density_magnesium = get_field(fc, "metal_density_magnesium") + my_fields.metal_density_silicon = get_field(fc, "metal_density_silicon") + my_fields.metal_density_iron = get_field(fc, "metal_density_iron") my_fields.dust_density = get_field(fc, "dust_density") + my_fields.dust_density_silicate = get_field(fc, "dust_density_silicate") + my_fields.dust_density_mg_silicate = get_field(fc, "dust_density_mg_silicate") + my_fields.dust_density_fe_silicate = get_field(fc, "dust_density_fe_silicate") + my_fields.dust_density_carbonaceous = get_field(fc, "dust_density_carbonaceous") my_fields.e_density = get_field(fc, "e_density") my_fields.HI_density = get_field(fc, "HI_density") @@ -764,6 +773,9 @@ cdef c_field_data setup_field_data(object fc, int[::1] buf, my_fields.vol_org_dust_temperature = get_field(fc, "vol_org_dust_temperature") my_fields.H2O_ice_dust_temperature = get_field(fc, "H2O_ice_dust_temperature") + my_fields.sne_rate = get_field(fc, "sne_rate") + my_fields.tau_dest = get_field(fc, "tau_dest") + return my_fields def solve_chemistry(fc, my_dt): diff --git a/src/python/gracklepy/utilities/convenience.py b/src/python/gracklepy/utilities/convenience.py index 2cee3f2f3..1663078ce 100644 --- a/src/python/gracklepy/utilities/convenience.py +++ b/src/python/gracklepy/utilities/convenience.py @@ -20,6 +20,7 @@ FluidContainer from gracklepy.utilities.atomic import \ + atomic_mass, \ atomic_number, \ mass_number, \ primordial_elements, \ @@ -28,6 +29,117 @@ mass_hydrogen_cgs, \ sec_per_Myr +_DUST_SPECIES_ELEMENT_FIELDS = { + "C": "metal_density_carbon", + "O": "metal_density_oxygen", + "Mg": "metal_density_magnesium", + "Si": "metal_density_silicon", + "Fe": "metal_density_iron", +} + +# Default Mg-silicate / Fe-silicate / carbonaceous mass split for IC seeding. +# REF: Draine 2003 ARA&A 41, 241; Zubko, Dwek & Arendt 2004 ApJS 152, 211 — +# canonical MW diffuse-ISM split is ~0.6-0.7 silicate, ~0.3-0.4 carbonaceous. +# The Mg/Fe split follows COLIBRE's equal-number Mg2SiO4/Fe2SiO4 seed, which +# is not an equal-mass split. +_DUST_SPECIES_FRACTIONS = { + "silicate": 0.65, + "mg_silicate": 0.408428, + "fe_silicate": 0.591572, + "carbonaceous": 0.35, +} + + +def solar_metal_mass_fractions(elements=("C", "O", "Mg", "Si", "Fe")): + """ + Mass fractions of given elements relative to total solar metal mass. + + Converts number-density ratios n_X/n_H from atomic.solar_abundance into + mass fractions of the total metal mass via + f_X = (n_X * A_X) / sum_Y (n_Y * A_Y) + where the sum runs over all metals (everything except H, He) so the + fractions sum to <= 1 (the 5 dust-relevant elements account for ~80% + of solar metal mass). + """ + metals = [el for el in solar_abundance if el not in ("H", "He")] + total = sum(solar_abundance[el] * atomic_mass[el] for el in metals) + return {el: solar_abundance[el] * atomic_mass[el] / total for el in elements} + + +def seed_dust_species_dust(fc): + """ + Split fc['dust_density'] into Mg-silicate / Fe-silicate / carbonaceous + reservoirs using canonical MW diffuse-ISM mass fractions (Draine 2003) + and the COLIBRE equal-molecule Mg2SiO4/Fe2SiO4 seed split. + """ + silicate = ( + _DUST_SPECIES_FRACTIONS["silicate"] * fc["dust_density"] + ) + fc["dust_density_mg_silicate"][:] = ( + _DUST_SPECIES_FRACTIONS["mg_silicate"] * silicate + ) + fc["dust_density_fe_silicate"][:] = ( + _DUST_SPECIES_FRACTIONS["fe_silicate"] * silicate + ) + fc["dust_density_silicate"][:] = ( + fc["dust_density_mg_silicate"] + fc["dust_density_fe_silicate"] + ) + fc["dust_density_carbonaceous"][:] = ( + _DUST_SPECIES_FRACTIONS["carbonaceous"] * fc["dust_density"] + ) + + +def seed_dust_species_metal_elements(fc): + """ + Seed the gas-phase per-element reservoirs so that for each tracked + element X in {C, O, Mg, Si, Fe} the total elemental budget is + conserved: + + gas_X + dust_X = solar_X_share * metal_density_total + + where ``metal_density_total`` is the input fc['metal_density'] + (the Z-scaled solar metal share) before any subtraction. Must run + AFTER seed_dust_species_dust(): reads the dust species fields and + the chemistry_data stoichiometric f_X attributes to compute dust_X. + + Also reduces fc['metal_density'] to its gas-phase-only value + (subtracting total dust mass) so the C++ make_consistent invariant + ``metal_density = sum(tracked gas X) + untracked metals`` holds. + + Why: prior implementation set gas_X = solar_X_share * metal_density + while dust was seeded independently. Total elemental budget exceeded + solar by (1 + DGR/metal_mass_fraction) ≈ 1.77x at MW conditions, + putting IC DTM at ~0.55 instead of the physically motivated ~0.15. + """ + chem = fc.chemistry_data + fractions = solar_metal_mass_fractions(_DUST_SPECIES_ELEMENT_FIELDS.keys()) + + metal_density_total = np.asarray(fc["metal_density"]).copy() + dust_mg_sil = np.asarray(fc["dust_density_mg_silicate"]) + dust_fe_sil = np.asarray(fc["dust_density_fe_silicate"]) + dust_carb = np.asarray(fc["dust_density_carbonaceous"]) + + dust_per_element = { + "C": dust_carb, + "O": (chem.dust_mg_silicate_f_O * dust_mg_sil + + chem.dust_fe_silicate_f_O * dust_fe_sil), + "Mg": chem.dust_mg_silicate_f_Mg * dust_mg_sil, + "Si": (chem.dust_mg_silicate_f_Si * dust_mg_sil + + chem.dust_fe_silicate_f_Si * dust_fe_sil), + "Fe": chem.dust_fe_silicate_f_Fe * dust_fe_sil, + } + + total_dust = np.zeros_like(metal_density_total) + for el, field in _DUST_SPECIES_ELEMENT_FIELDS.items(): + total_X = fractions[el] * metal_density_total + fc[field][:] = np.maximum(total_X - dust_per_element[el], 0.0) + total_dust = total_dust + dust_per_element[el] + + fc["metal_density"][:] = np.maximum( + metal_density_total - total_dust, 0.0 + ) + + solar_mass_abundance = {element: solar_abundance[element] * mass_number[element] for element in solar_abundance @@ -313,6 +425,13 @@ def setup_fluid_container(my_chemistry, for field in fc.density_fields: fc[field][:] = state_vals.get(field, tiny_density) + if my_chemistry.dust_species_track == 1: + # Order matters: seed dust species first so that + # seed_dust_species_metal_elements can subtract their mass from + # the per-element gas-phase budgets (elemental conservation). + seed_dust_species_dust(fc) + seed_dust_species_metal_elements(fc) + fc.calculate_mean_molecular_weight() fc["internal_energy"] = temperature / \ fc.chemistry_data.temperature_units / \