From adf069f11205ce423f102720f7f6553c28f89648 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <matthieu.schaller@durham.ac.uk> Date: Thu, 18 Jan 2018 17:43:56 +0100 Subject: [PATCH] Add the possibility to read basic chemistry information from the parameter file and propagate this onto the particle at initialisation time. More uniform naming convention for the first initialisation of properties. --- examples/CoolingBox/coolingBox.yml | 14 +++++++ examples/main.c | 6 ++- examples/parameter_example.yml | 18 +++++++++ src/chemistry.c | 12 ++++-- src/chemistry.h | 5 ++- src/chemistry/EAGLE/chemistry.h | 52 ++++++++++++++++++++------ src/chemistry/EAGLE/chemistry_struct.h | 20 +++++++++- src/chemistry/gear/chemistry.h | 13 +++++-- src/chemistry/gear/chemistry_struct.h | 5 +++ src/chemistry/none/chemistry.h | 13 +++++-- src/chemistry/none/chemistry_struct.h | 7 ++++ src/cooling/const_du/cooling.h | 6 ++- src/cooling/const_lambda/cooling.h | 6 ++- src/cooling/grackle/cooling.h | 6 ++- src/cooling/none/cooling.h | 6 ++- src/engine.c | 15 ++++++++ src/engine.h | 5 +++ src/space.c | 23 +++++------- src/space.h | 6 +++ 19 files changed, 189 insertions(+), 49 deletions(-) diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml index ea84eeac27..c2ada783a3 100644 --- a/examples/CoolingBox/coolingBox.yml +++ b/examples/CoolingBox/coolingBox.yml @@ -46,3 +46,17 @@ GrackleCooling: UVbackground: 0 # Enable or not the UV background GrackleRedshift: 0 # Redshift to use (-1 means time based redshift) GrackleHSShieldingDensityThreshold: 1.1708e-26 # self shielding threshold in internal system of units + +EAGLEChemistry: + InitMetallicity: 0. + InitAbundance_Hydrogen: 0.752 + InitAbundance_Helium: 0.248 + InitAbundance_Carbon: 0.000 + InitAbundance_Nitrogen: 0.000 + InitAbundance_Oxygen: 0.000 + InitAbundance_Neon: 0.000 + InitAbundance_Magnesium: 0.000 + InitAbundance_Silicon: 0.000 + InitAbundance_Iron: 0.000 + CalciumOverSilicon: 0.0941736 + SulphurOverSilicon: 0.6054160 diff --git a/examples/main.c b/examples/main.c index bdd559c9b8..42abab6240 100644 --- a/examples/main.c +++ b/examples/main.c @@ -613,7 +613,9 @@ int main(int argc, char *argv[]) { if (with_cooling && myrank == 0) cooling_print(&cooling_func); /* Initialise the chemistry */ - chemistry_print(); + struct chemistry_data chemistry; + chemistry_init(params, &us, &prog_const, &chemistry); + if (myrank == 0) chemistry_print(&chemistry); /* Initialise the feedback properties */ struct sourceterms sourceterms; @@ -639,7 +641,7 @@ int main(int argc, char *argv[]) { engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, N_total[0], N_total[1], with_aff, engine_policies, talking, &reparttype, &us, &prog_const, &hydro_properties, &gravity_properties, &potential, - &cooling_func, &sourceterms); + &cooling_func, &chemistry, &sourceterms); if (myrank == 0) { clocks_gettime(&toc); message("engine_init took %.3f %s.", clocks_diff(&tic, &toc), diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index efaa5bf2d1..19c5c95be1 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -151,3 +151,21 @@ GrackleCooling: UVbackground: 1 # Enable or not the UV background GrackleRedshift: 0 # Redshift to use (-1 means time based redshift) GrackleHSShieldingDensityThreshold: 1.1708e-26 # self shielding threshold in internal system of units + +# Parameters related to chemistry models ----------------------------------------------- + +# EAGLE model +EAGLEChemistry: + InitMetallicity: 0. # Inital fraction of particle mass in *all* metals + InitAbundance_Hydrogen: 0.752 # Inital fraction of particle mass in Hydrogen + InitAbundance_Helium: 0.248 # Inital fraction of particle mass in Helium + InitAbundance_Carbon: 0.000 # Inital fraction of particle mass in Carbon + InitAbundance_Nitrogen: 0.000 # Inital fraction of particle mass in Nitrogen + InitAbundance_Oxygen: 0.000 # Inital fraction of particle mass in Oxygen + InitAbundance_Neon: 0.000 # Inital fraction of particle mass in Neon + 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 + diff --git a/src/chemistry.c b/src/chemistry.c index ae9fcc31f5..35a30d8301 100644 --- a/src/chemistry.c +++ b/src/chemistry.c @@ -31,17 +31,23 @@ * @param parameter_file The parsed parameter file. * @param us The current internal system of units. * @param phys_const The physical constants in internal units. + * @param data The properties to initialise. */ void chemistry_init(const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const) { + const struct phys_const* phys_const, + struct chemistry_data* data) { - chemistry_init_backend(parameter_file, us, phys_const); + chemistry_init_backend(parameter_file, us, phys_const, data); } /** * @brief Prints the properties of the chemistry model to stdout. * * Calls chemistry_print_backend for the chosen chemistry model. + * + * @brief The #chemistry_data containing information about the current model. */ -void chemistry_print() { chemistry_print_backend(); } +void chemistry_print(const struct chemistry_data* data) { + chemistry_print_backend(data); +} diff --git a/src/chemistry.h b/src/chemistry.h index 57c1fc877d..b7a17bdb31 100644 --- a/src/chemistry.h +++ b/src/chemistry.h @@ -42,8 +42,9 @@ /* Common functions */ void chemistry_init(const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const); + const struct phys_const* phys_const, + struct chemistry_data* data); -void chemistry_print(); +void chemistry_print(const struct chemistry_data* data); #endif /* SWIFT_CHEMISTRY_H */ diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h index d268486631..df0a4599a0 100644 --- a/src/chemistry/EAGLE/chemistry.h +++ b/src/chemistry/EAGLE/chemistry.h @@ -44,8 +44,8 @@ __attribute__((always_inline)) INLINE static const char* chemistry_get_element_name(enum chemistry_element elem) { static const char* chemistry_element_names[chemistry_element_count] = { - "Hydrogen", "Helium", "Carbon", "Nitrogen", "Oxygen", - "Neon", "Magnesium", "Silicium", "Iron"}; + "Hydrogen", "Helium", "Carbon", "Nitrogen", "Oxygen", + "Neon", "Magnesium", "Silicon", "Iron"}; return chemistry_element_names[elem]; } @@ -54,33 +54,63 @@ chemistry_get_element_name(enum chemistry_element elem) { * @brief Sets the chemistry properties of the (x-)particles to a valid start * state. * - * Nothing to do here. - * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param data The global chemistry information. */ -__attribute__((always_inline)) INLINE static void chemistry_init_part( - const struct part* restrict p, struct xpart* restrict xp) {} +__attribute__((always_inline)) INLINE static void chemistry_first_init_part( + struct part* restrict p, struct xpart* restrict xp, + const struct chemistry_data* data) { + + p->chemistry_data.metal_mass_fraction_total = + data->initial_metal_mass_fraction_total; + for (int elem = 0; elem < chemistry_element_count; ++elem) + p->chemistry_data.metal_mass_fraction[elem] = + data->initial_metal_mass_fraction[elem]; +} /** * @brief Initialises the chemistry properties. * - * Nothing to do here. - * * @param parameter_file The parsed parameter file. * @param us The current internal system of units. * @param phys_const The physical constants in internal units. + * @param data The properties to initialise. */ static INLINE void chemistry_init_backend( const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const) {} + const struct phys_const* phys_const, struct chemistry_data* data) { + + /* Read the total metallicity */ + data->initial_metal_mass_fraction_total = + parser_get_param_float(parameter_file, "EAGLEChemistry:InitMetallicity"); + + /* Read the individual mass fractions */ + for (int elem = 0; elem < chemistry_element_count; ++elem) { + char buffer[50]; + sprintf(buffer, "EAGLEChemistry:InitAbundance_%s", + chemistry_get_element_name(elem)); + + 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"); +} /** * @brief Prints the properties of the chemistry model to stdout. + * + * @brief The #chemistry_data containing information about the current model. */ -static INLINE void chemistry_print_backend() { +static INLINE void chemistry_print_backend(const struct chemistry_data* data) { - message("Chemistry model is 'EAGLE'."); + message("Chemistry model is 'EAGLE' tracking %d elements.", + chemistry_element_count); } #endif /* SWIFT_CHEMISTRY_EAGLE_H */ diff --git a/src/chemistry/EAGLE/chemistry_struct.h b/src/chemistry/EAGLE/chemistry_struct.h index 0e7afbd02c..e76e082d9b 100644 --- a/src/chemistry/EAGLE/chemistry_struct.h +++ b/src/chemistry/EAGLE/chemistry_struct.h @@ -36,7 +36,25 @@ enum chemistry_element { }; /** - * @brief Chemical abundances traced in the EAGLE model. + * @brief Global chemical abundance information in the EAGLE model. + */ +struct chemistry_data { + + /*! Fraction of the particle mass in given elements at the start of the run */ + float initial_metal_mass_fraction[chemistry_element_count]; + + /*! 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; +}; + +/** + * @brief Chemical abundances traced by the #part in the EAGLE model. */ struct chemistry_part_data { diff --git a/src/chemistry/gear/chemistry.h b/src/chemistry/gear/chemistry.h index de8badb08e..0718cfd59f 100644 --- a/src/chemistry/gear/chemistry.h +++ b/src/chemistry/gear/chemistry.h @@ -45,9 +45,11 @@ * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param data The global chemistry information. */ -__attribute__((always_inline)) INLINE static void chemistry_init_part( - const struct part* restrict p, struct xpart* restrict xp) {} +__attribute__((always_inline)) INLINE static void chemistry_first_init_part( + const struct part* restrict p, struct xpart* restrict xp, + const struct chemistry_data* data) {} /** * @brief Initialises the chemistry properties. @@ -57,15 +59,18 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part( * @param parameter_file The parsed parameter file. * @param us The current internal system of units. * @param phys_const The physical constants in internal units. + * @param data The properties to initialise. */ static INLINE void chemistry_init_backend( const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const) {} + const struct phys_const* phys_const, struct chemistry_data* data) {} /** * @brief Prints the properties of the chemistry model to stdout. + * + * @brief The #chemistry_data containing information about the current model. */ -static INLINE void chemistry_print_backend() { +static INLINE void chemistry_print_backend(const struct chemistry_data* data) { message("Chemistry function is 'gear'."); } diff --git a/src/chemistry/gear/chemistry_struct.h b/src/chemistry/gear/chemistry_struct.h index d0aa9c43b9..401a68cdbf 100644 --- a/src/chemistry/gear/chemistry_struct.h +++ b/src/chemistry/gear/chemistry_struct.h @@ -24,6 +24,11 @@ */ enum chemistry_element { chemistry_element_count = 0 }; +/** + * @brief Global chemical abundance information. + */ +struct chemistry_data {}; + /** * @brief Properties of the chemistry function. */ diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h index d4aa7fb10a..ec183f7701 100644 --- a/src/chemistry/none/chemistry.h +++ b/src/chemistry/none/chemistry.h @@ -45,9 +45,11 @@ * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param data The global chemistry information used for this run. */ -__attribute__((always_inline)) INLINE static void chemistry_init_part( - const struct part* restrict p, struct xpart* restrict xp) {} +__attribute__((always_inline)) INLINE static void chemistry_first_init_part( + const struct part* restrict p, struct xpart* restrict xp, + const struct chemistry_data* data) {} /** * @brief Initialises the chemistry properties. @@ -57,15 +59,18 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part( * @param parameter_file The parsed parameter file. * @param us The current internal system of units. * @param phys_const The physical constants in internal units. + * @param data The global chemistry information (to be filled). */ static INLINE void chemistry_init_backend( const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const) {} + const struct phys_const* phys_const, struct chemistry_data* data) {} /** * @brief Prints the properties of the chemistry model to stdout. + * + * @brief The #chemistry_data containing information about the current model. */ -static INLINE void chemistry_print_backend() { +static INLINE void chemistry_print_backend(const struct chemistry_data* data) { message("Chemistry function is 'No chemistry'."); } diff --git a/src/chemistry/none/chemistry_struct.h b/src/chemistry/none/chemistry_struct.h index af56214858..8ac5f6924c 100644 --- a/src/chemistry/none/chemistry_struct.h +++ b/src/chemistry/none/chemistry_struct.h @@ -29,6 +29,13 @@ */ enum chemistry_element { chemistry_element_count = 0 }; +/** + * @brief Global chemical abundance information. + * + * Nothing here. + */ +struct chemistry_data {}; + /** * @brief Chemistry properties carried by the #part. * diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h index aae8e1aa46..6881e7d009 100644 --- a/src/cooling/const_du/cooling.h +++ b/src/cooling/const_du/cooling.h @@ -131,9 +131,11 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param cooling The properties of the cooling function. */ -__attribute__((always_inline)) INLINE static void cooling_init_part( - const struct part* restrict p, struct xpart* restrict xp) { +__attribute__((always_inline)) INLINE static void cooling_first_init_part( + const struct part* restrict p, struct xpart* restrict xp, + const struct cooling_function_data* cooling) { xp->cooling_data.radiated_energy = 0.f; } diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index 04f7e2ef6b..f4fe4ac4df 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -142,9 +142,11 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param cooling The properties of the cooling function. */ -__attribute__((always_inline)) INLINE static void cooling_init_part( - const struct part* restrict p, struct xpart* restrict xp) { +__attribute__((always_inline)) INLINE static void cooling_first_init_part( + const struct part* restrict p, struct xpart* restrict xp, + const struct cooling_function_data* cooling) { xp->cooling_data.radiated_energy = 0.f; } diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h index 2ec103ae36..29e00435e1 100644 --- a/src/cooling/grackle/cooling.h +++ b/src/cooling/grackle/cooling.h @@ -58,9 +58,11 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param cooling The properties of the cooling function. */ -__attribute__((always_inline)) INLINE static void cooling_init_part( - const struct part* restrict p, struct xpart* restrict xp) { +__attribute__((always_inline)) INLINE static void cooling_first_init_part( + const struct part* restrict p, struct xpart* restrict xp, + const struct cooling_function_data* cooling) { xp->cooling_data.radiated_energy = 0.f; } diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h index f24986c7f1..861159ec10 100644 --- a/src/cooling/none/cooling.h +++ b/src/cooling/none/cooling.h @@ -91,9 +91,11 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. + * @param cooling The properties of the cooling function. */ -__attribute__((always_inline)) INLINE static void cooling_init_part( - const struct part* restrict p, struct xpart* restrict xp) {} +__attribute__((always_inline)) INLINE static void cooling_first_init_part( + const struct part* restrict p, struct xpart* restrict xp, + const struct cooling_function_data* cooling) {} /** * @brief Returns the total radiated energy by this particle. diff --git a/src/engine.c b/src/engine.c index cef4e3e6f3..b8a6e12d10 100644 --- a/src/engine.c +++ b/src/engine.c @@ -4151,6 +4151,18 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs, struct clocks_time time1, time2; clocks_gettime(&time1); + /* Start by setting the particles in a good state */ + ticks tic = getticks(); + if (e->nodeID == 0) message("first init..."); + /* Set the particles in a state where they are ready for a run */ + space_first_init_parts(s, e->chemistry); + space_first_init_xparts(s, e->cooling_func); + space_first_init_gparts(s); + space_first_init_sparts(s); + if (e->verbose) + message("took %.3f %s.", clocks_from_ticks(getticks() - tic), + clocks_getunit()); + if (e->nodeID == 0) message("Computing initial gas densities."); /* Initialise the softening lengths */ @@ -5101,6 +5113,7 @@ void engine_unpin() { * @param gravity The #gravity_props used for this run. * @param potential The properties of the external potential. * @param cooling_func The properties of the cooling function. + * @param chemistry The chemistry information. * @param sourceterms The properties of the source terms function. */ void engine_init(struct engine *e, struct space *s, @@ -5113,6 +5126,7 @@ void engine_init(struct engine *e, struct space *s, const struct gravity_props *gravity, const struct external_potential *potential, const struct cooling_function_data *cooling_func, + const struct chemistry_data *chemistry, struct sourceterms *sourceterms) { /* Clean-up everything */ @@ -5173,6 +5187,7 @@ void engine_init(struct engine *e, struct space *s, e->gravity_properties = gravity; e->external_potential = potential; e->cooling_func = cooling_func; + e->chemistry = chemistry; e->sourceterms = sourceterms; e->parameter_file = params; #ifdef WITH_MPI diff --git a/src/engine.h b/src/engine.h index a571f2b24d..02103022c6 100644 --- a/src/engine.h +++ b/src/engine.h @@ -34,6 +34,7 @@ /* Includes. */ #include "barrier.h" +#include "chemistry_struct.h" #include "clocks.h" #include "collectgroup.h" #include "cooling_struct.h" @@ -277,6 +278,9 @@ struct engine { /* Properties of the cooling scheme */ const struct cooling_function_data *cooling_func; + /* Properties of the chemistry model */ + const struct chemistry_data *chemistry; + /* Properties of source terms */ struct sourceterms *sourceterms; @@ -307,6 +311,7 @@ void engine_init(struct engine *e, struct space *s, const struct gravity_props *gravity, const struct external_potential *potential, const struct cooling_function_data *cooling_func, + const struct chemistry_data *chemistry, struct sourceterms *sourceterms); void engine_launch(struct engine *e); void engine_prepare(struct engine *e); diff --git a/src/space.c b/src/space.c index 68d24ae689..c772f86a08 100644 --- a/src/space.c +++ b/src/space.c @@ -41,6 +41,7 @@ /* Local headers. */ #include "atomic.h" +#include "chemistry.h" #include "const.h" #include "cooling.h" #include "engine.h" @@ -2618,7 +2619,8 @@ void space_synchronize_particle_positions(struct space *s) { * * Calls hydro_first_init_part() on all the particles */ -void space_first_init_parts(struct space *s) { +void space_first_init_parts(struct space *s, + const struct chemistry_data *chemistry) { const size_t nr_parts = s->nr_parts; struct part *restrict p = s->parts; @@ -2638,6 +2640,9 @@ void space_first_init_parts(struct space *s) { hydro_first_init_part(&p[i], &xp[i]); + /* Also initialise the chemistry */ + chemistry_first_init_part(&p[i], &xp[i], chemistry); + #ifdef SWIFT_DEBUG_CHECKS p[i].ti_drift = 0; p[i].ti_kick = 0; @@ -2650,7 +2655,8 @@ void space_first_init_parts(struct space *s) { * * Calls cooling_init_xpart() on all the particles */ -void space_first_init_xparts(struct space *s) { +void space_first_init_xparts(struct space *s, + const struct cooling_function_data *cool_func) { const size_t nr_parts = s->nr_parts; struct part *restrict p = s->parts; @@ -2658,7 +2664,7 @@ void space_first_init_xparts(struct space *s) { for (size_t i = 0; i < nr_parts; ++i) { - cooling_init_part(&p[i], &xp[i]); + cooling_first_init_part(&p[i], &xp[i], cool_func); } } @@ -3002,17 +3008,6 @@ void space_init(struct space *s, const struct swift_params *params, hydro_space_init(&s->hs, s); - ticks tic = getticks(); - if (verbose) message("first init..."); - /* Set the particles in a state where they are ready for a run */ - space_first_init_parts(s); - space_first_init_xparts(s); - space_first_init_gparts(s); - space_first_init_sparts(s); - if (verbose) - message("took %.3f %s.", clocks_from_ticks(getticks() - tic), - clocks_getunit()); - /* Init the space lock. */ if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock."); diff --git a/src/space.h b/src/space.h index 2d2a8e8ff4..3edf81b50e 100644 --- a/src/space.h +++ b/src/space.h @@ -223,6 +223,12 @@ void space_synchronize_particle_positions(struct space *s); void space_do_parts_sort(); void space_do_gparts_sort(); void space_do_sparts_sort(); +void space_first_init_parts(struct space *s, + const struct chemistry_data *chemistry); +void space_first_init_xparts(struct space *s, + const struct cooling_function_data *cool_func); +void space_first_init_gparts(struct space *s); +void space_first_init_sparts(struct space *s); void space_init_parts(struct space *s, int verbose); void space_init_gparts(struct space *s, int verbose); void space_convert_quantities(struct space *s, int verbose); -- GitLab