/* 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 "error.h" #include "star_formation.h" #include float dust_get_local_dust_to_gas_ratio(void) { return local_dust_to_gas_ratio; } /** * @brief Return the ratio of grain sizes assumed in the model */ double dust_get_large_to_small_ratio(void) { return radius_scaling_large / radius_scaling_small; } /** * @brief recompute composition of dust species, particularly due to the * balance of Magnesium to Iron grains. * * @param *comp_array Array to store computed constituent mass fractions. * @param mg_specfrac Silicate mass fraction in Magnesium grains. * @param fe_specfrac Silicate mass fraction in Iron grains. * @param dp Global dust properties. */ void dust_recompute_composition(double *comp_array, const double mg_specfrac, const double fe_specfrac, const struct dustevo_props *dp) { double spec_mfracs[dust_grain_molecule_count] = {1., mg_specfrac, fe_specfrac}; double element_mfrac[chemistry_element_count] = {0.}; for (int mol = 0; mol < dust_grain_molecule_count; mol++) { for (int elem = 0; elem < dp->grain_element_count[mol]; elem++) { const int eldx = dp->grain_element_indices[mol][elem]; element_mfrac[eldx] += spec_mfracs[mol] * dp->grain_element_mfrac[mol][elem]; } } for (int elem = 0; elem < dp->cgrain_element_count[dust_cgrain_species_silicate_l]; elem++) { const int eldx = dp->cgrain_element_indices[dust_cgrain_species_silicate_l][elem]; comp_array[eldx] = element_mfrac[eldx]; } comp_array[chemistry_element_C] = 1.; } /** * @brief redistribute any dust mass back to element abundances * for use on feedback heated gas in the cooling step, and re-setting * particle dust properties to 0. * * @param p The gas particles. * @param dp Global dust parameters for initialisation. */ void dust_redistribute_masses_cooling_step(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 */ if (dp->pair_to_cooling) { 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 */ p->chemistry_data.metal_mass_fraction[eldx] += grain_mass * dp->grain_element_mfrac[grain][elem]; } /* set dust abundances to 0 */ p->dust_data.grain_mass_fraction[grain] = 0; } /* also add back depleted SNIa Iron to total SNIa Iron fraction */ p->chemistry_data.iron_mass_fraction_from_SNIa = min(p->dust_data.grain_iron_mass_fraction_from_SNIa + p->chemistry_data.iron_mass_fraction_from_SNIa, p->chemistry_data.metal_mass_fraction[chemistry_element_Fe]); /* set dust tracer to 0. */ p->dust_data.grain_iron_mass_fraction_from_SNIa = 0.; } } /** * @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 */ if (dp->pair_to_cooling) { 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) */ sp->chemistry_data.metal_mass_fraction[eldx] += grain_mass * dp->grain_element_mfrac[grain][elem]; } } /* also add back depleted SNIa Iron to total SNIa Iron fraction */ sp->chemistry_data.iron_mass_fraction_from_SNIa = min(p->dust_data.grain_iron_mass_fraction_from_SNIa + sp->chemistry_data.iron_mass_fraction_from_SNIa, sp->chemistry_data.metal_mass_fraction[chemistry_element_Fe]); } } /** * @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 (multi grain " "size)."); } /** * @brief Prints the dust evolution model to stdout. * * @param depletion_value the reference depletion factor for a given element * @param reduction_factor the reduction factor to the leftover diffuse-phase * mass fraction */ float dustevo_new_depfac(float depletion_value, float reduction_factor) { return (1. - depletion_value) * reduction_factor; } /** * @brief checks validity of initial abundances for dust given chemistry * * @param dp Global dust parameters for initialisation * @param params The parsed parameter file. */ void dustevo_check_init_abundance(struct dustevo_props *dp, struct swift_params *params) { float initial_depleted_mfrac[chemistry_element_count] = {0.}; for (int grain = 0; grain < dust_grain_species_count; grain++) { const float grain_mass = dp->initial_grain_mass_fraction[grain]; for (int elem = 0; elem < dp->grain_element_count[grain]; elem++) { const int eldx = dp->grain_element_indices[grain][elem]; initial_depleted_mfrac[eldx] += grain_mass * dp->grain_element_mfrac[grain][elem]; } } float init_element_abundance = 0.; for (int elem = 0; elem < chemistry_element_count; ++elem) { char buffer[50]; sprintf(buffer, "COLIBREChemistry:init_abundance_%s", chemistry_get_element_name((enum chemistry_element)elem)); init_element_abundance = parser_get_param_float(params, buffer); if ((init_element_abundance - initial_depleted_mfrac[elem]) < 0.) { error( "Initial abundance error: dust abundance requires more %s than " "total avaialable (would require %f but only %f available)", chemistry_get_element_name((enum chemistry_element)elem), initial_depleted_mfrac[elem], init_element_abundance); } } } /** * @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_param_int(params, "DustEvolution:use_subgrid_props"); dp->with_sputtering = parser_get_param_int(params, "DustEvolution:use_sputtering"); dp->with_SNII_destruction = parser_get_param_int(params, "DustEvolution:use_SNII_destruction"); dp->with_accretion = parser_get_param_int(params, "DustEvolution:use_accretion"); dp->with_coagulation = parser_get_param_int(params, "DustEvolution:use_coagulation"); dp->with_shattering = parser_get_param_int(params, "DustEvolution:use_shattering"); char temp[40]; parser_get_opt_param_string(params, "DustEvolution:clumping_factor_mode", temp, "chimes_synced"); if (strcmp(temp, "constant") == 0) dp->clumping_factor_mode = clumping_mode_constant; else if (strcmp(temp, "chimes_synced") == 0) dp->clumping_factor_mode = clumping_mode_syncedchimes; else if (strcmp(temp, "variable") == 0) dp->clumping_factor_mode = clumping_mode_variable; else error( "Dust Error: Invalid clumping factor mode, can choose from \n " "- 'constant' [constant clumping factor] \n - 'chimes_synced' " " [synchronised clumping with cooling dust boost] \n - 'variable' " "[compute density dependent clumping] \n "); /* initial abundances */ dp->initial_grain_mass_fraction[dust_grain_species_graphite_l] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_graphite_l", 0.); dp->initial_grain_mass_fraction[dust_grain_species_mgsilicate_l] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_mgsilicate_l", 0.); dp->initial_grain_mass_fraction[dust_grain_species_fesilicate_l] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_fesilicate_l", 0.); dp->initial_grain_mass_fraction[dust_grain_species_graphite_s] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_graphite_s", 0.); dp->initial_grain_mass_fraction[dust_grain_species_mgsilicate_s] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_mgsilicate_s", 0.); dp->initial_grain_mass_fraction[dust_grain_species_fesilicate_s] = parser_get_opt_param_float( params, "DustEvolution:initial_abundance_fesilicate_s", 0.); /* Grain size mass fractions at injection */ dp->grain_size_mass_fraction[dust_grain_species_graphite_l] = parser_get_opt_param_float(params, "DustEvolution:graphite_large_grain_mfrac", default_injection_fraction_large); dp->grain_size_mass_fraction[dust_grain_species_mgsilicate_l] = parser_get_opt_param_float(params, "DustEvolution:mgsilicate_large_grain_mfrac", default_injection_fraction_large); dp->grain_size_mass_fraction[dust_grain_species_fesilicate_l] = parser_get_opt_param_float(params, "DustEvolution:fesilicate_large_grain_mfrac", default_injection_fraction_large); dp->grain_size_mass_fraction[dust_grain_species_graphite_s] = parser_get_opt_param_float(params, "DustEvolution:graphite_small_grain_mfrac", default_injection_fraction_small); dp->grain_size_mass_fraction[dust_grain_species_mgsilicate_s] = parser_get_opt_param_float(params, "DustEvolution:mgsilicate_small_grain_mfrac", default_injection_fraction_small); dp->grain_size_mass_fraction[dust_grain_species_fesilicate_s] = parser_get_opt_param_float(params, "DustEvolution:fesilicate_small_grain_mfrac", default_injection_fraction_small); /* 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->coagulation_timescale_norm_prefactor = parser_get_opt_param_float( params, "DustEvolution:coagulation_timescale_norm_prefactor", 1.0); dp->nu = parser_get_opt_param_float( params, "DustEvolution:silicate_fe_grain_fraction", -1.); if (dp->clumping_factor_mode == clumping_mode_variable) { dp->clumping_factor_nH_min_cgs = parser_get_opt_param_float( params, "DustEvolution:clumping_factor_nH_min_cgs", 1.); dp->clumping_factor_nH_max_cgs = parser_get_opt_param_float( params, "DustEvolution:clumping_factor_nH_max_cgs", 10.); } /* raise errors for some forbidden parametrisations */ if (dp->nu != -1.) { error( "Dust Error: silicate_fe_grain_fraction is set but is a " "deprecated parameter and no longer has any effect on " "the dust model. Instead, use silicate_molecule_subscript_X " "parameters to modify silicate composition.\n"); } if ((dp->clumping_factor_mode < 0) || (dp->clumping_factor_mode > 2)) { } #if !defined(COOLING_CHIMES) if (dp->clumping_factor_mode == clumping_mode_syncedchimes) { error( "Dust Error: Chose dust clumping factor mode 'constant' (synchronise " "clumping with cooling dust boost) but this is only valid for " "CHIMES cooling. Either recompile with this cooling type " "or choose a different dust clumping factor mode."); } #endif dp->o_sub = parser_get_opt_param_float( params, "DustEvolution:silicate_molecule_subscript_oxygen", default_subscript_silicate_O); dp->mg_sub = parser_get_opt_param_float( params, "DustEvolution:silicate_molecule_subscript_magnesium", default_subscript_silicate_Mg); dp->si_sub = parser_get_opt_param_float( params, "DustEvolution:silicate_molecule_subscript_silicon", default_subscript_silicate_Si); dp->fe_sub = parser_get_opt_param_float( params, "DustEvolution:silicate_molecule_subscript_iron", default_subscript_silicate_Fe); dp->limit_depletion = parser_get_opt_param_float(params, "DustEvolution:limit_depletion", 1); dp->diff_redfactor = parser_get_opt_param_float( params, "DustEvolution:undepleted_min_reduction_factor", default_diffuse_fraction_scaling); 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 Mg to Fe mass ratio */ dp->mg_to_fe_aweight = atomic_weight_Fe / atomic_weight_Mg; /* initialise max depletions to 100%, then limit depletion if requested */ for (int el = 0; el < chemistry_element_count; el++) { dp->grain_element_min_diffuse[el] = 1.; } if (dp->limit_depletion) { /* Set depletion limits for each grain constituent element assuming * Ploeckinger+20 patterns */ dp->grain_element_min_diffuse[chemistry_element_C] = dustevo_new_depfac(fiducial_diffuse_C, dp->diff_redfactor); dp->grain_element_min_diffuse[chemistry_element_O] = dustevo_new_depfac(fiducial_diffuse_O, dp->diff_redfactor); dp->grain_element_min_diffuse[chemistry_element_Mg] = dustevo_new_depfac(fiducial_diffuse_Mg, dp->diff_redfactor); dp->grain_element_min_diffuse[chemistry_element_Si] = dustevo_new_depfac(fiducial_diffuse_Si, dp->diff_redfactor); dp->grain_element_min_diffuse[chemistry_element_Fe] = dustevo_new_depfac(fiducial_diffuse_Fe, dp->diff_redfactor); } /* set element count contributing to each grain */ dp->grain_element_count[dust_grain_species_graphite_l] = atomicity_graphite; dp->grain_element_count[dust_grain_species_mgsilicate_l] = atomicity_mgsilicate; dp->grain_element_count[dust_grain_species_fesilicate_l] = atomicity_fesilicate; dp->grain_element_count[dust_grain_species_graphite_s] = atomicity_graphite; dp->grain_element_count[dust_grain_species_mgsilicate_s] = atomicity_mgsilicate; dp->grain_element_count[dust_grain_species_fesilicate_s] = atomicity_fesilicate; /* set element count contributing to each composite grain */ dp->cgrain_element_count[dust_cgrain_species_graphite_l] = atomicity_graphite; dp->cgrain_element_count[dust_cgrain_species_silicate_l] = atomicity_silicate; dp->cgrain_element_count[dust_cgrain_species_graphite_s] = atomicity_graphite; dp->cgrain_element_count[dust_cgrain_species_silicate_s] = atomicity_silicate; /* set radius scale factor relative to representative_grain_size */ dp->grain_size_factors[dust_grain_species_graphite_l] = radius_scaling_large; dp->grain_size_factors[dust_grain_species_mgsilicate_l] = radius_scaling_large; dp->grain_size_factors[dust_grain_species_fesilicate_l] = radius_scaling_large; dp->grain_size_factors[dust_grain_species_graphite_s] = radius_scaling_small; dp->grain_size_factors[dust_grain_species_mgsilicate_s] = radius_scaling_small; dp->grain_size_factors[dust_grain_species_fesilicate_s] = radius_scaling_small; /* set condensation fractions */ dp->condensation_frac[dust_cgrain_species_graphite_l] = graphite_condensation_efficiency; dp->condensation_frac[dust_cgrain_species_silicate_l] = silicate_condensation_efficiency; dp->condensation_frac[dust_cgrain_species_graphite_s] = graphite_condensation_efficiency; dp->condensation_frac[dust_cgrain_species_silicate_s] = silicate_condensation_efficiency; /* set element count contributing to each grain */ dp->accretion_coeff[dust_cgrain_species_graphite_l] = graphite_accretion_timescale_norm; dp->accretion_coeff[dust_cgrain_species_silicate_l] = silicate_accretion_timescale_norm; dp->accretion_coeff[dust_cgrain_species_graphite_s] = graphite_accretion_timescale_norm; dp->accretion_coeff[dust_cgrain_species_silicate_s] = silicate_accretion_timescale_norm; /* set indices between composite and subspecies grain arrays */ dp->grain_to_cgrain_indices[dust_grain_species_graphite_l] = dust_cgrain_species_graphite_l; dp->grain_to_cgrain_indices[dust_grain_species_mgsilicate_l] = dust_cgrain_species_silicate_l; dp->grain_to_cgrain_indices[dust_grain_species_fesilicate_l] = dust_cgrain_species_silicate_l; dp->grain_to_cgrain_indices[dust_grain_species_graphite_s] = dust_cgrain_species_graphite_s; dp->grain_to_cgrain_indices[dust_grain_species_mgsilicate_s] = dust_cgrain_species_silicate_s; dp->grain_to_cgrain_indices[dust_grain_species_fesilicate_s] = dust_cgrain_species_silicate_s; /* compute subspecies constituent mass frations */ double init_mgsubspec_mfrac = ((2.0 * dp->atomic_weight[chemistry_element_Fe]) - (dp->mg_sub * dp->atomic_weight[chemistry_element_Mg] + dp->fe_sub * dp->atomic_weight[chemistry_element_Fe])) / (2.0 * (dp->atomic_weight[chemistry_element_Fe] - dp->atomic_weight[chemistry_element_Mg])); dp->grain_species_mass_frac[dust_grain_species_graphite_l] = 1.; dp->grain_species_mass_frac[dust_grain_species_mgsilicate_l] = init_mgsubspec_mfrac; dp->grain_species_mass_frac[dust_grain_species_fesilicate_l] = 1. - init_mgsubspec_mfrac; dp->grain_species_mass_frac[dust_grain_species_graphite_s] = 1.; dp->grain_species_mass_frac[dust_grain_species_mgsilicate_s] = init_mgsubspec_mfrac; dp->grain_species_mass_frac[dust_grain_species_fesilicate_s] = 1. - init_mgsubspec_mfrac; /* 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_grain_arrays(dp); dustevo_set_jagged_cgrain_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) **/ /* check validity of dust abundance initialisation */ dustevo_check_init_abundance(dp, params); 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_grain_arrays(dp); dustevo_set_jagged_cgrain_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 grain arrays in the dustevo_props structure * * @param dp the struct */ void dustevo_set_jagged_grain_arrays(struct dustevo_props *dp) { /* effective grain molecular weights in amu */ double Agra = dp->atomic_weight[chemistry_element_C]; double Afesil = 2 * dp->atomic_weight[chemistry_element_Fe] + dp->si_sub * dp->atomic_weight[chemistry_element_Si] + dp->o_sub * dp->atomic_weight[chemistry_element_O]; double Amgsil = 2 * dp->atomic_weight[chemistry_element_Mg] + dp->si_sub * dp->atomic_weight[chemistry_element_Si] + dp->o_sub * dp->atomic_weight[chemistry_element_O]; /* set grain composition */ double gcomp1[1] = {dp->atomic_weight[chemistry_element_C] / Agra}; double gcomp2[3] = { dp->o_sub * dp->atomic_weight[chemistry_element_O] / Amgsil, dp->si_sub * dp->atomic_weight[chemistry_element_Si] / Amgsil, 2 * dp->atomic_weight[chemistry_element_Mg] / Amgsil}; double gcomp3[3] = { dp->o_sub * dp->atomic_weight[chemistry_element_O] / Afesil, dp->si_sub * dp->atomic_weight[chemistry_element_Si] / Afesil, 2 * dp->atomic_weight[chemistry_element_Fe] / Afesil}; /* allocate memory */ dp->grain_element_mfrac[dust_grain_species_graphite_l] = malloc( dp->grain_element_count[dust_grain_species_graphite_l] * sizeof(double)); dp->grain_element_mfrac[dust_grain_species_mgsilicate_l] = malloc(dp->grain_element_count[dust_grain_species_mgsilicate_l] * sizeof(double)); dp->grain_element_mfrac[dust_grain_species_fesilicate_l] = malloc(dp->grain_element_count[dust_grain_species_fesilicate_l] * sizeof(double)); dp->grain_element_mfrac[dust_grain_species_graphite_s] = malloc( dp->grain_element_count[dust_grain_species_graphite_s] * sizeof(double)); dp->grain_element_mfrac[dust_grain_species_mgsilicate_s] = malloc(dp->grain_element_count[dust_grain_species_mgsilicate_s] * sizeof(double)); dp->grain_element_mfrac[dust_grain_species_fesilicate_s] = malloc(dp->grain_element_count[dust_grain_species_fesilicate_s] * sizeof(double)); /* deep copy grain arrays */ memcpy( dp->grain_element_mfrac[dust_grain_species_graphite_l], &gcomp1, dp->grain_element_count[dust_grain_species_graphite_l] * sizeof(double)); memcpy(dp->grain_element_mfrac[dust_grain_species_mgsilicate_l], &gcomp2, dp->grain_element_count[dust_grain_species_mgsilicate_l] * sizeof(double)); memcpy(dp->grain_element_mfrac[dust_grain_species_fesilicate_l], &gcomp3, dp->grain_element_count[dust_grain_species_fesilicate_l] * sizeof(double)); memcpy( dp->grain_element_mfrac[dust_grain_species_graphite_s], &gcomp1, dp->grain_element_count[dust_grain_species_graphite_s] * sizeof(double)); memcpy(dp->grain_element_mfrac[dust_grain_species_mgsilicate_s], &gcomp2, dp->grain_element_count[dust_grain_species_mgsilicate_s] * sizeof(double)); memcpy(dp->grain_element_mfrac[dust_grain_species_fesilicate_s], &gcomp3, dp->grain_element_count[dust_grain_species_fesilicate_s] * sizeof(double)); /* set chemistry indices composition */ int gidx1[1] = {chemistry_element_C}; int gidx2[3] = {chemistry_element_O, chemistry_element_Si, chemistry_element_Mg}; int gidx3[3] = {chemistry_element_O, chemistry_element_Si, chemistry_element_Fe}; /* allocate memory */ dp->grain_element_indices[dust_grain_species_graphite_l] = malloc( dp->grain_element_count[dust_grain_species_graphite_l] * sizeof(int)); dp->grain_element_indices[dust_grain_species_mgsilicate_l] = malloc( dp->grain_element_count[dust_grain_species_mgsilicate_l] * sizeof(int)); dp->grain_element_indices[dust_grain_species_fesilicate_l] = malloc( dp->grain_element_count[dust_grain_species_fesilicate_l] * sizeof(int)); dp->grain_element_indices[dust_grain_species_graphite_s] = malloc( dp->grain_element_count[dust_grain_species_graphite_s] * sizeof(int)); dp->grain_element_indices[dust_grain_species_mgsilicate_s] = malloc( dp->grain_element_count[dust_grain_species_mgsilicate_s] * sizeof(int)); dp->grain_element_indices[dust_grain_species_fesilicate_s] = malloc( dp->grain_element_count[dust_grain_species_fesilicate_s] * sizeof(int)); /* deep copy grain arrays */ memcpy(dp->grain_element_indices[dust_grain_species_graphite_l], &gidx1, dp->grain_element_count[dust_grain_species_graphite_l] * sizeof(int)); memcpy( dp->grain_element_indices[dust_grain_species_mgsilicate_l], &gidx2, dp->grain_element_count[dust_grain_species_mgsilicate_l] * sizeof(int)); memcpy( dp->grain_element_indices[dust_grain_species_fesilicate_l], &gidx3, dp->grain_element_count[dust_grain_species_fesilicate_l] * sizeof(int)); memcpy(dp->grain_element_indices[dust_grain_species_graphite_s], &gidx1, dp->grain_element_count[dust_grain_species_graphite_s] * sizeof(int)); memcpy( dp->grain_element_indices[dust_grain_species_mgsilicate_s], &gidx2, dp->grain_element_count[dust_grain_species_mgsilicate_s] * sizeof(int)); memcpy( dp->grain_element_indices[dust_grain_species_fesilicate_s], &gidx3, dp->grain_element_count[dust_grain_species_fesilicate_s] * sizeof(int)); } /** * @brief Set jagged composite grain arrays in the dustevo_props structure * * @param dp the struct */ void dustevo_set_jagged_cgrain_arrays(struct dustevo_props *dp) { /* effective grain molecular weights in amu */ double Agra = dp->atomic_weight[chemistry_element_C]; double Asil = dp->fe_sub * dp->atomic_weight[chemistry_element_Fe] + dp->mg_sub * dp->atomic_weight[chemistry_element_Mg] + dp->si_sub * dp->atomic_weight[chemistry_element_Si] + dp->o_sub * dp->atomic_weight[chemistry_element_O]; /* set grain composition */ double gcomp1[1] = {dp->atomic_weight[chemistry_element_C] / Agra}; double gcomp2[4] = { dp->o_sub * dp->atomic_weight[chemistry_element_O] / Asil, dp->mg_sub * dp->atomic_weight[chemistry_element_Mg] / Asil, dp->si_sub * dp->atomic_weight[chemistry_element_Si] / Asil, dp->fe_sub * dp->atomic_weight[chemistry_element_Fe] / Asil}; /* allocate memory */ dp->cgrain_element_mfrac[dust_cgrain_species_graphite_l] = malloc(dp->cgrain_element_count[dust_cgrain_species_graphite_l] * sizeof(double)); dp->cgrain_element_mfrac[dust_cgrain_species_silicate_l] = malloc(dp->cgrain_element_count[dust_cgrain_species_silicate_l] * sizeof(double)); dp->cgrain_element_mfrac[dust_cgrain_species_graphite_s] = malloc(dp->cgrain_element_count[dust_cgrain_species_graphite_s] * sizeof(double)); dp->cgrain_element_mfrac[dust_cgrain_species_silicate_s] = malloc(dp->cgrain_element_count[dust_cgrain_species_silicate_s] * sizeof(double)); /* deep copy cgrain arrays */ memcpy(dp->cgrain_element_mfrac[dust_cgrain_species_graphite_l], &gcomp1, dp->cgrain_element_count[dust_cgrain_species_graphite_l] * sizeof(double)); memcpy(dp->cgrain_element_mfrac[dust_cgrain_species_silicate_l], &gcomp2, dp->cgrain_element_count[dust_cgrain_species_silicate_l] * sizeof(double)); memcpy(dp->cgrain_element_mfrac[dust_cgrain_species_graphite_s], &gcomp1, dp->cgrain_element_count[dust_cgrain_species_graphite_s] * sizeof(double)); memcpy(dp->cgrain_element_mfrac[dust_cgrain_species_silicate_s], &gcomp2, dp->cgrain_element_count[dust_cgrain_species_silicate_s] * sizeof(double)); /* 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->cgrain_element_indices[dust_cgrain_species_graphite_l] = malloc( dp->cgrain_element_count[dust_cgrain_species_graphite_l] * sizeof(int)); dp->cgrain_element_indices[dust_cgrain_species_silicate_l] = malloc( dp->cgrain_element_count[dust_cgrain_species_silicate_l] * sizeof(int)); dp->cgrain_element_indices[dust_cgrain_species_graphite_s] = malloc( dp->cgrain_element_count[dust_cgrain_species_graphite_s] * sizeof(int)); dp->cgrain_element_indices[dust_cgrain_species_silicate_s] = malloc( dp->cgrain_element_count[dust_cgrain_species_silicate_s] * sizeof(int)); /* deep copy cgrain arrays */ memcpy( dp->cgrain_element_indices[dust_cgrain_species_graphite_l], &gidx1, dp->cgrain_element_count[dust_cgrain_species_graphite_l] * sizeof(int)); memcpy( dp->cgrain_element_indices[dust_cgrain_species_silicate_l], &gidx2, dp->cgrain_element_count[dust_cgrain_species_silicate_l] * sizeof(int)); memcpy( dp->cgrain_element_indices[dust_cgrain_species_graphite_s], &gidx1, dp->cgrain_element_count[dust_cgrain_species_graphite_s] * sizeof(int)); memcpy( dp->cgrain_element_indices[dust_cgrain_species_silicate_s], &gidx2, dp->cgrain_element_count[dust_cgrain_species_silicate_s] * 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; /* set up species mass fraction arrays and check for dust before evolution * step */ double grain_mfracs[dust_cgrain_species_count] = {0.}; float D_pre = 0.; int gspec = 0; dust_check_for_NaNs(p) for (int grain = 0; grain < dust_grain_species_count; grain++) { /* sum for total mass fraction in dust before step */ D_pre += p->dust_data.grain_mass_fraction[grain]; grain_mfracs[gspec] = p->dust_data.grain_mass_fraction[grain]; if ((grain % dust_grain_molecule_count) == dust_grain_molecule_mgsilicate) { grain++; grain_mfracs[gspec] += p->dust_data.grain_mass_fraction[grain]; } gspec++; } /* no dust, no accretion or destruction. return. */ if (D_pre <= FLT_MIN) return; /* compute dust fractions in magnesium species */ double mg_specfrac_dust = 0.; const double totsil_mfrac = grain_mfracs[dust_cgrain_species_silicate_l] + grain_mfracs[dust_cgrain_species_silicate_s]; if (totsil_mfrac > 0) { mg_specfrac_dust = (p->dust_data.grain_mass_fraction[dust_grain_species_mgsilicate_l] + p->dust_data.grain_mass_fraction[dust_grain_species_mgsilicate_s]) / totsil_mfrac; } const double fe_specfrac_dust = (1. - mg_specfrac_dust); /* default subspecies composition in dust change */ double mg_specfrac_delta = mg_specfrac_dust; double fe_specfrac_delta = fe_specfrac_dust; /* compute current composition of silicate grains in particle */ double dust_composition[chemistry_element_count] = {0.}; dust_recompute_composition(dust_composition, mg_specfrac_dust, fe_specfrac_dust, dp); /* initialise array for storing diffuse metal fractions */ double diffuse_mfrac[chemistry_element_count]; for (int el = 0; el < chemistry_element_count; el++) { diffuse_mfrac[el] = p->chemistry_data.metal_mass_fraction[el]; } /* 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 --------------- */ /* base grain radius in cgs, 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); /* Decide what clumping factor to use... */ float clumping_factor; if (dp->clumping_factor_mode == clumping_mode_constant) { /* use a constant clumping factor */ clumping_factor = dp->clumping_factor; } else if (dp->clumping_factor_mode == clumping_mode_syncedchimes) { /* synchronise clumping factor with the cooling */ clumping_factor = cooling_get_dust_clumping_boost_factor(p, xp); } else if (dp->clumping_factor_mode == clumping_mode_variable) { /* compute an independent variable clumping factor */ if (n_H_cgs <= dp->clumping_factor_nH_min_cgs) clumping_factor = 1.0f; else if (n_H_cgs >= dp->clumping_factor_nH_max_cgs) clumping_factor = dp->clumping_factor; else { const float log_nH = log10f(n_H_cgs); const float log_nH_min = log10f(dp->clumping_factor_nH_min_cgs); const float log_nH_max = log10f(dp->clumping_factor_nH_max_cgs); const float log_max_boost = log10f(dp->clumping_factor); clumping_factor = exp10f( ((log_nH - log_nH_min) / max(log_nH_max - log_nH_min, FLT_MIN)) * log_max_boost); } } else { error("Dust Error: invalid clumping_factor_mode"); clumping_factor = -1.0f; } /*------------------------- 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.0); const double dfdt_cgs = drdt_grain_cgs * dfdr_cgs; deltf_sput = dfdt_cgs * dt_cgs; } else { deltf_sput = 0.; } /*------------------------- SNII shocks -------------------------*/ double deltf_SNII; const float SFR = star_formation_get_SFR(p, xp); if (dp->with_SNII_destruction && (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 = SFR * dp->specific_numSNII / unit_time_cgs; /* Yamasawa+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_cgrain_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; if (dp->pair_to_cooling) { /* paired mode */ for (int cgrain_l = 0; cgrain_l < dust_grain_chemistry_count; cgrain_l++) { const int sgrain_l = cgrain_l; const int cgrain_s = cgrain_l + dust_grain_chemistry_count; const int sgrain_s = cgrain_l + dust_grain_molecule_count; float min_abundance_ratio = FLT_MAX; for (int elem = 0; elem < dp->cgrain_element_count[cgrain_l]; elem++) { const int eldx = dp->cgrain_element_indices[cgrain_l][elem]; /* find the minimum solar abundance ratio among constituent elements * indicating element that bottlenecks accretion of effective grain */ if (eldx == chemistry_element_Mg || eldx == chemistry_element_Fe) { min_abundance_ratio = min( min_abundance_ratio, ((dp->mg_to_fe_aweight * p->chemistry_data.metal_mass_fraction[chemistry_element_Mg]) + p->chemistry_data.metal_mass_fraction[chemistry_element_Fe]) / ((dp->mg_to_fe_aweight * dp->abundance_pattern[chemistry_element_Mg]) + dp->abundance_pattern[chemistry_element_Fe])); } else { min_abundance_ratio = min(min_abundance_ratio, p->chemistry_data.metal_mass_fraction[eldx] / dp->abundance_pattern[eldx]); } } /* base accretion timescale for grain */ deltf_acc[cgrain_l] = dt_cgs * clumping_factor * min_abundance_ratio / (dp->accretion_coeff[cgrain_l] * mult_terms * dp->grain_size_factors[sgrain_l]); deltf_acc[cgrain_s] = dt_cgs * clumping_factor * min_abundance_ratio / (dp->accretion_coeff[cgrain_s] * mult_terms * dp->grain_size_factors[sgrain_s]); } } else { /* unpaired mode - need to find diffuse metal mass fraction */ for (int cgrain_l = 0; cgrain_l < dust_grain_chemistry_count; cgrain_l++) { const int cgrain_s = cgrain_l + dust_grain_chemistry_count; const int sgrain_s = cgrain_l + dust_grain_molecule_count; const int sgrain_l = cgrain_l; float min_abundance_ratio = FLT_MAX; if (grain_mfracs[cgrain_l] > 0) { for (int elem = 0; elem < dp->cgrain_element_count[cgrain_l]; elem++) { const int eldx = dp->cgrain_element_indices[cgrain_l][elem]; /* find the minimum solar abundance ratio among constituent elements * indicating element that bottlenecks accretion of effective grain */ if (eldx == chemistry_element_Mg || eldx == chemistry_element_Fe) { diffuse_mfrac[chemistry_element_Mg] -= dust_composition[chemistry_element_Mg] * (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s]); diffuse_mfrac[chemistry_element_Fe] -= dust_composition[chemistry_element_Fe] * (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s]); min_abundance_ratio = min(min_abundance_ratio, ((dp->mg_to_fe_aweight * diffuse_mfrac[chemistry_element_Mg]) + diffuse_mfrac[chemistry_element_Fe]) / ((dp->mg_to_fe_aweight * dp->abundance_pattern[chemistry_element_Mg]) + dp->abundance_pattern[chemistry_element_Fe])); } else { diffuse_mfrac[eldx] -= dust_composition[eldx] * (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s]); min_abundance_ratio = min(min_abundance_ratio, diffuse_mfrac[eldx] / dp->abundance_pattern[eldx]); } } if (min_abundance_ratio > FLT_MIN) { /* only in this case modify deltf[grain] to non-zero value */ deltf_acc[cgrain_l] = dt_cgs * clumping_factor * min_abundance_ratio / (dp->accretion_coeff[cgrain_l] * mult_terms * dp->grain_size_factors[sgrain_l]); deltf_acc[cgrain_s] = dt_cgs * clumping_factor * min_abundance_ratio / (dp->accretion_coeff[cgrain_s] * mult_terms * dp->grain_size_factors[sgrain_s]); } } } } } dust_check_for_NaNs(p) /* ------------------------- Coagulation ------------------------- */ double deltf_coag[dust_grain_chemistry_count] = {0.}; if (dp->with_coagulation) { for (int gchem = 0; gchem < dust_grain_chemistry_count; gchem++) { if (grain_mfracs[gchem + dust_grain_chemistry_count] > FLT_MIN) { deltf_coag[gchem] = dt_cgs * clumping_factor / ((1e-2 / grain_mfracs[gchem + dust_grain_chemistry_count]) * (1.e3 / n_H_cgs) * dp->coagulation_timescale_norm_prefactor * coagulation_timescale_norm); } else { deltf_coag[gchem] = 0.; } } } dust_check_for_NaNs(p) /* ------------------------- Shattering -------------------------- */ double deltf_shat[dust_grain_chemistry_count] = {0.}; if (dp->with_shattering) { for (int gchem = 0; gchem < dust_grain_chemistry_count; gchem++) { if (grain_mfracs[gchem] > FLT_MIN) { if (n_H_cgs < 1.) { deltf_shat[gchem] = dt_cgs / ((1e-2 / grain_mfracs[gchem]) * (1. / n_H_cgs) * shattering_timescale_norm); } else { deltf_shat[gchem] = dt_cgs / ((1e-2 / grain_mfracs[gchem]) * (1. / cbrt(n_H_cgs)) * shattering_timescale_norm); } } else { deltf_shat[gchem] = 0.; } } } dust_check_for_NaNs(p); /* ------------------------- Combine ------------------------- */ double delta_dust_composition[chemistry_element_count] = {0.}; delta_dust_composition[chemistry_element_C] = 1.; for (int species = 0; species < dust_grain_chemistry_count; species++) { double dfrac_renorm = 1.; double delta_dfrac_max = 1.; const int cgrain_l = species; const int sgrain_l = species; const int cgrain_s = species + dust_grain_chemistry_count; const int sgrain_s = species + dust_grain_molecule_count; /* compute mass transfer between small and large grains */ const double delta_dfrac_shat = grain_mfracs[cgrain_l] * (expm1(deltf_shat[species])); const double delta_dfrac_coag = grain_mfracs[cgrain_s] * (expm1(deltf_coag[species])); double delta_dfrac_s_to_l = delta_dfrac_coag - delta_dfrac_shat; /* combine all mass change contributions */ double delta_dfrac_l = min( FLT_MAX, grain_mfracs[cgrain_l] * (expm1((deltf_sput / dp->grain_size_factors[sgrain_l]) + deltf_SNII + deltf_acc[cgrain_l]))); double delta_dfrac_s = grain_mfracs[cgrain_s] * min(FLT_MAX, expm1((deltf_sput / dp->grain_size_factors[sgrain_s]) + deltf_SNII + deltf_acc[cgrain_s])); double delta_dfrac_tot = delta_dfrac_l + delta_dfrac_s; /* limit grain growth */ if (delta_dfrac_tot > 0) { /* for net growing silicates, recompute composition for new grains */ if (species == dust_grain_chemistry_silicate) { const double specfrac_denominator = (diffuse_mfrac[chemistry_element_Mg] / dp->grain_element_mfrac[dust_grain_species_mgsilicate_l] [dust_mgsilicate_constituent_Mg]) + (diffuse_mfrac[chemistry_element_Fe] / dp->grain_element_mfrac[dust_grain_species_fesilicate_l] [dust_fesilicate_constituent_Fe]); if (specfrac_denominator > 0) { mg_specfrac_delta = (diffuse_mfrac[chemistry_element_Mg] / dp->grain_element_mfrac[dust_grain_species_mgsilicate_l] [dust_mgsilicate_constituent_Mg]) / specfrac_denominator; fe_specfrac_delta = 1. - mg_specfrac_delta; } dust_recompute_composition(delta_dust_composition, mg_specfrac_delta, fe_specfrac_delta, dp); } if (dp->pair_to_cooling) { for (int elem = 0; elem < dp->cgrain_element_count[cgrain_l]; elem++) { const int eldx = dp->cgrain_element_indices[cgrain_l][elem]; /* renormalise grain growth limiting available mass in a consitituent * element */ const double total_mfrac = diffuse_mfrac[eldx] + (dust_composition[eldx] * (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s])); const double diffuse_limit = dp->grain_element_min_diffuse[eldx] * total_mfrac; if (delta_dust_composition[eldx] > 0) { delta_dfrac_max = min(delta_dfrac_max, (diffuse_mfrac[eldx] - diffuse_limit) / delta_dust_composition[eldx]); } } } else { for (int elem = 0; elem < dp->cgrain_element_count[cgrain_l]; elem++) { if (grain_mfracs[cgrain_l] > 0) { const int eldx = dp->cgrain_element_indices[cgrain_l][elem]; const double diffuse_limit = dp->grain_element_min_diffuse[eldx] * p->chemistry_data.metal_mass_fraction[eldx]; if (delta_dust_composition[eldx] > 0) { delta_dfrac_max = min(delta_dfrac_max, (diffuse_mfrac[eldx] - diffuse_limit) / delta_dust_composition[eldx]); } } } } dfrac_renorm = min(dfrac_renorm, delta_dfrac_max / delta_dfrac_tot); if (dfrac_renorm != 0.) { delta_dfrac_l *= dfrac_renorm; delta_dfrac_s *= dfrac_renorm; } else { // delta_dfrac_s = delta_dfrac_max * (delta_dfrac_s / delta_dfrac_tot); // delta_dfrac_l = delta_dfrac_max * (delta_dfrac_l / delta_dfrac_tot); delta_dfrac_l = 0.; delta_dfrac_s = 0.; } } else { /* if not growing, composition of grain change is current * grain composition */ memcpy(delta_dust_composition, dust_composition, sizeof(dust_composition)); } const double new_grain_mfrac_l = grain_mfracs[cgrain_l] + delta_dfrac_l; const double new_grain_mfrac_s = grain_mfracs[cgrain_s] + delta_dfrac_s; const double new_grain_mfrac_tot = new_grain_mfrac_l + new_grain_mfrac_s; delta_dfrac_tot = delta_dfrac_l + delta_dfrac_s; /* ensure mass transfer beetween grain size is always constrained (wrt * rounding errors) */ if (delta_dfrac_s_to_l > new_grain_mfrac_s) { delta_dfrac_s_to_l = new_grain_mfrac_s; } else if (delta_dfrac_s_to_l < -new_grain_mfrac_l) { delta_dfrac_s_to_l = -new_grain_mfrac_l; } // MATTHIEU: Temporary check to track down possbile NaNs in dust. if (isnan(new_grain_mfrac_tot) || isnan(delta_dfrac_s_to_l) || isnan(grain_mfracs[cgrain_l]) || isnan(grain_mfracs[cgrain_s])) { error("nans appeared, crashing..."); } if (species == dust_grain_chemistry_silicate) { /* need the diffuse iron mass fraction from SNIa for both modes */ double diffuse_iron_mass_fraction_from_SNIa = p->chemistry_data.iron_mass_fraction_from_SNIa; if (dp->pair_to_cooling == 0) { diffuse_iron_mass_fraction_from_SNIa -= p->dust_data.grain_iron_mass_fraction_from_SNIa; } /* need ratio of diffuse mass in iron to mass in iron from SNIa */ const double delta_depFefrac_tot = delta_dfrac_tot * delta_dust_composition[chemistry_element_Fe]; double ratio_iron_SNIa = 0; /* use initial zero value if ratio denominator go below min float limit */ if ((delta_depFefrac_tot > 0) && diffuse_mfrac[chemistry_element_Fe] > FLT_MIN) { /* transferring diffuse -> depleted iron, use diffuse SNIa Fe ratio */ ratio_iron_SNIa = diffuse_iron_mass_fraction_from_SNIa / diffuse_mfrac[chemistry_element_Fe]; } else if ((delta_depFefrac_tot < 0) && (dust_composition[chemistry_element_Fe] * (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s]) > FLT_MIN)) { /* transferring depleted -> diffuse iron, use depleted SNIa Fe ratio */ ratio_iron_SNIa = p->dust_data.grain_iron_mass_fraction_from_SNIa / (dust_composition[chemistry_element_Fe] * (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s])); } /* knowing the fraction of iron in SNIa, mass fraction change in * SNIa iron */ const double delta_depleted_iron_SNIa = delta_dfrac_tot * ratio_iron_SNIa * delta_dust_composition[chemistry_element_Fe]; if (dp->pair_to_cooling) { /* in paired mode, tracer represents diffuse */ p->chemistry_data.iron_mass_fraction_from_SNIa -= delta_depleted_iron_SNIa; p->dust_data.grain_iron_mass_fraction_from_SNIa += delta_depleted_iron_SNIa; } } const double species_mfrac_delta = (max(0., new_grain_mfrac_l + delta_dfrac_s_to_l) + max(0., new_grain_mfrac_s - delta_dfrac_s_to_l)) - (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s]); if (dp->pair_to_cooling) { for (int elem = 0; elem < dp->cgrain_element_count[cgrain_l]; elem++) { const int eldx = dp->cgrain_element_indices[cgrain_l][elem]; /* update constituent element abundances in paired mode */ p->chemistry_data.metal_mass_fraction[eldx] -= species_mfrac_delta * delta_dust_composition[eldx]; } } if (species != dust_grain_chemistry_silicate) { /* update particle grain mass fractions */ p->dust_data.grain_mass_fraction[sgrain_l] = max(0., new_grain_mfrac_l + delta_dfrac_s_to_l); p->dust_data.grain_mass_fraction[sgrain_s] = max(0., new_grain_mfrac_s - delta_dfrac_s_to_l); } else { double mg_specfrac_new = mg_specfrac_dust; double fe_specfrac_new = fe_specfrac_dust; if (delta_dfrac_tot > 0) { /* in net accreting case need to calculate new silicate composition */ mg_specfrac_new = (delta_dfrac_tot * mg_specfrac_delta + (grain_mfracs[cgrain_l] + grain_mfracs[cgrain_s]) * mg_specfrac_dust) / new_grain_mfrac_tot; fe_specfrac_new = 1. - mg_specfrac_new; } p->dust_data.grain_mass_fraction[sgrain_l] = max(0., new_grain_mfrac_l + delta_dfrac_s_to_l) * mg_specfrac_new; p->dust_data.grain_mass_fraction[sgrain_s] = max(0., new_grain_mfrac_s - delta_dfrac_s_to_l) * mg_specfrac_new; p->dust_data.grain_mass_fraction[sgrain_l + 1] = max(0., new_grain_mfrac_l + delta_dfrac_s_to_l) * fe_specfrac_new; p->dust_data.grain_mass_fraction[sgrain_s + 1] = max(0., new_grain_mfrac_s - delta_dfrac_s_to_l) * fe_specfrac_new; } } dust_check_for_NaNs(p); }