/* Local includes. */ #include "dust.h" #include "chemistry.h" #include "cooling.h" #include "dust_constants.h" #include "dust_properties.h" #include "dust_yield_tables.h" #include /** * @brief redistribute any dust mass back to element abundances * on star particle formation according to dust composition, to * represent astration * * @param p The gas particles. * @param dp Global dust parameters for initialisation. */ void dust_redistribute_masses_astration(struct spart *sp, const struct part *p, const struct dustevo_props *dp) { /** * iterate through grain species and element types and redistribute dust * abundances to element abundances, according to composition array */ for (int grain = 0; grain < dust_grain_species_count; grain++) { const float grain_mass = p->dust_data.grain_mass_fraction[grain]; for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { const int eldx = dp->grain_element_indices[grain][elem]; /* put fraction of grain mass back into constituent element (astration) */ if (dp->pair_to_cooling) { sp->chemistry_data.metal_mass_fraction[eldx] += grain_mass * dp->grain_element_mfrac[grain][elem]; } } } } /** * @brief Prints the dust evolution model to stdout. * * @param dust #dustevo_props struct. */ void dustevo_print_backend(const struct dustevo_props *dp) { message("Running with a Trayford et al. (2020) dust evolution model."); } /** * @brief initialise structure housing global dust parametrisation. * In particular, flags and values set in the parameter file, * and any hard-coded properties * @brief Prints the dust evolution model to stdout. * * @param dp Global dust parameters for initialisation * @param params The parsed parameter file. * @param phys_const The physical constants in internal units. * @param us The current internal system of units. * @param dust #dustevo_props struct. */ void dustevo_props_init_backend(struct dustevo_props *dp, struct swift_params *params, struct feedback_props *fp, struct cooling_function_data *cooling, const struct phys_const *phys_const, const struct unit_system *us) { /** * initialise structure housing global dust parametrisation. * First, set global parameters: * * initialise dustevo_props struct object and store: * - Parameter file flags (e.g. with_accretion, with_sputtering) * - Parameter file values (e.g diffusion boost, clumping factor) * - Hard coded values (e.g. grain radii, density, composition) * * * Then, using the yield tables from feedback_props, compute * and set dust yields for each of the dust species * * Finally apply corections to the cooling and heating rates * to remove implicit dust, via the scale_out_table_depletion * function. **/ /* read some parameters */ /* operation modes */ dp->pair_to_cooling = parser_get_opt_param_int(params, "DustEvolution:pair_to_cooling", 0); dp->with_subgrid_props = parser_get_opt_param_int(params, "DustEvolution:use_subgrid_props", 1); dp->with_sputtering = parser_get_opt_param_int(params, "DustEvolution:use_sputtering", 1); dp->with_SNII_destruction = parser_get_opt_param_int(params, "DustEvolution:use_SNII_destruction", 1); dp->with_accretion = parser_get_opt_param_int(params, "DustEvolution:use_accretion", 1); /* initial abundances */ dp->initial_grain_mass_fraction[dust_grain_species_graphite] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_graphite", 0.); dp->initial_grain_mass_fraction[dust_grain_species_silicate] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_silicate", 0.); /* global parameters */ dp->clumping_factor = parser_get_opt_param_float( params, "DustEvolution:clumping_factor", default_clumping_factor); dp->diffusion_rate_boost = parser_get_opt_param_float( params, "DustEvolution:diffusion_boost_factor", 1.0); dp->nu = parser_get_opt_param_float( params, "DustEvolution:silicate_fe_grain_fraction", default_fe_grain_fraction); parser_get_opt_param_string(params, "DustEvolution:dust_yields_path", dp->AGB_dyield_path, "./dust_yields"); /* initialise 1D arrays */ memset(dp->abundance_pattern, 0, sizeof(dp->abundance_pattern)); memset(dp->atomic_weight, 0, sizeof(dp->atomic_weight)); memset(dp->condensation_frac, 0, sizeof(dp->condensation_frac)); /* set assumed abundance patterns to Wiersma et al (2009a) */ dp->abundance_pattern[chemistry_element_H] = solar_abundance_H; dp->abundance_pattern[chemistry_element_He] = solar_abundance_He; dp->abundance_pattern[chemistry_element_C] = solar_abundance_C; dp->abundance_pattern[chemistry_element_N] = solar_abundance_N; dp->abundance_pattern[chemistry_element_O] = solar_abundance_O; dp->abundance_pattern[chemistry_element_Ne] = solar_abundance_Ne; dp->abundance_pattern[chemistry_element_Mg] = solar_abundance_Mg; dp->abundance_pattern[chemistry_element_Si] = solar_abundance_Si; dp->abundance_pattern[chemistry_element_Fe] = solar_abundance_Fe; dp->solar_metallicity = solar_metallicity_wiersma; /* set element atomic weights */ dp->atomic_weight[chemistry_element_H] = atomic_weight_H; dp->atomic_weight[chemistry_element_He] = atomic_weight_He; dp->atomic_weight[chemistry_element_C] = atomic_weight_C; dp->atomic_weight[chemistry_element_N] = atomic_weight_N; dp->atomic_weight[chemistry_element_O] = atomic_weight_O; dp->atomic_weight[chemistry_element_Ne] = atomic_weight_Ne; dp->atomic_weight[chemistry_element_Mg] = atomic_weight_Mg; dp->atomic_weight[chemistry_element_Si] = atomic_weight_Si; dp->atomic_weight[chemistry_element_Fe] = atomic_weight_Fe; /* set element count contributing to each grain */ dp->grain_element_count[dust_grain_species_graphite] = atomicity_graphite; dp->grain_element_count[dust_grain_species_silicate] = atomicity_silicate; /* set condensation fractions */ dp->condensation_frac[dust_grain_species_graphite] = graphite_condensation_efficiency; dp->condensation_frac[dust_grain_species_silicate] = silicate_condensation_efficiency; /* set element count contributing to each grain */ dp->accretion_coeff[dust_grain_species_graphite] = graphite_accretion_timescale_norm; dp->accretion_coeff[dust_grain_species_silicate] = silicate_accretion_timescale_norm; /* number of SNII per unit stellar mass */ dp->specific_numSNII = SNII_mass_fraction / phys_const->const_solar_mass; /* set 2d grain arrays */ dustevo_set_jagged_arrays(dp); /** NOTE 1: total metallicity yields untouched here, so Z represents the * conserved dust + gas phase metals **/ /** NOTE 2: only the IMF resampled tables are modified in fp, while plain * .yield arrays are unchanged (the original yield tables are only used to * compute these and are already modified via the SNII yield factors) **/ initialise_dyield_tables(fp, dp); compute_AGB_dyield(fp, dp); compute_SNII_dyield(fp, dp); // print_dyield_tables(fp, dp); } /** * @brief Write dustevo props to the given FILE as a stream of bytes. * * @param dp the struct * @param stream the file stream */ void dustevo_props_dump(const struct dustevo_props *dp, FILE *stream) { /* Zero all pointers in the dustevo_props struct. */ struct dustevo_props dp_copy = *dp; dp_copy.logfD = NULL; dp_copy.dyield_AGB.yield = NULL; dp_copy.dyield_AGB.yield_IMF_resampled = NULL; dp_copy.dyield_SNII.yield = NULL; dp_copy.dyield_SNII.yield_IMF_resampled = NULL; /* iterate through pointer arrays and zero */ for (int grain = 0; grain < dust_grain_species_count; grain++) { dp_copy.grain_element_indices[grain] = NULL; dp_copy.grain_element_mfrac[grain] = NULL; } restart_write_blocks((void *)&dp_copy, sizeof(struct dustevo_props), 1, stream, "dp", "dustevo properties"); } /** * @brief Restore a dustevo_props struct from the given FILE as a stream of * bytes. * * @param dp the struct * @param fp #feedback_props structure * @param stream the file stream */ void dustevo_props_restore(struct dustevo_props *dp, struct feedback_props *fp, FILE *stream) { restart_read_blocks((void *)dp, sizeof(struct dustevo_props), 1, stream, NULL, "dustevo properties"); dustevo_set_jagged_arrays(dp); initialise_dyield_tables(fp, dp); compute_AGB_dyield(fp, dp); compute_SNII_dyield(fp, dp); // print_dyield_tables(fp, dp); } /** * @brief Set jagged arrays in the dustevo_props structure * * @param dp the struct */ void dustevo_set_jagged_arrays(struct dustevo_props *dp) { /* effective grain molecular weights in amu */ float Agra = dp->atomic_weight[chemistry_element_C]; float Asil = 2 * dp->nu * dp->atomic_weight[chemistry_element_Fe] + (2 - 2 * dp->nu) * dp->atomic_weight[chemistry_element_Mg] + dp->atomic_weight[chemistry_element_Si] + 4 * dp->atomic_weight[chemistry_element_O]; /* set grain composition */ float gcomp1[1] = {dp->atomic_weight[chemistry_element_C] / Agra}; float gcomp2[4] = { 4 * dp->atomic_weight[chemistry_element_O] / Asil, (2 - 2 * dp->nu) * dp->atomic_weight[chemistry_element_Mg] / Asil, dp->atomic_weight[chemistry_element_Si] / Asil, 2 * dp->nu * dp->atomic_weight[chemistry_element_Fe] / Asil}; dp->grain_element_mfrac[dust_grain_species_graphite] = gcomp1; dp->grain_element_mfrac[dust_grain_species_silicate] = gcomp2; /* allocate memory */ dp->grain_element_mfrac[dust_grain_species_graphite] = malloc( dp->grain_element_count[dust_grain_species_graphite] * sizeof(float)); dp->grain_element_mfrac[dust_grain_species_silicate] = malloc( dp->grain_element_count[dust_grain_species_silicate] * sizeof(float)); /* deep copy grain arrays */ memcpy(dp->grain_element_mfrac[dust_grain_species_graphite], &gcomp1, dp->grain_element_count[dust_grain_species_graphite] * sizeof(float)); memcpy(dp->grain_element_mfrac[dust_grain_species_silicate], &gcomp2, dp->grain_element_count[dust_grain_species_silicate] * sizeof(float)); /* set chemistry indices composition */ int gidx1[1] = {chemistry_element_C}; int gidx2[4] = {chemistry_element_O, chemistry_element_Mg, chemistry_element_Si, chemistry_element_Fe}; /* allocate memory */ dp->grain_element_indices[dust_grain_species_graphite] = malloc( dp->grain_element_count[dust_grain_species_graphite] * sizeof(int)); dp->grain_element_indices[dust_grain_species_silicate] = malloc( dp->grain_element_count[dust_grain_species_silicate] * sizeof(int)); /* deep copy grain arrays */ memcpy(dp->grain_element_indices[dust_grain_species_graphite], &gidx1, dp->grain_element_count[dust_grain_species_graphite] * sizeof(int)); memcpy(dp->grain_element_indices[dust_grain_species_silicate], &gidx2, dp->grain_element_count[dust_grain_species_silicate] * sizeof(int)); } void dust_evolve_part(const struct phys_const *phys_const, const struct unit_system *us, const struct cosmology *cosmo, const struct hydro_props *hydro_properties, const struct entropy_floor_properties *floor_props, const struct cooling_function_data *cooling, const struct dustevo_props *dp, struct part *restrict p, struct xpart *restrict xp, const float dt, const float dt_therm, const double time) { const float X = p->chemistry_data.metal_mass_fraction[chemistry_element_H]; const float Z = p->chemistry_data.metal_mass_fraction_total; int eldx; /* metal mass fraction in dust before evolution step */ float D_pre = 0.; for (int grain = 0; grain < dust_grain_species_count; grain++) { D_pre += p->dust_data.grain_mass_fraction[grain]; } /* no dust, no accretion or destruction. return. */ if (D_pre <= 0.) return; /* set temperature and density corresponding to either subgrid or hydro */ float rho, temp; if (dp->with_subgrid_props) { rho = p->cooling_data.subgrid_dens; temp = p->cooling_data.subgrid_temp; } else { rho = hydro_get_physical_density(p, cosmo); temp = cooling_get_temperature(phys_const, hydro_properties, us, cosmo, cooling, p, xp); } /* --------------- Common Parameters --------------- */ /* grain radius in cgs, hard coded for now, 0.1 micron. */ const double grain_radius_cgs = representative_grain_size; /* convert Hydrogen mass fraction into Hydrogen number density */ const double n_H = rho * X / phys_const->const_proton_mass; /* convert some properties to cgs */ const double n_H_cgs = n_H * units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY); const double dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME); /*------------------------- Sputtering -------------------------*/ double deltf_sput; if (dp->with_sputtering) { /* compute fractional change in grain species mass this timestep */ const double dfdr_cgs = 3. / grain_radius_cgs; const double drdt_grain_cgs = -sputtering_timescale_norm * n_H_cgs / (pow(2e6 / temp, 2.5) + 1); const double dfdt_cgs = drdt_grain_cgs * dfdr_cgs; deltf_sput = dfdt_cgs * dt_cgs; } else deltf_sput = 0.; /*------------------------- SNII shocks -------------------------*/ double deltf_SNII; if (dp->with_SNII_destruction && (xp->sf_data.SFR > 0.)) { const float eps_eff = 0.1; const double unit_mass_cgs = units_cgs_conversion_factor(us, UNIT_CONV_MASS); const double unit_time_cgs = units_cgs_conversion_factor(us, UNIT_CONV_TIME); /* SNII rate assuming constant SFR, in cgs units */ const float rate_SNII = xp->sf_data.SFR * dp->specific_numSNII / unit_time_cgs; /* Yamasawa et. al (2011): the cgs mass of material swept up in SNII shocks */ const float m_swept = snii_destruction_coeff * pow(n_H_cgs, -0.202) * pow((Z / dp->solar_metallicity) + 0.039, -0.298); /* destruction timescale of dust in cgs */ const float tau_dest = p->mass * unit_mass_cgs / (eps_eff * m_swept * rate_SNII); /* fractional change in dust due to dust destruction */ deltf_SNII = -dt_cgs / tau_dest; } else { deltf_SNII = 0.; } /* ------------------------- Accretion ------------------------- */ double deltf_acc[dust_grain_species_count] = {0.}; if (dp->with_accretion) { /* accretion timescale multiplicative terms (Hirashita & Voshchinnikov 2013), implicitly assuming S_acc = 0.3 */ const float grain_radius_term = grain_radius_cgs / 1e-5; const float n_H_term = 1.e3 / n_H_cgs; const float T_term = sqrt(10. / temp); const float mult_terms = grain_radius_term * n_H_term * T_term; float min_abundance_ratio = FLT_MAX; if (dp->pair_to_cooling) { for (int grain = 0; grain < dust_grain_species_count; grain++) { min_abundance_ratio = FLT_MAX; for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { eldx = dp->grain_element_indices[grain][elem]; /* find the minimum solar abundance ratio among constituent elements indicating element that bottlenecks accretion of effective grain */ min_abundance_ratio = min( min_abundance_ratio, p->chemistry_data.metal_mass_fraction[eldx] / dp->abundance_pattern[eldx]); } /* i.e. tau_acc = dp->accretion_coeff[grain] * mult_terms / * min_abundance_ratio */ deltf_acc[grain] = dt_cgs * dp->clumping_factor * min_abundance_ratio / (dp->accretion_coeff[grain] * mult_terms); } } else { float diffuse_mfrac; /* unpaired mode - need to find diffuse metal mass fraction*/ for (int grain = 0; grain < dust_grain_species_count; grain++) { min_abundance_ratio = 1e30; // FLT_MAX; for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { eldx = dp->grain_element_indices[grain][elem]; /* find the minimum solar abundance ratio among constituent elements indicating element that bottlenecks accretion of effective grain */ diffuse_mfrac = p->chemistry_data.metal_mass_fraction[eldx] - (dp->grain_element_mfrac[grain][elem] * p->dust_data.grain_mass_fraction[grain]); min_abundance_ratio = min( min_abundance_ratio, diffuse_mfrac / dp->abundance_pattern[eldx]); } /* i.e. tau_acc = dp->accretion_coeff[grain] * mult_terms / * min_abundance_ratio */ deltf_acc[grain] = dt_cgs * dp->clumping_factor * min_abundance_ratio / (dp->accretion_coeff[grain] * mult_terms); } } } /* ------------------------- Combine ------------------------- */ float delta_dfrac; float new_grain_mfrac; float D_post = 0.; for (int grain = 0; grain < dust_grain_species_count; grain++) { delta_dfrac = p->dust_data.grain_mass_fraction[grain] * (exp(deltf_sput + deltf_SNII + deltf_acc[grain]) - 1.); /* limit grain growth */ if (dp->pair_to_cooling) { for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { eldx = dp->grain_element_indices[grain][elem]; /* grain growth limited by available mass in a consitituent element */ delta_dfrac = min(delta_dfrac, p->chemistry_data.metal_mass_fraction[eldx] / dp->grain_element_mfrac[grain][elem]); } } else { float diffuse_mfrac; for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { eldx = dp->grain_element_indices[grain][elem]; /* grain growth limited by available mass in a consitituent element */ diffuse_mfrac = p->chemistry_data.metal_mass_fraction[eldx] - (dp->grain_element_mfrac[grain][elem] * p->dust_data.grain_mass_fraction[grain]); delta_dfrac = min(delta_dfrac, diffuse_mfrac / dp->grain_element_mfrac[grain][elem]); } } new_grain_mfrac = p->dust_data.grain_mass_fraction[grain] + delta_dfrac; if (new_grain_mfrac > FLT_MIN) { /* update particle grain mass fractions */ p->dust_data.grain_mass_fraction[grain] += delta_dfrac; if (dp->pair_to_cooling) { for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { eldx = dp->grain_element_indices[grain][elem]; /* update constituent element abundances in paired mode */ p->chemistry_data.metal_mass_fraction[eldx] -= delta_dfrac * dp->grain_element_mfrac[grain][elem]; } } D_post += p->dust_data.grain_mass_fraction[grain]; } /* limit dust mass fractions positive to mitigate very short destruction timescales */ else { p->dust_data.grain_mass_fraction[grain] = FLT_MIN; if (dp->pair_to_cooling) { for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { eldx = dp->grain_element_indices[grain][elem]; /* update constituent element abundances in paired mode */ p->chemistry_data.metal_mass_fraction[eldx] -= new_grain_mfrac * dp->grain_element_mfrac[grain][elem]; } } D_post += FLT_MIN; } } if (0) { message("delta dust: %e, T %e rho %e", (D_post - D_pre) / D_pre, temp, rho); } }