### Merge branch 'master' into timestep_limiter_upgrade

parents 593763c1 5b64c24a
 .. Basic sub-grid model Matthieu Schaller, 20th December 2018 Basic model =========== Cooling: Analytic models ~~~~~~~~~~~~~~~~~~~~~~~~ Currently, we have 3 different simple cooling models (const-lambda, const-du and Compton). These are all based on analytic formulas and can be used to quickly understand how the cooling interacts with the rest of the code before moving to more complex models. Equations --------- The first table compares the different analytical cooling while the next ones are specific to a given cooling. The quantities are the internal energy (\$$u \$$), the density \$$rho \$$, the element mass fraction (\$$X_i \$$), the cooling function (\$$\\Lambda\$$, the proton mass (\$$m_H \$$) and the time step condition (\$$t\_\\text{step}\$$). If not specified otherwise, all cooling contains a temperature floor avoiding negative temperature. .. csv-table:: Analytical Cooling :header: "Variable", "Const-Lambda", "Const-du" "\$$\\frac{ \\mathrm{d}u }{ \\mathrm{d}t } \$$", "\$$-\\Lambda \\frac{\\rho^2 X_H^2}{\\rho m_H^2} \$$", "const" "\$$\\Delta t\_\\text{max} \$$", "\$$t\_\\text{step} \\frac{u}{\\left|\\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \$$", "\$$t\_\\text{step} \\frac{u}{\\ \\left| \\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \$$" TODO: Add description of the parameters and units. TODO: Add Compton cooling model How to Implement a New Cooling ------------------------------ The developer should provide at least one function for: * writing the cooling name in HDF5 * cooling a particle * the maximal time step possible * initializing a particle * computing the total energy radiated by a particle * initializing the cooling parameters * printing the cooling type For implementation details, see src/cooling/none/cooling.h See :ref:new_option for the full list of changes required.
 .. EAGLE sub-grid model Matthieu Schaller, 20th December 2018 EAGLE model =========== Chemical tracers ~~~~~~~~~~~~~~~~ The gas particles in the EAGLE model carry metal abundance information in the form of metal mass fractions. We follow the following 9 elements: H, He, C, N, O, Ne, Mg, Si and Fe. We additionally follow the total metal mass fraction Z. This is typically larger than the sum of the 7 metals that are individually traced since this will also contain the contribution of all the elements that are not individually followed. As part of the diagnostics, we additionally trace the elements coming from the different stellar evolution channels. We store for each particle the total mass coming from all the SNIa that enriched that particle and the metal mass fraction from SNIa. This is the fraction of the *total* gas mass that is in the form of metals originating from SNIa stars. By construction this fraction will be smaller than the total metal mass fraction. The same tracers exist for the SNII and AGB channels. Finally, we also compute the iron gas fraction from SNIa. This it the fraction of the *total* gas mass that is made of iron originating from SNIa explosions. We finally also compute the smoothed versions of the individual element mass fractions, of the total metal mass fractions, and of the iron gas fraction from SNIa. The chemistry module in src/chemistry/EAGLE includes all the arrays that are added to the particles and the functions used to compute the smoothed elements. In the snapshots, we output: +------------------------------+--------------------------------+-----------+-----------------------------+ | Name | Description | Units | Comments | +==============================+================================+===========+=============================+ | ElementAbundance | | Fraction of the gas mass | [-] | | Array of length | | | | in the different elements | | | 9 for each particle | +------------------------------+--------------------------------+-----------+-----------------------------+ | SmoothedElementAbundance | | Fraction of the gas mass | [-] | | Array of length | | | | in the different elements | | | 9 for each particle | | | | smoothed over SPH neighbours | | | +------------------------------+--------------------------------+-----------+-----------------------------+ | Metallicity | | Fraction of the gas mass | [-] | | | | | in *all* metals | | | +------------------------------+--------------------------------+-----------+-----------------------------+ | SmoothedMetallicity | | Fraction of the gas mass | [-] | | | | | in *all* metals | | | | | | smoothed over SPH neighbours | | | +------------------------------+--------------------------------+-----------+-----------------------------+ Cooling: Wiersma+2008a ~~~~~~~~~~~~~~~~~~~~~~ Particle tracers ~~~~~~~~~~~~~~~~ Star formation: Schaye+2008 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Stellar enrichment: Wiersma+2008b ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Supernova feedback: Schaye+2012 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Black-hole creation ~~~~~~~~~~~~~~~~~~~ Black-hole accretion ~~~~~~~~~~~~~~~~~~~~ AGN feedback ~~~~~~~~~~~~
 .. Equation of State Loic Hausammann, 7th April 2018 .. GEAR sub-grid model Matthieu Schaller, 20th December 2018 .. _cooling: Cooling ======= GEAR model =========== Currently, we have 5 different cooling (EAGLE, Grackle, const-lambda, const-du and none). Three of them are easily solved analytically (const-lambda, const-du and none) while the two last requires complex chemical networks. Equations --------- The first table compares the different analytical cooling while the next ones are specific to a given cooling. The quantities are the internal energy (\$$u \$$), the density \$$rho \$$, the element mass fraction (\$$X_i \$$), the cooling function (\$$\\Lambda\$$, the proton mass (\$$m_H \$$) and the time step condition (\$$t\_\\text{step}\$$). If not specified otherwise, all cooling contains a temperature floor avoiding negative temperature. .. csv-table:: Analytical Cooling :header: "Variable", "Const-Lambda", "Const-du", "None" "\$$\\frac{ \\mathrm{d}u }{ \\mathrm{d}t } \$$", "\$$-\\Lambda \\frac{\\rho^2 X_H^2}{\\rho m_H^2} \$$", "const", "0" "\$$\\Delta t\_\\text{max} \$$", "\$$t\_\\text{step} \\frac{u}{\\left|\\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \$$", "\$$t\_\\text{step} \\frac{u}{\\ \\left| \\frac{ \\mathrm{d}u }{ \\mathrm{d}t }\\right|} \$$", "None" Grackle ~~~~~~~ Cooling: Grackle ~~~~~~~~~~~~~~~~ Grackle is a chemistry and cooling library presented in B. Smith et al. 2016 _ (do not forget to cite if used). Four different modes are available: ... ... @@ -57,24 +35,3 @@ Finish with make machine-linux-gnu; make && make install. If you encounter any problem, you can look at the Grackle documentation _ You can now provide the path given for MACH_INSTALL_PREFIX to with-grackle. Eagle ~~~~~ TODO How to Implement a New Cooling ------------------------------ The developer should provide at least one function for: * writing the cooling name in HDF5 * cooling a particle * the maximal time step possible * initializing a particle * computing the total energy radiated by a particle * initializing the cooling parameters * printing the cooling type For implementation details, see src/cooling/none/cooling.h See :ref:new_option for the full list of changes required.
 .. Subgrid Models Matthieu Schaller, 20th December 2018 .. _subgrid: Galaxy Formation Subgrid Models =============================== Multiple models are available in SWIFT. The 'Basic' model can be use as an empty canvas to be copied to create additional models. .. toctree:: :maxdepth: 2 :caption: Available models: Basic/index EAGLE/index GEAR/index
 ... ... @@ -19,7 +19,7 @@ difference is the parameter file that will need to be adapted for SWIFT. ParameterFiles/index InitialConditions/index HydroSchemes/index Cooling/index SubgridModels/index EquationOfState/index ExternalPotentials/index NewOption/index ... ...
 ... ... @@ -70,8 +70,6 @@ EAGLEChemistry: InitAbundance_Magnesium: 5.907064e-4 InitAbundance_Silicon: 6.825874e-4 InitAbundance_Iron: 1.1032152e-3 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 GearChemistry: InitialMetallicity: 0.01295
 ... ... @@ -26,9 +26,6 @@ EAGLEChemistry: InitAbundance_Magnesium: 5.907064e-4 InitAbundance_Silicon: 6.825874e-4 InitAbundance_Iron: 1.1032152e-3 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 EagleCooling: filename: ./coolingtables/ ... ...
 ... ... @@ -68,8 +68,6 @@ EAGLEChemistry: # Solar abundances InitAbundance_Magnesium: 5.907064e-4 InitAbundance_Silicon: 6.825874e-4 InitAbundance_Iron: 1.1032152e-3 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 EagleCooling: filename: ./coolingtables/ ... ...
 ... ... @@ -77,8 +77,6 @@ EAGLEChemistry: # Solar abundances InitAbundance_Magnesium: 5.907064e-4 InitAbundance_Silicon: 6.825874e-4 InitAbundance_Iron: 1.1032152e-3 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 EagleCooling: filename: ./coolingtables/ ... ...
 ... ... @@ -70,8 +70,6 @@ EAGLEChemistry: # Solar abundances InitAbundance_Magnesium: 5.907064e-4 InitAbundance_Silicon: 6.825874e-4 InitAbundance_Iron: 1.1032152e-3 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 EagleCooling: filename: ./coolingtables/ ... ...
 ... ... @@ -81,8 +81,6 @@ EAGLEChemistry: # Solar abundances InitAbundance_Magnesium: 5.907064e-4 InitAbundance_Silicon: 6.825874e-4 InitAbundance_Iron: 1.1032152e-3 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 EagleCooling: filename: ./coolingtables/ ... ...
 ... ... @@ -70,8 +70,6 @@ EAGLEChemistry: InitAbundance_Magnesium: 0.0 InitAbundance_Silicon: 0.0 InitAbundance_Iron: 0.0 CalciumOverSilicon: 0.0 SulphurOverSilicon: 0.0 EagleCooling: filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/ ... ...
 ... ... @@ -82,5 +82,3 @@ EAGLEChemistry: InitAbundance_Magnesium: 0. InitAbundance_Silicon: 0. InitAbundance_Iron: 0. CalciumOverSilicon: 0. SulphurOverSilicon: 0.
 ... ... @@ -467,6 +467,13 @@ int main(int argc, char *argv[]) { error("VEOCIraptor not yet enabled over MPI."); #endif /* Temporary early aborts for modes not supported with hand-vec. */ #if defined(WITH_VECTORIZATION) && !defined(CHEMISTRY_NONE) error( "Cannot run with chemistry and hand-vectorization (yet). " "Use --disable-hand-vec at configure time."); #endif /* Check that we can write the snapshots by testing if the output * directory exists and is searchable and writable. */ char basename[PARSER_MAX_LINE_SIZE]; ... ...
 ... ... @@ -248,12 +248,15 @@ LambdaCooling: lambda_nH2_cgs: 1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3]) cooling_tstep_mult: 1.0 # (Optional) Dimensionless pre-factor for the time-step condition. # Parameters of the EAGLE cooling model (Wiersma+08 cooling tables). EagleCooling: filename: ./coolingtables/ # Location of the Wiersma+08 cooling tables reionisation_redshift: 11.5 # Redshift of Hydrogen re-ionization He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian He_reion_ev_pH: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom filename: ./coolingtables/ # Location of the Wiersma+08 cooling tables reionisation_redshift: 8.5 # Redshift of Hydrogen re-ionization He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian He_reion_ev_pH: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom CalciumOverSiliconInSolar: 1. # (Optional) Ratio of Ca/Si to use in units of solar. If set to 1, the code uses [Ca/Si] = 0, i.e. Ca/Si = 0.0941736. SulphurOverSiliconInSolar: 1. # (Optional) Ratio of S/Si to use in units of solar. If set to 1, the code uses [S/Si] = 0, i.e. S/Si = 0.6054160. # Cooling with Grackle 3.0 GrackleCooling: ... ... @@ -282,8 +285,6 @@ EAGLEChemistry: InitAbundance_Magnesium: 0.000 # Inital fraction of particle mass in Magnesium InitAbundance_Silicon: 0.000 # Inital fraction of particle mass in Silicon InitAbundance_Iron: 0.000 # Inital fraction of particle mass in Iron CalciumOverSilicon: 0.0941736 # Constant ratio of Calcium over Silicon abundance SulphurOverSilicon: 0.6054160 # Constant ratio of Sulphur over Silicon abundance # Structure finding options (requires velociraptor) StructureFinding: ... ...
 ... ... @@ -3438,6 +3438,8 @@ int cell_has_tasks(struct cell *c) { */ void cell_drift_part(struct cell *c, const struct engine *e, int force) { const int periodic = e->s->periodic; const double dim = {e->s->dim, e->s->dim, e->s->dim}; const int with_cosmology = (e->policy & engine_policy_cosmology); const float hydro_h_max = e->hydro_properties->h_max; const integertime_t ti_old_part = c->hydro.ti_old_part; ... ... @@ -3550,35 +3552,27 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { } #endif #ifdef PLANETARY_SPH /* Remove particles that cross the non-periodic box edge */ if (!(e->s->periodic)) { for (int i = 0; i < 3; i++) { if ((p->x[i] - xp->v_full[i] * dt_drift > e->s->dim[i]) || (p->x[i] - xp->v_full[i] * dt_drift < 0.f) || ((p->mass != 0.f) && ((p->x[i] < 0.01f * e->s->dim[i]) || (p->x[i] > 0.99f * e->s->dim[i])))) { /* (TEMPORARY) Crudely stop the particle manually */ message( "Particle %lld hit a box edge. \n" " pos=%.4e %.4e %.4e vel=%.2e %.2e %.2e", p->id, p->x, p->x, p->x, p->v, p->v, p->v); for (int j = 0; j < 3; j++) { p->v[j] = 0.f; p->gpart->v_full[j] = 0.f; xp->v_full[j] = 0.f; } p->h = hydro_h_max; p->time_bin = time_bin_inhibited; p->gpart->time_bin = time_bin_inhibited; hydro_part_has_no_neighbours(p, xp, e->cosmology); p->mass = 0.f; p->gpart->mass = 0.f; break; } /* In non-periodic BC runs, remove particles that crossed the border */ if (!periodic) { /* Did the particle leave the box? */ if ((p->x > dim) || (p->x < 0.) || // x (p->x > dim) || (p->x < 0.) || // y (p->x > dim) || (p->x < 0.)) { // z /* One last action before death? */ hydro_remove_part(p, xp); /* Remove the particle entirely */ struct gpart *gp = p->gpart; cell_remove_part(e, c, p, xp); /* and it's gravity friend */ if (gp != NULL) cell_remove_gpart(e, c, gp); continue; } } #endif /* Limit h to within the allowed range */ p->h = min(p->h, hydro_h_max); ... ... @@ -3633,6 +3627,9 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { */ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) { const int periodic = e->s->periodic; const double dim = {e->s->dim, e->s->dim, e->s->dim}; const int with_cosmology = (e->policy & engine_policy_cosmology); const float stars_h_max = e->stars_properties->h_max; const integertime_t ti_old_gpart = c->grav.ti_old_part; const integertime_t ti_current = e->ti_current; ... ... @@ -3699,11 +3696,12 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) { /* Drift from the last time the cell was drifted to the current time */ double dt_drift; if (e->policy & engine_policy_cosmology) if (with_cosmology) { dt_drift = cosmology_get_drift_factor(e->cosmology, ti_old_gpart, ti_current); else } else { dt_drift = (ti_current - ti_old_gpart) * e->time_base; } /* Loop over all the g-particles in the cell */ const size_t nr_gparts = c->grav.count; ... ... @@ -3727,25 +3725,20 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) { } #endif #ifdef PLANETARY_SPH /* Remove particles that cross the non-periodic box edge */ if (!(e->s->periodic)) { for (int i = 0; i < 3; i++) { if ((gp->x[i] - gp->v_full[i] * dt_drift > e->s->dim[i]) || (gp->x[i] - gp->v_full[i] * dt_drift < 0.f) || ((gp->mass != 0.f) && ((gp->x[i] < 0.01f * e->s->dim[i]) || (gp->x[i] > 0.99f * e->s->dim[i])))) { /* (TEMPORARY) Crudely stop the particle manually */ for (int j = 0; j < 3; j++) { gp->v_full[j] = 0.f; } gp->time_bin = time_bin_inhibited; gp->mass = 0.f; break; } /* In non-periodic BC runs, remove particles that crossed the border */ if (!periodic) { /* Did the particle leave the box? */ if ((gp->x > dim) || (gp->x < 0.) || // x (gp->x > dim) || (gp->x < 0.) || // y (gp->x > dim) || (gp->x < 0.)) { // z /* Remove the particle entirely */ if (gp->type == swift_type_dark_matter) cell_remove_gpart(e, c, gp); continue; } } #endif /* Init gravity force fields. */ if (gpart_is_active(gp, e)) { ... ... @@ -3775,6 +3768,25 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) { } #endif /* In non-periodic BC runs, remove particles that crossed the border */ if (!periodic) { /* Did the particle leave the box? */ if ((sp->x > dim) || (sp->x < 0.) || // x (sp->x > dim) || (sp->x < 0.) || // y (sp->x > dim) || (sp->x < 0.)) { // z /* Remove the particle entirely */ struct gpart *gp = sp->gpart; cell_remove_spart(e, c, sp); /* and it's gravity friend */ cell_remove_gpart(e, c, gp); continue; } } /* Limit h to within the allowed range */ sp->h = min(sp->h, stars_h_max); ... ...
 ... ... @@ -188,12 +188,6 @@ static INLINE void chemistry_init_backend(struct swift_params* parameter_file, data->initial_metal_mass_fraction[elem] = parser_get_param_float(parameter_file, buffer); } /* Read the constant ratios */ data->calcium_over_silicon_ratio = parser_get_param_float( parameter_file, "EAGLEChemistry:CalciumOverSilicon"); data->sulphur_over_silicon_ratio = parser_get_param_float( parameter_file, "EAGLEChemistry:SulphurOverSilicon"); } } ... ...
 ... ... @@ -45,12 +45,6 @@ struct chemistry_global_data { /*! Fraction of the particle mass in *all* metals at the start of the run */ float initial_metal_mass_fraction_total; /*! Constant ratio of Calcium over Silicium */ float calcium_over_silicon_ratio; /*! Constant ratio of Calcium over Silicium */ float sulphur_over_silicon_ratio; }; /** ... ...
 ... ... @@ -466,7 +466,10 @@ void cooling_cool_part(const struct phys_const *restrict phys_const, (See cosmology theory document for the derivation) */ const double delta_redshift = -dt * cosmo->H * cosmo->a_inv; /* Get this particle's abundance ratios */ /* Get this particle's abundance ratios compared to solar * Note that we need to add S and Ca that are in the tables but not tracked * by the particles themselves. * The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */ float abundance_ratio[chemistry_element_count + 2]; abundance_ratio_to_solar(p, cooling, abundance_ratio); ... ... @@ -739,10 +742,10 @@ void cooling_init_backend(struct swift_params *parameter_file, cooling->cooling_table_path); cooling->H_reion_z = parser_get_param_float( parameter_file, "EagleCooling:reionisation_redshift"); cooling->calcium_over_silicon_ratio = parser_get_param_float( parameter_file, "EAGLEChemistry:CalciumOverSilicon"); cooling->sulphur_over_silicon_ratio = parser_get_param_float( parameter_file, "EAGLEChemistry:SulphurOverSilicon"); cooling->Ca_over_Si_ratio_in_solar = parser_get_opt_param_float( parameter_file, "EAGLECooling:CalciumOverSiliconInSolar", 1.f); cooling->S_over_Si_ratio_in_solar = parser_get_opt_param_float( parameter_file, "EAGLECooling:SulphurOverSiliconInSolar", 1.f); cooling->He_reion_z_centre = parser_get_param_float(parameter_file, "EagleCooling:He_reion_z_centre"); cooling->He_reion_z_sigma = ... ... @@ -866,6 +869,7 @@ void cooling_clean(struct cooling_function_data *cooling) { free(cooling->HeFrac); free(cooling->Therm); free(cooling->SolarAbundances); free(cooling->SolarAbundances_inv); /* Free the tables */ free(cooling->table.metal_heating); ... ...
 ... ... @@ -29,50 +29,68 @@ #include "interpolate.h" /** * @brief Calculate ratio of particle element abundances * to solar abundance. * Multiple if statements are necessary because order of elements * in tables is different from chemistry_element enum. * Tables: H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe * Enum: H, He, C, N, O, Ne, Mg, Si, Fe * The order in ratio_solar is: * H, He, C, N, O, Ne, Mg, Si, Fe, S, Ca * Hence Fe, S, Ca need to be treated separately to be put in the * correct place in the output array. * * @param p Pointer to #part struct * @param cooling #cooling_function_data struct * @param ratio_solar Pointer to array or ratios * @brief Compute ratio of mass fraction to solar mass fraction * for each element carried by a given particle. * * The solar abundances are taken from the tables themselves. * * The EAGLE chemistry model does not track S and Ca. We assume * that their abundance with respect to solar is the same as * the ratio for Si. * We optionally apply a correction if the user asked for a different * ratio. * * We also re-order the elements such that they match the order of the * tables. This is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe]. * * The solar abundances table is arranged as * [H, He, C, N, O, Ne, Mg, Si, S, C, Fe]. * * @param p Pointer to #part struct. * @param cooling #cooling_function_data struct. * @param ratio_solar (return) Array of ratios to solar abundances. */ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar( const struct part *p, const struct cooling_function_data *cooling, float ratio_solar[chemistry_element_count + 2]) { /* compute ratios for all elements */ for (enum chemistry_element elem = chemistry_element_H; elem < chemistry_element_count; elem++) { if (elem == chemistry_element_Fe) { /* NOTE: solar abundances have iron last with calcium and sulphur directly * before, hence +2 */ ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /