diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml index ea84eeac27603d0c5d8185c44ecf1c0374ca7c76..c2ada783a316a08ee729907ba69e9a66aebe6910 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 bdd559c9b8aab718962398f6f6e06c43f94f9df6..42abab62409b2ce61a5bd12bc745d63a62cb79b9 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 efaa5bf2d172f49954a4e8e118b8b8d555a6b380..19c5c95be1f1a9fcaa4e790fb24fd0f60d1ef4c2 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 ae9fcc31f5f8bd91d00744d7c96c7b4095356026..35a30d8301d92120498da3aa358c1cc3d617eee4 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 57c1fc877dcfc512d735d64bacfa88e8fd693d23..b7a17bdb313f3b5ddc2416cdc7a729e4f8916ffe 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 d268486631acd21522b7036259d4d3431e0440d1..df0a4599a0f81eaf619b078b86c0ade4e1eba7b8 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 0e7afbd02c0dfd0202158b867281ae2c66840e7f..e76e082d9b92bcde800b084feea44515b9579dc8 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 de8badb08eb3742cc7105b84f673a517c50b81ee..0718cfd59fd2e86c4eb968241d5418dfd52045d9 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 d0aa9c43b95c0e48d7d9fdba768af1ded557873e..401a68cdbf09cc59028c66d10833e19a95b5a297 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 d4aa7fb10a5cb1d872dde201e8ec6caa3c9c62f0..ec183f7701dd6da0bcea12eff23cad2ee30ab97e 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 af56214858a4773583cc6ef8c224fbf8097cc9cc..8ac5f6924c6893717033cd5f144fe873f6f17e24 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 aae8e1aa46453ef8c5d952f40db500a124100e98..6881e7d0093583b8ba60b3d9723a2214ef677f8a 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 04f7e2ef6b12fc6e5f078127b6f49c18906ac31c..f4fe4ac4df16af53c88a901edc97e5d70f4ae806 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 2ec103ae36969532e8246328ab32cb5428ea842a..29e00435e16751f01bdb58fe26af151b422bd5a1 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 f24986c7f1e2c326bef663ded0578c347c32da84..861159ec1022fcf6f7617fc6ca0c9b470aa6e7fa 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 cef4e3e6f3aa00c40ba92e3980b9fc472b6898bc..b8a6e12d1085a8fbc950bc36ceff93d2d5873754 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 a571f2b24d57b2720c3c77ebd7600a3830e4d2a3..02103022c6d7e38ae5e8d48054b45d60edf1a45b 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 68d24ae689d1b0c36541b6ec96545276b2191363..c772f86a083230fcc3a6cc1c9659dfd815b2fd2c 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 2d2a8e8ff43b5bfb586beada0a40d0af4d4e15fb..3edf81b50ea14bc9400fa9091ec369fffee15e82 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);