diff --git a/configure.ac b/configure.ac index e730619ea17ce87df7b5c62225094e1d7d0e34db..58f25b1eedbab44f9e3f0ecbf157f6c394674f78 100644 --- a/configure.ac +++ b/configure.ac @@ -962,7 +962,7 @@ esac # Cooling function AC_ARG_WITH([cooling], [AS_HELP_STRING([--with-cooling=], - [cooling function @<:@none, const-du, const-lambda, grackle default: none@:>@] + [cooling function @<:@none, const-du, const-lambda, EAGLE, grackle, grackle1, grackle2, grackle3 default: none@:>@] )], [with_cooling="$withval"], [with_cooling="none"] @@ -979,7 +979,20 @@ case "$with_cooling" in ;; grackle) AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) + AC_DEFINE([COOLING_GRACKLE_MODE], [0], [Grackle chemistry network, mode 0]) + ;; + grackle1) + AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) + AC_DEFINE([COOLING_GRACKLE_MODE], [1], [Grackle chemistry network, mode 1]) ;; + grackle2) + AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) + AC_DEFINE([COOLING_GRACKLE_MODE], [2], [Grackle chemistry network, mode 2]) + ;; + grackle3) + AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library]) + AC_DEFINE([COOLING_GRACKLE_MODE], [3], [Grackle chemistry network, mode 3]) + ;; EAGLE) AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model]) ;; diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml index c2ada783a316a08ee729907ba69e9a66aebe6910..2bd2f19f6d78388ae638521f590255d410bc8697 100644 --- a/examples/CoolingBox/coolingBox.yml +++ b/examples/CoolingBox/coolingBox.yml @@ -42,11 +42,17 @@ LambdaCooling: # Cooling with Grackle 2.0 GrackleCooling: - GrackleCloudyTable: ../CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) - 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 - + CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) + WithUVbackground: 0 # Enable or not the UV background + Redshift: 0 # Redshift to use (-1 means time based redshift) + WithMetalCooling: 1 # Enable or not the metal cooling + ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates + ProvideSpecificHeatingRates: 0 # User provide specific heating rates + SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method + OutputMode: 1 # Write in output corresponding primordial chemistry mode + MaxSteps: 1000 + ConvergenceLimit: 1e-2 + EAGLEChemistry: InitMetallicity: 0. InitAbundance_Hydrogen: 0.752 @@ -60,3 +66,6 @@ EAGLEChemistry: InitAbundance_Iron: 0.000 CalciumOverSilicon: 0.0941736 SulphurOverSilicon: 0.6054160 + +GearChemistry: + InitialMetallicity: 0.01295 diff --git a/examples/main.c b/examples/main.c index ca91c9e718cca9f7f8d82fce2f5d237f9f5b98e4..39177560981c8d292ff4e69d871db7515645431c 100644 --- a/examples/main.c +++ b/examples/main.c @@ -124,7 +124,7 @@ int main(int argc, char *argv[]) { /* Structs used by the engine. Declare now to make sure these are always in * scope. */ - struct chemistry_data chemistry; + struct chemistry_global_data chemistry; struct cooling_function_data cooling_func; struct cosmology cosmo; struct external_potential potential; diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 0a81716a3e34bbb305916b1148a9d8199b84fa5a..edf1afe51a1169c672eedc6fc629e044dda0d966 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -179,12 +179,18 @@ LambdaCooling: hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless) cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition -# Cooling with Grackle 2.0 +# Cooling with Grackle 3.0 GrackleCooling: - GrackleCloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) - 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 + CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository) + WithUVbackground: 1 # Enable or not the UV background + Redshift: 0 # Redshift to use (-1 means time based redshift) + WithMetalCooling: 1 # Enable or not the metal cooling + ProvideVolumetricHeatingRates: 0 # (optional) User provide volumetric heating rates + ProvideSpecificHeatingRates: 0 # (optional) User provide specific heating rates + SelfShieldingMethod: 0 # (optional) Grackle (<= 3) or Gear self shielding method + OutputMode: 0 # (optional) Write in output corresponding primordial chemistry mode + MaxSteps: 10000 # (optional) Max number of step when computing the initial composition + ConvergenceLimit: 1e-2 # (optional) Convergence threshold (relative) for initial composition # Parameters related to chemistry models ----------------------------------------------- diff --git a/src/Makefile.am b/src/Makefile.am index 1bd80b669396d4b4597f6149acafe6b1a4165f43..e21cb6539d2b7ae360ea07c92c426e99342bc2ec 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -43,8 +43,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \ common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \ physical_constants.h physical_constants_cgs.h potential.h version.h \ - hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \ - sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \ + hydro_properties.h riemann.h threadpool.h cooling_io.h cooling.h cooling_struct.h \ + sourceterms.h sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \ gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \ chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h @@ -118,10 +118,15 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h potential/isothermal/potential.h potential/disc_patch/potential.h \ potential/sine_wave/potential.h \ cooling/none/cooling.h cooling/none/cooling_struct.h \ + cooling/none/cooling_io.h \ cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \ + cooling/const_du/cooling_io.h \ cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \ + cooling/const_lambda/cooling_io.h \ cooling/grackle/cooling.h cooling/grackle/cooling_struct.h \ + cooling/grackle/cooling_io.h \ cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h \ + cooling/EAGLE/cooling_io.h \ chemistry/none/chemistry.h \ chemistry/none/chemistry_io.h \ chemistry/none/chemistry_struct.h \ diff --git a/src/chemistry.c b/src/chemistry.c index 091a0b363d16507894f3dcbba2ca79860f3cbaa9..44cbea1361d96c4cf1d4d3d21c3c91e5225640a5 100644 --- a/src/chemistry.c +++ b/src/chemistry.c @@ -36,7 +36,7 @@ void chemistry_init(const struct swift_params* parameter_file, const struct unit_system* us, const struct phys_const* phys_const, - struct chemistry_data* data) { + struct chemistry_global_data* data) { chemistry_init_backend(parameter_file, us, phys_const, data); } @@ -46,9 +46,10 @@ void chemistry_init(const struct swift_params* parameter_file, * * Calls chemistry_print_backend for the chosen chemistry model. * - * @brief The #chemistry_data containing information about the current model. + * @brief The #chemistry_global_data containing information about the current + * model. */ -void chemistry_print(const struct chemistry_data* data) { +void chemistry_print(const struct chemistry_global_data* data) { chemistry_print_backend(data); } @@ -58,10 +59,10 @@ void chemistry_print(const struct chemistry_data* data) { * @param chemistry the struct * @param stream the file stream */ -void chemistry_struct_dump(const struct chemistry_data* chemistry, +void chemistry_struct_dump(const struct chemistry_global_data* chemistry, FILE* stream) { - restart_write_blocks((void*)chemistry, sizeof(struct chemistry_data), 1, - stream, "chemistry", "chemistry function"); + restart_write_blocks((void*)chemistry, sizeof(struct chemistry_global_data), + 1, stream, "chemistry", "chemistry function"); } /** @@ -71,8 +72,8 @@ void chemistry_struct_dump(const struct chemistry_data* chemistry, * @param chemistry the struct * @param stream the file stream */ -void chemistry_struct_restore(const struct chemistry_data* chemistry, +void chemistry_struct_restore(const struct chemistry_global_data* chemistry, FILE* stream) { - restart_read_blocks((void*)chemistry, sizeof(struct chemistry_data), 1, + restart_read_blocks((void*)chemistry, sizeof(struct chemistry_global_data), 1, stream, NULL, "chemistry function"); } diff --git a/src/chemistry.h b/src/chemistry.h index a5cbd77efbaab99e98ca5f031512e7e5bef0a613..77b8d5b5fa5314d692fb8f46c8e6229b57119a6c 100644 --- a/src/chemistry.h +++ b/src/chemistry.h @@ -46,14 +46,14 @@ void chemistry_init(const struct swift_params* parameter_file, const struct unit_system* us, const struct phys_const* phys_const, - struct chemistry_data* data); + struct chemistry_global_data* data); -void chemistry_print(const struct chemistry_data* data); +void chemistry_print(const struct chemistry_global_data* data); /* Dump/restore. */ -void chemistry_struct_dump(const struct chemistry_data* chemistry, +void chemistry_struct_dump(const struct chemistry_global_data* chemistry, FILE* stream); -void chemistry_struct_restore(const struct chemistry_data* chemistry, +void chemistry_struct_restore(const struct chemistry_global_data* chemistry, FILE* stream); #endif /* SWIFT_CHEMISTRY_H */ diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h index ae333252426f5fb79cc1fd1cee818754483214f8..fa63c8cc877861d443d6cdb3cc1c8f0c94befefe 100644 --- a/src/chemistry/EAGLE/chemistry.h +++ b/src/chemistry/EAGLE/chemistry.h @@ -57,10 +57,10 @@ chemistry_get_element_name(enum chemistry_element elem) { * the various smooth metallicity tasks * * @param p The particle to act upon - * @param cd #chemistry_data containing chemistry informations. + * @param cd #chemistry_global_data containing chemistry informations. */ __attribute__((always_inline)) INLINE static void chemistry_init_part( - struct part* restrict p, const struct chemistry_data* cd) {} + struct part* restrict p, const struct chemistry_global_data* cd) {} /** * @brief Finishes the smooth metal calculation. @@ -71,11 +71,11 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part( * This function requires the #hydro_end_density to have been called. * * @param p The particle to act upon. - * @param cd #chemistry_data containing chemistry informations. + * @param cd #chemistry_global_data containing chemistry informations. * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void chemistry_end_density( - struct part* restrict p, const struct chemistry_data* cd, + struct part* restrict p, const struct chemistry_global_data* cd, const struct cosmology* cosmo) {} /** @@ -87,8 +87,11 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( * @param data The global chemistry information. */ __attribute__((always_inline)) INLINE static void chemistry_first_init_part( - struct part* restrict p, struct xpart* restrict xp, - const struct chemistry_data* data) { + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct chemistry_global_data* data, struct part* restrict p, + struct xpart* restrict xp) { p->chemistry_data.metal_mass_fraction_total = data->initial_metal_mass_fraction_total; @@ -107,7 +110,7 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_part( */ static INLINE void chemistry_init_backend( const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const, struct chemistry_data* data) { + const struct phys_const* phys_const, struct chemistry_global_data* data) { /* Read the total metallicity */ data->initial_metal_mass_fraction_total = @@ -133,9 +136,11 @@ static INLINE void chemistry_init_backend( /** * @brief Prints the properties of the chemistry model to stdout. * - * @brief The #chemistry_data containing information about the current model. + * @brief The #chemistry_global_data containing information about the current + * model. */ -static INLINE void chemistry_print_backend(const struct chemistry_data* data) { +static INLINE void chemistry_print_backend( + const struct chemistry_global_data* data) { message("Chemistry model is 'EAGLE' tracking %d elements.", chemistry_element_count); diff --git a/src/chemistry/EAGLE/chemistry_struct.h b/src/chemistry/EAGLE/chemistry_struct.h index e76e082d9b92bcde800b084feea44515b9579dc8..9093709e62d0af638ae485bd1154a4537791e84a 100644 --- a/src/chemistry/EAGLE/chemistry_struct.h +++ b/src/chemistry/EAGLE/chemistry_struct.h @@ -38,7 +38,7 @@ enum chemistry_element { /** * @brief Global chemical abundance information in the EAGLE model. */ -struct chemistry_data { +struct chemistry_global_data { /*! Fraction of the particle mass in given elements at the start of the run */ float initial_metal_mass_fraction[chemistry_element_count]; diff --git a/src/chemistry/gear/chemistry.h b/src/chemistry/gear/chemistry.h index e9d3d00febe2438c52728f94d0136fb72ff8dc3b..9653e7eb05005977c39bca20eab961af83f40b5f 100644 --- a/src/chemistry/gear/chemistry.h +++ b/src/chemistry/gear/chemistry.h @@ -29,6 +29,7 @@ #include /* Local includes. */ +#include "chemistry_io.h" #include "chemistry_struct.h" #include "error.h" #include "hydro.h" @@ -38,16 +39,16 @@ #include "units.h" /** - * @brief Return a string containing the name of a given #chemistry_element. + * @brief Compute the metal mass fraction + * + * @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 const char* -chemistry_get_element_name(enum chemistry_element elem) { - - static const char* chemistry_element_names[chemistry_element_count] = { - "Oxygen", "Magnesium", "Sulfur", "Iron", "Zinc", - "Strontium", "Yttrium", "Barium", "Europium"}; - - return chemistry_element_names[elem]; +__attribute__((always_inline)) INLINE static float +chemistry_metal_mass_fraction(const struct part* restrict p, + const struct xpart* restrict xp) { + return p->chemistry_data.Z; } /** @@ -62,16 +63,10 @@ chemistry_get_element_name(enum chemistry_element elem) { */ static INLINE void chemistry_init_backend( const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const, struct chemistry_data* data) {} + const struct phys_const* phys_const, struct chemistry_global_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(const struct chemistry_data* data) { - - message("Chemistry function is 'Gear'."); + /* read parameters */ + chemistry_read_parameters(parameter_file, us, phys_const, data); } /** @@ -81,10 +76,10 @@ static INLINE void chemistry_print_backend(const struct chemistry_data* data) { * the various smooth metallicity tasks * * @param p The particle to act upon - * @param cd #chemistry_data containing chemistry informations. + * @param cd #chemistry_global_data containing chemistry informations. */ __attribute__((always_inline)) INLINE static void chemistry_init_part( - struct part* restrict p, const struct chemistry_data* cd) { + struct part* restrict p, const struct chemistry_global_data* cd) { struct chemistry_part_data* cpd = &p->chemistry_data; @@ -102,11 +97,11 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part( * This function requires the #hydro_end_density to have been called. * * @param p The particle to act upon. - * @param cd #chemistry_data containing chemistry informations. + * @param cd #chemistry_global_data containing chemistry informations. * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void chemistry_end_density( - struct part* restrict p, const struct chemistry_data* cd, + struct part* restrict p, const struct chemistry_global_data* cd, const struct cosmology* cosmo) { /* Some smoothing length multiples. */ @@ -138,8 +133,13 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( * @param data The global chemistry information. */ __attribute__((always_inline)) INLINE static void chemistry_first_init_part( - struct part* restrict p, struct xpart* restrict xp, - const struct chemistry_data* data) { + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct chemistry_global_data* data, struct part* restrict p, + struct xpart* restrict xp) { + + p->chemistry_data.Z = data->initial_metallicity; chemistry_init_part(p, data); } diff --git a/src/chemistry/gear/chemistry_io.h b/src/chemistry/gear/chemistry_io.h index 9b0852fa2d23c8dd3e4d36406b55f4353d23b807..2343b6ba18ada5d04825a12d53a6f9cc7e984ef2 100644 --- a/src/chemistry/gear/chemistry_io.h +++ b/src/chemistry/gear/chemistry_io.h @@ -19,9 +19,38 @@ #ifndef SWIFT_CHEMISTRY_IO_GEAR_H #define SWIFT_CHEMISTRY_IO_GEAR_H -#include "chemistry.h" #include "chemistry_struct.h" +#include "error.h" #include "io_properties.h" +#include "parser.h" +#include "part.h" +#include "physical_constants.h" +#include "units.h" + +/** + * @brief Return a string containing the name of a given #chemistry_element. + */ +__attribute__((always_inline)) INLINE static const char* +chemistry_get_element_name(enum chemistry_element elem) { + + static const char* chemistry_element_names[chemistry_element_count] = { + "Oxygen", "Magnesium", "Sulfur", "Iron", "Zinc", + "Strontium", "Yttrium", "Barium", "Europium"}; + + return chemistry_element_names[elem]; +} + +/** + * @brief Prints the properties of the chemistry model to stdout. + * + * @brief The #chemistry_global_data containing information about the current + * model. + */ +static INLINE void chemistry_print_backend( + const struct chemistry_global_data* data) { + + message("Chemistry function is 'Gear'."); +} /** * @brief Specifies which particle fields to read from a dataset @@ -31,13 +60,25 @@ * * @return Returns the number of fields to read. */ -int chemistry_read_particles(struct part* parts, struct io_props* list) { +__attribute__((always_inline)) INLINE static int chemistry_read_particles( + struct part* parts, struct io_props* list) { /* List what we want to read */ list[0] = io_make_input_field( "ElementAbundance", FLOAT, chemistry_element_count, OPTIONAL, UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass_fraction); - return 1; + list[1] = io_make_input_field("Z", FLOAT, 1, OPTIONAL, UNIT_CONV_NO_UNITS, + parts, chemistry_data.Z); + + return 2; +} + +__attribute__((always_inline)) INLINE static void chemistry_read_parameters( + const struct swift_params* parameter_file, const struct unit_system* us, + const struct phys_const* phys_const, struct chemistry_global_data* data) { + + data->initial_metallicity = parser_get_opt_param_float( + parameter_file, "GearChemistry:InitialMetallicity", -1); } /** @@ -48,18 +89,21 @@ int chemistry_read_particles(struct part* parts, struct io_props* list) { * * @return Returns the number of fields to write. */ -int chemistry_write_particles(const struct part* parts, struct io_props* list) { +__attribute__((always_inline)) INLINE static int chemistry_write_particles( + const struct part* parts, struct io_props* list) { /* List what we want to write */ list[0] = io_make_output_field( "SmoothedElementAbundance", FLOAT, chemistry_element_count, UNIT_CONV_NO_UNITS, parts, chemistry_data.smoothed_metal_mass_fraction); + list[1] = io_make_output_field("Z", FLOAT, 1, UNIT_CONV_NO_UNITS, parts, + chemistry_data.Z); - list[1] = io_make_output_field("ElementAbundance", FLOAT, + list[2] = io_make_output_field("ElementAbundance", FLOAT, chemistry_element_count, UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass_fraction); - return 2; + return 3; } #ifdef HAVE_HDF5 @@ -68,7 +112,8 @@ int chemistry_write_particles(const struct part* parts, struct io_props* list) { * @brief Writes the current model of SPH to the file * @param h_grp The HDF5 group in which to write */ -void chemistry_write_flavour(hid_t h_grp) { +__attribute__((always_inline)) INLINE static void chemistry_write_flavour( + hid_t h_grp) { io_write_attribute_s(h_grp, "Chemistry Model", "GEAR"); for (enum chemistry_element i = chemistry_element_O; diff --git a/src/chemistry/gear/chemistry_struct.h b/src/chemistry/gear/chemistry_struct.h index b0ddee2b37c2d1126523ff7a4283f5c18ae3a7b4..c3443126cf0ac72c76774d9000fe218378cc6663 100644 --- a/src/chemistry/gear/chemistry_struct.h +++ b/src/chemistry/gear/chemistry_struct.h @@ -38,7 +38,11 @@ enum chemistry_element { /** * @brief Global chemical abundance information. */ -struct chemistry_data {}; +struct chemistry_global_data { + + /* Initial metallicity Z */ + float initial_metallicity; +}; /** * @brief Properties of the chemistry function. @@ -50,6 +54,8 @@ struct chemistry_part_data { /*! Smoothed fraction of the particle mass in a given element */ float smoothed_metal_mass_fraction[chemistry_element_count]; + + float Z; }; #endif /* SWIFT_CHEMISTRY_STRUCT_GEAR_H */ diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h index 7fb83339bed6fc2db007e58a5b8372f003186db5..3ca51660ddfeead2b7ad0010979b719e59c4934e 100644 --- a/src/chemistry/none/chemistry.h +++ b/src/chemistry/none/chemistry.h @@ -61,14 +61,16 @@ chemistry_get_element_name(enum chemistry_element elem) { */ static INLINE void chemistry_init_backend( const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const, struct chemistry_data* data) {} + const struct phys_const* phys_const, struct chemistry_global_data* data) {} /** * @brief Prints the properties of the chemistry model to stdout. * - * @brief The #chemistry_data containing information about the current model. + * @brief The #chemistry_global_data containing information about the current + * model. */ -static INLINE void chemistry_print_backend(const struct chemistry_data* data) { +static INLINE void chemistry_print_backend( + const struct chemistry_global_data* data) { message("Chemistry function is 'No chemistry'."); } @@ -81,7 +83,7 @@ static INLINE void chemistry_print_backend(const struct chemistry_data* data) { * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void chemistry_end_density( - struct part* restrict p, const struct chemistry_data* cd, + struct part* restrict p, const struct chemistry_global_data* cd, const struct cosmology* cosmo) {} /** @@ -90,13 +92,19 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( * * Nothing to do here. * + * @param phys_const The physical constant in internal units. + * @param us The unit system. + * @param cosmo The current cosmological model. + * @param data The global chemistry information used for this run. * @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_first_init_part( - const struct part* restrict p, struct xpart* restrict xp, - const struct chemistry_data* data) {} + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct chemistry_global_data* data, const struct part* restrict p, + struct xpart* restrict xp) {} /** * @brief Sets the chemistry properties of the (x-)particles to a valid start @@ -108,6 +116,6 @@ __attribute__((always_inline)) INLINE static void chemistry_first_init_part( * @param data The global chemistry information. */ __attribute__((always_inline)) INLINE static void chemistry_init_part( - struct part* restrict p, const struct chemistry_data* data) {} + struct part* restrict p, const struct chemistry_global_data* data) {} #endif /* SWIFT_CHEMISTRY_NONE_H */ diff --git a/src/chemistry/none/chemistry_struct.h b/src/chemistry/none/chemistry_struct.h index 8ac5f6924c6893717033cd5f144fe873f6f17e24..a80699d66cccd2e290555006d45216ff53d9c99c 100644 --- a/src/chemistry/none/chemistry_struct.h +++ b/src/chemistry/none/chemistry_struct.h @@ -34,7 +34,7 @@ enum chemistry_element { chemistry_element_count = 0 }; * * Nothing here. */ -struct chemistry_data {}; +struct chemistry_global_data {}; /** * @brief Chemistry properties carried by the #part. diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h index 62d1cc9b72681fe96c901185a10efe5523b206e3..bdf3801887256cb97ae1d5b6a3095250764aa822 100644 --- a/src/cooling/EAGLE/cooling.h +++ b/src/cooling/EAGLE/cooling.h @@ -20,10 +20,13 @@ #define SWIFT_COOLING_EAGLE_H /** - * @file src/cooling/none/cooling.h - * @brief Empty infrastructure for the cases without cooling function + * @file src/cooling/EAGLE/cooling.h + * @brief EAGLE cooling function */ +/* Config parameters. */ +#include "../config.h" + /* Some standard headers. */ #include #include @@ -31,25 +34,11 @@ /* Local includes. */ #include "error.h" #include "hydro.h" -#include "io_properties.h" #include "parser.h" #include "part.h" #include "physical_constants.h" #include "units.h" -#ifdef HAVE_HDF5 - -/** - * @brief Writes the current model of SPH to the file - * @param h_grpsph The HDF5 group in which to write - */ -__attribute__((always_inline)) INLINE static void cooling_write_flavour( - hid_t h_grpsph) { - - io_write_attribute_s(h_grpsph, "Cooling Model", "EAGLE"); -} -#endif - /** * @brief Apply the cooling function to a particle. * @@ -95,8 +84,11 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * @param cooling The properties of the cooling function. */ __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) {} + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, struct xpart* restrict xp) {} /** * @brief Returns the total radiated energy by this particle. diff --git a/src/cooling/EAGLE/cooling_io.h b/src/cooling/EAGLE/cooling_io.h new file mode 100644 index 0000000000000000000000000000000000000000..d6ed4f122863be7b591b095b11803a3d2729f694 --- /dev/null +++ b/src/cooling/EAGLE/cooling_io.h @@ -0,0 +1,56 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + ******************************************************************************/ +#ifndef SWIFT_COOLING_EAGLE_IO_H +#define SWIFT_COOLING_EAGLE_IO_H + +/* Config parameters. */ +#include "../config.h" + +/* Local includes */ +#include "io_properties.h" + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +__attribute__((always_inline)) INLINE static void cooling_write_flavour( + hid_t h_grpsph) { + + io_write_attribute_s(h_grpsph, "Cooling Model", "EAGLE"); +} +#endif + +/** + * @brief Specifies which particle fields to write to a dataset + * + * @param parts The particle array. + * @param list The list of i/o properties to write. + * @param cooling The #cooling_function_data + * + * @return Returns the number of fields to write. + */ +__attribute__((always_inline)) INLINE static int cooling_write_particles( + const struct xpart* xparts, struct io_props* list, + const struct cooling_function_data* cooling) { + return 0; +} + +#endif /* SWIFT_COOLING_EAGLE_IO_H */ diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h index 25aa974e47b277d0dfcc21077e5041a49e8168c8..ba8211174919419c37856dc1fcbdaa73b23e319e 100644 --- a/src/cooling/const_du/cooling.h +++ b/src/cooling/const_du/cooling.h @@ -30,33 +30,21 @@ * realistic functions. */ +/* Config parameters. */ +#include "../config.h" + /* Some standard headers. */ #include /* Local includes. */ #include "const.h" -#include "cooling_struct.h" #include "error.h" #include "hydro.h" -#include "io_properties.h" #include "parser.h" #include "part.h" #include "physical_constants.h" #include "units.h" -#ifdef HAVE_HDF5 - -/** - * @brief Writes the current model of SPH to the file - * @param h_grpsph The HDF5 group in which to write - */ -__attribute__((always_inline)) INLINE static void cooling_write_flavour( - hid_t h_grpsph) { - - io_write_attribute_s(h_grpsph, "Cooling Model", "Constant du/dt"); -} -#endif - /** * @brief Apply the cooling function to a particle. * @@ -141,8 +129,11 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * @param cooling The properties of the cooling function. */ __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) { + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, struct xpart* restrict xp) { xp->cooling_data.radiated_energy = 0.f; } diff --git a/src/cooling/const_du/cooling_io.h b/src/cooling/const_du/cooling_io.h new file mode 100644 index 0000000000000000000000000000000000000000..52a943aca86e51665fd1841d7bcb8a100b046ed8 --- /dev/null +++ b/src/cooling/const_du/cooling_io.h @@ -0,0 +1,58 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * Richard Bower (r.g.bower@durham.ac.uk) + * Stefan Arridge (stefan.arridge@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + ******************************************************************************/ +#ifndef SWIFT_COOLING_CONST_DU_IO_H +#define SWIFT_COOLING_CONST_DU_IO_H + +/* Config parameters. */ +#include "../config.h" + +/* Local includes */ +#include "io_properties.h" + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +__attribute__((always_inline)) INLINE static void cooling_write_flavour( + hid_t h_grpsph) { + + io_write_attribute_s(h_grpsph, "Cooling Model", "Constant du/dt"); +} +#endif + +/** + * @brief Specifies which particle fields to write to a dataset + * + * @param parts The particle array. + * @param list The list of i/o properties to write. + * @param cooling The #cooling_function_data + * + * @return Returns the number of fields to write. + */ +__attribute__((always_inline)) INLINE static int cooling_write_particles( + const struct xpart* xparts, struct io_props* list, + const struct cooling_function_data* cooling) { + return 0; +} + +#endif /* SWIFT_COOLING_CONST_DU_IO_H */ diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h index a402b749300e64ef79253764c402d41f5a214d82..43ca7ab75b0bce370d7405e52cea9b54335ae73c 100644 --- a/src/cooling/const_lambda/cooling.h +++ b/src/cooling/const_lambda/cooling.h @@ -23,6 +23,9 @@ #ifndef SWIFT_COOLING_CONST_LAMBDA_H #define SWIFT_COOLING_CONST_LAMBDA_H +/* Config parameters. */ +#include "../config.h" + /* Some standard headers. */ #include #include @@ -31,25 +34,11 @@ #include "const.h" #include "error.h" #include "hydro.h" -#include "io_properties.h" #include "parser.h" #include "part.h" #include "physical_constants.h" #include "units.h" -#ifdef HAVE_HDF5 - -/** - * @brief Writes the current model of SPH to the file - * @param h_grpsph The HDF5 group in which to write - */ -__attribute__((always_inline)) INLINE static void cooling_write_flavour( - hid_t h_grpsph) { - - io_write_attribute_s(h_grpsph, "Cooling Model", "Constant Lambda"); -} -#endif - /** * @brief Calculates du/dt in code units for a particle. * @@ -154,8 +143,11 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * @param cooling The properties of the cooling function. */ __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) { + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, struct xpart* restrict xp) { xp->cooling_data.radiated_energy = 0.f; } diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h new file mode 100644 index 0000000000000000000000000000000000000000..89c9471a291a4a6a5740a8c6c816913cbc6316a0 --- /dev/null +++ b/src/cooling/const_lambda/cooling_io.h @@ -0,0 +1,59 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk) + * Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * Richard Bower (r.g.bower@durham.ac.uk) + * Stefan Arridge (stefan.arridge@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + ******************************************************************************/ +#ifndef SWIFT_COOLING_CONST_LAMBDA_IO_H +#define SWIFT_COOLING_CONST_LAMBDA_IO_H + +/* Config parameters. */ +#include "../config.h" + +/* Local includes */ +#include "io_properties.h" + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +__attribute__((always_inline)) INLINE static void cooling_write_flavour( + hid_t h_grpsph) { + + io_write_attribute_s(h_grpsph, "Cooling Model", "Constant Lambda"); +} +#endif + +/** + * @brief Specifies which particle fields to write to a dataset + * + * @param parts The particle array. + * @param list The list of i/o properties to write. + * @param cooling The #cooling_function_data + * + * @return Returns the number of fields to write. + */ +__attribute__((always_inline)) INLINE static int cooling_write_particles( + const struct xpart* xparts, struct io_props* list, + const struct cooling_function_data* cooling) { + return 0; +} + +#endif /* SWIFT_COOLING_CONST_LAMBDA_IO_H */ diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h index d1da52300121ca56feea4114d74322120bd501ce..dd59e9af1431681a8c5bdc1e5cb0c22053063651 100644 --- a/src/cooling/grackle/cooling.h +++ b/src/cooling/grackle/cooling.h @@ -20,19 +20,24 @@ #define SWIFT_COOLING_GRACKLE_H /** - * @file src/cooling/none/cooling.h - * @brief Empty infrastructure for the cases without cooling function + * @file src/cooling/grackle/cooling.h + * @brief Cooling using the GRACKLE 3.0 library. */ +/* Config parameters. */ +#include "../config.h" + /* Some standard headers. */ #include -#include #include +/* The grackle library itself */ +#include + /* Local includes. */ +#include "chemistry.h" #include "error.h" #include "hydro.h" -#include "io_properties.h" #include "parser.h" #include "part.h" #include "physical_constants.h" @@ -42,18 +47,147 @@ #define GRACKLE_NPART 1 #define GRACKLE_RANK 3 -#ifdef HAVE_HDF5 +/* prototypes */ +static gr_float cooling_time( + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, struct xpart* restrict xp); + +static double cooling_rate(const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, + struct xpart* restrict xp, double dt); /** - * @brief Writes the current model of SPH to the file - * @param h_grpsph The HDF5 group in which to write + * @brief Print the chemical network + * + * @param xp The #xpart to print */ -__attribute__((always_inline)) INLINE static void cooling_write_flavour( - hid_t h_grpsph) { +__attribute__((always_inline)) INLINE static void cooling_print_fractions( + const struct xpart* restrict xp) { + + const struct cooling_xpart_data tmp = xp->cooling_data; +#if COOLING_GRACKLE_MODE > 0 + message("HI %g, HII %g, HeI %g, HeII %g, HeIII %g, e %g", tmp.HI_frac, + tmp.HII_frac, tmp.HeI_frac, tmp.HeII_frac, tmp.HeIII_frac, + tmp.e_frac); +#endif - io_write_attribute_s(h_grpsph, "Cooling Model", "Grackle"); +#if COOLING_GRACKLE_MODE > 1 + message("HM %g, H2I %g, H2II %g", tmp.HM_frac, tmp.H2I_frac, tmp.H2II_frac); +#endif + +#if COOLING_GRACKLE_MODE > 2 + message("DI %g, DII %g, HDI %g", tmp.DI_frac, tmp.DII_frac, tmp.HDI_frac); +#endif + message("Metal: %g", tmp.metal_frac); } + +/** + * @brief Check if the difference of a given field is lower than limit + * + * @param xp First #xpart + * @param old Second #xpart + * @param field The field to check + * @param limit Difference limit + * + * @return 0 if diff > limit + */ +#define cooling_check_field(xp, old, field, limit) \ + ({ \ + float tmp = xp->cooling_data.field - old->cooling_data.field; \ + tmp = fabs(tmp) / xp->cooling_data.field; \ + if (tmp > limit) return 0; \ + }) + +/** + * @brief Check if difference between two particles is lower than a given value + * + * @param xp One of the #xpart + * @param old The other #xpart + * @param limit The difference limit + */ +__attribute__((always_inline)) INLINE static int cooling_converged( + const struct xpart* restrict xp, const struct xpart* restrict old, + const float limit) { + +#if COOLING_GRACKLE_MODE > 0 + cooling_check_field(xp, old, HI_frac, limit); + cooling_check_field(xp, old, HII_frac, limit); + cooling_check_field(xp, old, HeI_frac, limit); + cooling_check_field(xp, old, HeII_frac, limit); + cooling_check_field(xp, old, HeIII_frac, limit); + cooling_check_field(xp, old, e_frac, limit); #endif +#if COOLING_GRACKLE_MODE > 1 + cooling_check_field(xp, old, HM_frac, limit); + cooling_check_field(xp, old, H2I_frac, limit); + cooling_check_field(xp, old, H2II_frac, limit); +#endif + +#if COOLING_GRACKLE_MODE > 2 + cooling_check_field(xp, old, DI_frac, limit); + cooling_check_field(xp, old, DII_frac, limit); + cooling_check_field(xp, old, HDI_frac, limit); +#endif + + return 1; +} + +/** + * @brief Compute equilibrium fraction + * + * @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_compute_equilibrium( + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, struct xpart* restrict xp) { + + /* get temporary data */ + struct part p_tmp = *p; + struct cooling_function_data cooling_tmp = *cooling; + cooling_tmp.chemistry.with_radiative_cooling = 0; + /* need density for computation, therefore quick estimate */ + p_tmp.rho = 0.2387 * p_tmp.mass / pow(p_tmp.h, 3); + + /* compute time step */ + const double alpha = 0.01; + double dt = + fabs(cooling_time(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp)); + cooling_rate(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp, dt); + dt = alpha * + fabs(cooling_time(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp)); + + /* init simple variables */ + int step = 0; + const int max_step = cooling_tmp.max_step; + const float conv_limit = cooling_tmp.convergence_limit; + struct xpart old; + + do { + /* update variables */ + step += 1; + old = *xp; + + /* update chemistry */ + cooling_rate(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp, dt); + } while (step < max_step && !cooling_converged(xp, &old, conv_limit)); + + if (step == max_step) + error( + "A particle element fraction failed to converge." + "You can change 'GrackleCooling:MaxSteps' or " + "'GrackleCooling:ConvergenceLimit' to avoid this problem"); +} /** * @brief Sets the cooling properties of the (x-)particles to a valid start @@ -64,10 +198,46 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour( * @param cooling The properties of the cooling function. */ __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) { + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* cooling, const struct part* restrict p, + struct xpart* restrict xp) { xp->cooling_data.radiated_energy = 0.f; + +#if COOLING_GRACKLE_MODE >= 1 + gr_float zero = 1.e-20; + + /* primordial chemistry >= 1 */ + xp->cooling_data.HI_frac = zero; + xp->cooling_data.HII_frac = grackle_data->HydrogenFractionByMass; + xp->cooling_data.HeI_frac = 1. - grackle_data->HydrogenFractionByMass; + xp->cooling_data.HeII_frac = zero; + xp->cooling_data.HeIII_frac = zero; + xp->cooling_data.e_frac = xp->cooling_data.HII_frac + + 0.25 * xp->cooling_data.HeII_frac + + 0.5 * xp->cooling_data.HeIII_frac; +#endif // MODE >= 1 + +#if COOLING_GRACKLE_MODE >= 2 + /* primordial chemistry >= 2 */ + xp->cooling_data.HM_frac = zero; + xp->cooling_data.H2I_frac = zero; + xp->cooling_data.H2II_frac = zero; +#endif // MODE >= 2 + +#if COOLING_GRACKLE_MODE >= 3 + /* primordial chemistry >= 3 */ + xp->cooling_data.DI_frac = grackle_data->DeuteriumToHydrogenRatio * + grackle_data->HydrogenFractionByMass; + xp->cooling_data.DII_frac = zero; + xp->cooling_data.HDI_frac = zero; +#endif // MODE >= 3 + +#if COOLING_GRACKLE_MODE > 0 + cooling_compute_equilibrium(phys_const, us, cosmo, cooling, p, xp); +#endif } /** @@ -90,17 +260,17 @@ __attribute__((always_inline)) INLINE static void cooling_print_backend( const struct cooling_function_data* cooling) { message("Cooling function is 'Grackle'."); - message("Using Grackle = %i", cooling->chemistry.use_grackle); - message("Chemical network = %i", - cooling->chemistry.primordial_chemistry); - message("Radiative cooling = %i", - cooling->chemistry.with_radiative_cooling); - message("Metal cooling = %i", cooling->chemistry.metal_cooling); - - message("CloudyTable = %s", cooling->cloudy_table); - message("UVbackground = %d", cooling->uv_background); - message("Redshift = %g", cooling->redshift); - message("Density Self Shielding = %g", cooling->density_self_shielding); + message("Using Grackle = %i", cooling->chemistry.use_grackle); + message("Chemical network = %i", cooling->chemistry.primordial_chemistry); + message("CloudyTable = %s", cooling->cloudy_table); + message("Redshift = %g", cooling->redshift); + message("UV background = %d", cooling->with_uv_background); + message("Metal cooling = %i", cooling->chemistry.metal_cooling); + message("Self Shielding = %i", cooling->self_shielding_method); + message("Specific Heating Rates = %i", + cooling->provide_specific_heating_rates); + message("Volumetric Heating Rates = %i", + cooling->provide_volumetric_heating_rates); message("Units:"); message("\tComoving = %i", cooling->units.comoving_coordinates); message("\tLength = %g", cooling->units.length_units); @@ -110,29 +280,221 @@ __attribute__((always_inline)) INLINE static void cooling_print_backend( } /** - * @brief Compute the cooling rate + * @brief copy a single field from the grackle data to a #xpart + * + * @param data The #grackle_field_data + * @param xp The #xpart + * @param rho Particle density + * @param field The field to copy + */ +#define cooling_copy_field_from_grackle(data, xp, rho, field) \ + xp->cooling_data.field##_frac = *data.field##_density / rho; + +/** + * @brief copy a single field from a #xpart to the grackle data + * + * @param data The #grackle_field_data + * @param xp The #xpart + * @param rho Particle density + * @param field The field to copy + */ +#define cooling_copy_field_to_grackle(data, xp, rho, field) \ + gr_float grackle_##field = xp->cooling_data.field##_frac * rho; \ + data.field##_density = &grackle_##field; + +/** + * @brief copy a #xpart to the grackle data + * + * Warning this function creates some variable, therefore the grackle call + * should be in a block that still has the variables. + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#if COOLING_GRACKLE_MODE > 0 +#define cooling_copy_to_grackle1(data, p, xp, rho) \ + cooling_copy_field_to_grackle(data, xp, rho, HI); \ + cooling_copy_field_to_grackle(data, xp, rho, HII); \ + cooling_copy_field_to_grackle(data, xp, rho, HeI); \ + cooling_copy_field_to_grackle(data, xp, rho, HeII); \ + cooling_copy_field_to_grackle(data, xp, rho, HeIII); \ + cooling_copy_field_to_grackle(data, xp, rho, e); +#else +#define cooling_copy_to_grackle1(data, p, xp, rho) \ + data.HI_density = NULL; \ + data.HII_density = NULL; \ + data.HeI_density = NULL; \ + data.HeII_density = NULL; \ + data.HeIII_density = NULL; \ + data.e_density = NULL; +#endif + +/** + * @brief copy a #xpart to the grackle data + * + * Warning this function creates some variable, therefore the grackle call + * should be in a block that still has the variables. + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#if COOLING_GRACKLE_MODE > 1 +#define cooling_copy_to_grackle2(data, p, xp, rho) \ + cooling_copy_field_to_grackle(data, xp, rho, HM); \ + cooling_copy_field_to_grackle(data, xp, rho, H2I); \ + cooling_copy_field_to_grackle(data, xp, rho, H2II); +#else +#define cooling_copy_to_grackle2(data, p, xp, rho) \ + data.HM_density = NULL; \ + data.H2I_density = NULL; \ + data.H2II_density = NULL; +#endif + +/** + * @brief copy a #xpart to the grackle data + * + * Warning this function creates some variable, therefore the grackle call + * should be in a block that still has the variables. + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#if COOLING_GRACKLE_MODE > 2 +#define cooling_copy_to_grackle3(data, p, xp, rho) \ + cooling_copy_field_to_grackle(data, xp, rho, DI); \ + cooling_copy_field_to_grackle(data, xp, rho, DII); \ + cooling_copy_field_to_grackle(data, xp, rho, HDI); +#else +#define cooling_copy_to_grackle3(data, p, xp, rho) \ + data.DI_density = NULL; \ + data.DII_density = NULL; \ + data.HDI_density = NULL; +#endif + +/** + * @brief copy the grackle data to a #xpart + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#if COOLING_GRACKLE_MODE > 0 +#define cooling_copy_from_grackle1(data, p, xp, rho) \ + cooling_copy_field_from_grackle(data, xp, rho, HI); \ + cooling_copy_field_from_grackle(data, xp, rho, HII); \ + cooling_copy_field_from_grackle(data, xp, rho, HeI); \ + cooling_copy_field_from_grackle(data, xp, rho, HeII); \ + cooling_copy_field_from_grackle(data, xp, rho, HeIII); \ + cooling_copy_field_from_grackle(data, xp, rho, e); +#else +#define cooling_copy_from_grackle1(data, p, xp, rho) +#endif + +/** + * @brief copy the grackle data to a #xpart + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#if COOLING_GRACKLE_MODE > 1 +#define cooling_copy_from_grackle2(data, p, xp, rho) \ + cooling_copy_field_from_grackle(data, xp, rho, HM); \ + cooling_copy_field_from_grackle(data, xp, rho, H2I); \ + cooling_copy_field_from_grackle(data, xp, rho, H2II); +#else +#define cooling_copy_from_grackle2(data, p, xp, rho) +#endif + +/** + * @brief copy the grackle data to a #xpart + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#if COOLING_GRACKLE_MODE > 2 +#define cooling_copy_from_grackle3(data, p, xp, rho) \ + cooling_copy_field_from_grackle(data, xp, rho, DI); \ + cooling_copy_field_from_grackle(data, xp, rho, DII); \ + cooling_copy_field_from_grackle(data, xp, rho, HDI); +#else +#define cooling_copy_from_grackle3(data, p, xp, rho) +#endif + +/** + * @brief copy a #xpart to the grackle data + * + * Warning this function creates some variable, therefore the grackle call + * should be in a block that still has the variables. + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#define cooling_copy_to_grackle(data, p, xp, rho) \ + cooling_copy_to_grackle1(data, p, xp, rho); \ + cooling_copy_to_grackle2(data, p, xp, rho); \ + cooling_copy_to_grackle3(data, p, xp, rho); \ + data.volumetric_heating_rate = NULL; \ + data.specific_heating_rate = NULL; \ + data.RT_heating_rate = NULL; \ + data.RT_HI_ionization_rate = NULL; \ + data.RT_HeI_ionization_rate = NULL; \ + data.RT_HeII_ionization_rate = NULL; \ + data.RT_H2_dissociation_rate = NULL; \ + gr_float metal_density = chemistry_metal_mass_fraction(p, xp) * rho; \ + data.metal_density = &metal_density; + +/** + * @brief copy a #xpart to the grackle data + * + * Warning this function creates some variable, therefore the grackle call + * should be in a block that still has the variables. + * + * @param data The #grackle_field_data + * @param p The #part + * @param xp The #xpart + * @param rho Particle density + */ +#define cooling_copy_from_grackle(data, p, xp, rho) \ + cooling_copy_from_grackle1(data, p, xp, rho); \ + cooling_copy_from_grackle2(data, p, xp, rho); \ + cooling_copy_from_grackle3(data, p, xp, rho); + +/** + * @brief Compute the cooling rate and update the particle chemistry data * * @param phys_const The physical constants in internal units. * @param us The internal system of units. * @param cooling The #cooling_function_data used in the run. * @param p Pointer to the particle data. + * @param xp Pointer to the particle extra data * @param dt The time-step of this particle. * * @return du / dt */ -__attribute__((always_inline)) INLINE static double cooling_rate( +__attribute__((always_inline)) INLINE static gr_float cooling_rate( const struct phys_const* restrict phys_const, const struct unit_system* restrict us, const struct cosmology* restrict cosmo, const struct cooling_function_data* restrict cooling, - struct part* restrict p, float dt) { - - if (cooling->chemistry.primordial_chemistry > 1) error("Not implemented"); + const struct part* restrict p, struct xpart* restrict xp, double dt) { /* set current time */ code_units units = cooling->units; if (cooling->redshift == -1) - error("TODO time dependant redshift"); + units.a_value = cosmo->a; else units.a_value = 1. / (1. + cooling->redshift); @@ -145,6 +507,7 @@ __attribute__((always_inline)) INLINE static double cooling_rate( int grid_start[GRACKLE_RANK] = {0, 0, 0}; int grid_end[GRACKLE_RANK] = {GRACKLE_NPART - 1, 0, 0}; + data.grid_dx = 0.; data.grid_rank = GRACKLE_RANK; data.grid_dimension = grid_dimension; data.grid_start = grid_start; @@ -154,70 +517,114 @@ __attribute__((always_inline)) INLINE static double cooling_rate( gr_float density = hydro_get_physical_density(p, cosmo); const double energy_before = hydro_get_physical_internal_energy(p, cosmo); gr_float energy = energy_before; - gr_float vx = 0; - gr_float vy = 0; - gr_float vz = 0; + /* initialize density */ data.density = &density; + + /* initialize energy */ data.internal_energy = &energy; - data.x_velocity = &vx; - data.y_velocity = &vy; - data.z_velocity = &vz; - - /* /\* primordial chemistry >= 1 *\/ */ - /* gr_float HI_density = density; */ - /* gr_float HII_density = 0.; */ - /* gr_float HeI_density = 0.; */ - /* gr_float HeII_density = 0.; */ - /* gr_float HeIII_density = 0.; */ - /* gr_float e_density = 0.; */ - - /* data.HI_density = &HI_density; */ - /* data.HII_density = &HII_density; */ - /* data.HeI_density = &HeI_density; */ - /* data.HeII_density = &HeII_density; */ - /* data.HeIII_density = &HeIII_density; */ - /* data.e_density = &e_density; */ - - /* /\* primordial chemistry >= 2 *\/ */ - /* gr_float HM_density = 0.; */ - /* gr_float H2I_density = 0.; */ - /* gr_float H2II_density = 0.; */ - - /* data.HM_density = &HM_density; */ - /* data.H2I_density = &H2I_density; */ - /* data.H2II_density = &H2II_density; */ - - /* /\* primordial chemistry >= 3 *\/ */ - /* gr_float DI_density = 0.; */ - /* gr_float DII_density = 0.; */ - /* gr_float HDI_density = 0.; */ - - /* data.DI_density = &DI_density; */ - /* data.DII_density = &DII_density; */ - /* data.HDI_density = &HDI_density; */ - - /* metal cooling = 1 */ - gr_float metal_density = density * grackle_data->SolarMetalFractionByMass; - data.metal_density = &metal_density; + /* grackle 3.0 doc: "Currently not used" */ + data.x_velocity = NULL; + data.y_velocity = NULL; + data.z_velocity = NULL; + + /* copy to grackle structure */ + cooling_copy_to_grackle(data, p, xp, density); + + /* solve chemistry */ + chemistry_data chemistry_grackle = cooling->chemistry; + chemistry_data_storage my_rates = grackle_rates; + int error_code = _solve_chemistry( + &chemistry_grackle, &my_rates, &units, dt, data.grid_dx, data.grid_rank, + data.grid_dimension, data.grid_start, data.grid_end, data.density, + data.internal_energy, data.x_velocity, data.y_velocity, data.z_velocity, + data.HI_density, data.HII_density, data.HM_density, data.HeI_density, + data.HeII_density, data.HeIII_density, data.H2I_density, + data.H2II_density, data.DI_density, data.DII_density, data.HDI_density, + data.e_density, data.metal_density, data.volumetric_heating_rate, + data.specific_heating_rate, data.RT_heating_rate, + data.RT_HI_ionization_rate, data.RT_HeI_ionization_rate, + data.RT_HeII_ionization_rate, data.RT_H2_dissociation_rate, NULL); + if (error_code == 0) error("Error in solve_chemistry."); + // if (solve_chemistry(&units, &data, dt) == 0) { + // error("Error in solve_chemistry."); + //} + + /* copy from grackle data to particle */ + cooling_copy_from_grackle(data, p, xp, density); + + /* compute rate */ + return (energy - energy_before) / dt; +} - /* /\* volumetric heating rate *\/ */ - /* gr_float volumetric_heating_rate = 0.; */ +/** + * @brief Compute the cooling time + * + * @param cooling The #cooling_function_data used in the run. + * @param p Pointer to the particle data. + * @param xp Pointer to the particle extra data + * + * @return cooling time + */ +__attribute__((always_inline)) INLINE static gr_float cooling_time( + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* restrict cooling, + const struct part* restrict p, struct xpart* restrict xp) { - /* data.volumetric_heating_rate = &volumetric_heating_rate; */ + /* set current time */ + code_units units = cooling->units; + if (cooling->redshift == -1) + error("TODO time dependant redshift"); + else + units.a_value = 1. / (1. + cooling->redshift); - /* /\* specific heating rate *\/ */ - /* gr_float specific_heating_rate = 0.; */ + /* initialize data */ + grackle_field_data data; + + /* set values */ + /* grid */ + int grid_dimension[GRACKLE_RANK] = {GRACKLE_NPART, 1, 1}; + int grid_start[GRACKLE_RANK] = {0, 0, 0}; + int grid_end[GRACKLE_RANK] = {GRACKLE_NPART - 1, 0, 0}; + + data.grid_rank = GRACKLE_RANK; + data.grid_dimension = grid_dimension; + data.grid_start = grid_start; + data.grid_end = grid_end; - /* data.specific_heating_rate = &specific_heating_rate; */ + /* general particle data */ + const gr_float energy_before = hydro_get_physical_internal_energy(p, cosmo); + gr_float density = hydro_get_physical_density(p, cosmo); + gr_float energy = energy_before; + + /* initialize density */ + data.density = &density; + + /* initialize energy */ + data.internal_energy = &energy; + + /* grackle 3.0 doc: "Currently not used" */ + data.x_velocity = NULL; + data.y_velocity = NULL; + data.z_velocity = NULL; + + /* copy data from particle to grackle data */ + cooling_copy_to_grackle(data, p, xp, density); - /* solve chemistry with table */ - if (solve_chemistry(&units, &data, dt) == 0) { - error("Error in solve_chemistry."); + /* Compute cooling time */ + gr_float cooling_time; + if (calculate_cooling_time(&units, &data, &cooling_time) == 0) { + error("Error in calculate_cooling_time."); } - return (energy - energy_before) / dt; + /* copy from grackle data to particle */ + cooling_copy_from_grackle(data, p, xp, density); + + /* compute rate */ + return cooling_time; } /** @@ -235,7 +642,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( const struct unit_system* restrict us, const struct cosmology* restrict cosmo, const struct cooling_function_data* restrict cooling, - struct part* restrict p, struct xpart* restrict xp, float dt) { + struct part* restrict p, struct xpart* restrict xp, double dt) { if (dt == 0.) return; @@ -243,7 +650,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( const float hydro_du_dt = hydro_get_internal_energy_dt(p); /* compute cooling rate */ - const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p, dt); + const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p, xp, dt); /* record energy lost */ xp->cooling_data.radiated_energy += -du_dt * dt * hydro_get_mass(p); @@ -273,37 +680,15 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( } /** - * @brief Initialises the cooling properties. + * @brief Initialises the cooling unit system. * - * @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 cooling The cooling properties to initialize */ -__attribute__((always_inline)) INLINE static void cooling_init_backend( - const struct swift_params* parameter_file, const struct unit_system* us, - const struct phys_const* phys_const, - struct cooling_function_data* cooling) { - - /* read parameters */ - parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable", - cooling->cloudy_table); - cooling->uv_background = - parser_get_param_int(parameter_file, "GrackleCooling:UVbackground"); - - cooling->redshift = - parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift"); - - cooling->density_self_shielding = parser_get_param_double( - parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold"); - -#ifdef SWIFT_DEBUG_CHECKS - /* enable verbose for grackle */ - grackle_verbose = 1; -#endif +__attribute__((always_inline)) INLINE static void cooling_init_units( + const struct unit_system* us, struct cooling_function_data* cooling) { - /* Set up the units system. - These are conversions from code units to cgs. */ + /* These are conversions from code units to cgs. */ /* first cosmo */ cooling->units.a_units = 1.0; // units for the expansion factor (1/1+zi) @@ -321,6 +706,20 @@ __attribute__((always_inline)) INLINE static void cooling_init_backend( cooling->units.velocity_units = cooling->units.a_units * cooling->units.length_units / cooling->units.time_units; +} + +/** + * @brief Initialises Grackle. + * + * @param cooling The cooling properties to initialize + */ +__attribute__((always_inline)) INLINE static void cooling_init_grackle( + struct cooling_function_data* cooling) { + +#ifdef SWIFT_DEBUG_CHECKS + /* enable verbose for grackle */ + grackle_verbose = 1; +#endif chemistry_data* chemistry = &cooling->chemistry; @@ -332,38 +731,60 @@ __attribute__((always_inline)) INLINE static void cooling_init_backend( // Set parameter values for chemistry. chemistry->use_grackle = 1; chemistry->with_radiative_cooling = 1; + /* molecular network with H, He, D From Cloudy table */ - chemistry->primordial_chemistry = 0; - chemistry->metal_cooling = 1; // metal cooling on - chemistry->UVbackground = cooling->uv_background; + chemistry->primordial_chemistry = COOLING_GRACKLE_MODE; + chemistry->metal_cooling = cooling->with_metal_cooling; + chemistry->UVbackground = cooling->with_uv_background; chemistry->grackle_data_file = cooling->cloudy_table; - chemistry->use_radiative_transfer = 0; - chemistry->use_volumetric_heating_rate = 0; - chemistry->use_specific_heating_rate = 0; + + /* radiative transfer */ + chemistry->use_radiative_transfer = cooling->provide_specific_heating_rates || + cooling->provide_volumetric_heating_rates; + chemistry->use_volumetric_heating_rate = + cooling->provide_volumetric_heating_rates; + chemistry->use_specific_heating_rate = + cooling->provide_specific_heating_rates; + + if (cooling->provide_specific_heating_rates && + cooling->provide_volumetric_heating_rates) + message( + "WARNING: You should specified either the specific or the volumetric " + "heating rates, not both"); + + /* self shielding */ + chemistry->self_shielding_method = cooling->self_shielding_method; /* Initialize the chemistry object. */ if (initialize_chemistry_data(&cooling->units) == 0) { error("Error in initialize_chemistry_data."); } +} + +/** + * @brief Initialises the cooling properties. + * + * @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 cooling The cooling properties to initialize + */ +__attribute__((always_inline)) INLINE static void cooling_init_backend( + const struct swift_params* parameter_file, const struct unit_system* us, + const struct phys_const* phys_const, + struct cooling_function_data* cooling) { -#ifdef SWIFT_DEBUG_CHECKS if (GRACKLE_NPART != 1) error("Grackle with multiple particles not implemented"); - float threshold = cooling->density_self_shielding; - threshold /= phys_const->const_proton_mass; - threshold /= pow(us->UnitLength_in_cgs, 3); + /* read parameters */ + cooling_read_parameters(parameter_file, cooling); - message("***************************************"); - message("initializing grackle cooling function"); - message(""); - cooling_print_backend(cooling); - message("Density Self Shielding = %g atom/cm3", threshold); + /* Set up the units system. */ + cooling_init_units(us, cooling); - message(""); - message("***************************************"); -#endif + cooling_init_grackle(cooling); } #endif /* SWIFT_COOLING_GRACKLE_H */ diff --git a/src/cooling/grackle/cooling_io.h b/src/cooling/grackle/cooling_io.h new file mode 100644 index 0000000000000000000000000000000000000000..5a6edb8f1c559a7b495351e256559f251b97c1cf --- /dev/null +++ b/src/cooling/grackle/cooling_io.h @@ -0,0 +1,169 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + ******************************************************************************/ +#ifndef SWIFT_COOLING_GRACKLE_IO_H +#define SWIFT_COOLING_GRACKLE_IO_H + +/* Config parameters. */ +#include "../config.h" + +/* Local includes */ +#include "cooling_struct.h" +#include "io_properties.h" + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +__attribute__((always_inline)) INLINE static void cooling_write_flavour( + hid_t h_grpsph) { + +#if COOLING_GRACKLE_MODE == 0 + io_write_attribute_s(h_grpsph, "Cooling Model", "Grackle"); +#elif COOLING_GRACKLE_MODE == 1 + io_write_attribute_s(h_grpsph, "Cooling Model", "Grackle1"); +#elif COOLING_GRACKLE_MODE == 2 + io_write_attribute_s(h_grpsph, "Cooling Model", "Grackle2"); +#elif COOLING_GRACKLE_MODE == 3 + io_write_attribute_s(h_grpsph, "Cooling Model", "Grackle3"); +#else + error("This function should be called only with one of the Grackle cooling."); +#endif +} +#endif + +/** + * @brief Specifies which particle fields to write to a dataset + * + * @param parts The particle array. + * @param list The list of i/o properties to write. + * @param cooling The #cooling_function_data + * + * @return Returns the number of fields to write. + */ +__attribute__((always_inline)) INLINE static int cooling_write_particles( + const struct xpart* xparts, struct io_props* list, + const struct cooling_function_data* cooling) { + + int num = 0; + + if (cooling->output_mode == 0) return num; + +#if COOLING_GRACKLE_MODE >= 1 + /* List what we want to write */ + list[0] = io_make_output_field("HI", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HI_frac); + + list[1] = io_make_output_field("HII", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HII_frac); + + list[2] = io_make_output_field("HeI", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HeI_frac); + + list[3] = io_make_output_field("HeII", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HeII_frac); + + list[4] = io_make_output_field("HeIII", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HeIII_frac); + + list[5] = io_make_output_field("e", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.e_frac); + + num += 6; +#endif + + if (cooling->output_mode == 1) return num; + +#if COOLING_GRACKLE_MODE >= 2 + list += num; + + list[0] = io_make_output_field("HM", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HM_frac); + + list[1] = io_make_output_field("H2I", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.H2I_frac); + + list[2] = io_make_output_field("H2II", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.H2II_frac); + + num += 3; +#endif + + if (cooling->output_mode == 2) return num; + +#if COOLING_GRACKLE_MODE >= 3 + list += num; + + list[0] = io_make_output_field("DI", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.DI_frac); + + list[1] = io_make_output_field("DII", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.DII_frac); + + list[2] = io_make_output_field("HDI", FLOAT, 1, UNIT_CONV_NO_UNITS, xparts, + cooling_data.HDI_frac); + + num += 3; +#endif + + return num; +} + +/** + * @brief Parser the parameter file and initialize the #cooling_function_data + * + * @param parameter_file The parser parameter file + * @param cooling The cooling properties to initialize + */ +__attribute__((always_inline)) INLINE static void cooling_read_parameters( + const struct swift_params* parameter_file, + struct cooling_function_data* cooling) { + + parser_get_param_string(parameter_file, "GrackleCooling:CloudyTable", + cooling->cloudy_table); + cooling->with_uv_background = + parser_get_param_int(parameter_file, "GrackleCooling:WithUVbackground"); + + cooling->redshift = + parser_get_param_double(parameter_file, "GrackleCooling:Redshift"); + + cooling->with_metal_cooling = + parser_get_param_int(parameter_file, "GrackleCooling:WithMetalCooling"); + + cooling->provide_volumetric_heating_rates = parser_get_opt_param_int( + parameter_file, "GrackleCooling:ProvideVolumetricHeatingRates", 0); + + cooling->provide_specific_heating_rates = parser_get_opt_param_int( + parameter_file, "GrackleCooling:ProvideSpecificHeatingRates", 0); + + cooling->self_shielding_method = parser_get_opt_param_int( + parameter_file, "GrackleCooling:SelfShieldingMethod", 0); + + cooling->output_mode = + parser_get_opt_param_int(parameter_file, "GrackleCooling:OutputMode", 0); + + cooling->max_step = parser_get_opt_param_int( + parameter_file, "GrackleCooling:MaxSteps", 10000); + + cooling->convergence_limit = parser_get_opt_param_double( + parameter_file, "GrackleCooling:ConvergenceLimit", 1e-2); +} + +#endif /* SWIFT_COOLING_GRACKLE_IO_H */ diff --git a/src/cooling/grackle/cooling_struct.h b/src/cooling/grackle/cooling_struct.h index 9ad1470922fdcbf892014f988389eccc1a4e9014..b714690ce4688268723748b29506e458cccc4be9 100644 --- a/src/cooling/grackle/cooling_struct.h +++ b/src/cooling/grackle/cooling_struct.h @@ -22,6 +22,8 @@ /* include grackle */ #include +#include "../config.h" + /** * @file src/cooling/none/cooling_struct.h * @brief Empty infrastructure for the cases without cooling function @@ -36,19 +38,40 @@ struct cooling_function_data { char cloudy_table[200]; /* Enable/Disable UV backgroud */ - int uv_background; + int with_uv_background; /* Redshift to use for the UV backgroud (-1 to use cosmological one) */ double redshift; - /* Density Threshold for the shielding */ - double density_self_shielding; - /* unit system */ code_units units; /* grackle chemistry data */ chemistry_data chemistry; + + /* Enable/Disable metal cooling */ + int with_metal_cooling; + + /* User provide volumetric heating rates */ + int provide_volumetric_heating_rates; + + /* User provide specific heating rates */ + int provide_specific_heating_rates; + + /* Self shielding method (<= 3) means grackle modes */ + int self_shielding_method; + + /* Output mode (correspond to primordial chemistry mode */ + int output_mode; + + /* convergence limit for first init */ + float convergence_limit; + + /* number of step max for first init */ + int max_step; + + /* over relaxation parameter */ + float omega; }; /** @@ -58,6 +81,36 @@ struct cooling_xpart_data { /*! Energy radiated away by this particle since the start of the run */ float radiated_energy; + +/* here all fractions are mass fraction */ +#if COOLING_GRACKLE_MODE >= 1 + /* primordial chemistry >= 1 */ + float HI_frac; + float HII_frac; + float HeI_frac; + float HeII_frac; + float HeIII_frac; + float e_frac; + +#if COOLING_GRACKLE_MODE >= 2 + /* primordial chemistry >= 2 */ + float HM_frac; + float H2I_frac; + float H2II_frac; + +#if COOLING_GRACKLE_MODE >= 3 + /* primordial chemistry >= 3 */ + float DI_frac; + float DII_frac; + float HDI_frac; +#endif // MODE >= 3 + +#endif // MODE >= 2 + +#endif // MODE >= 1 + + /* metal cooling = 1 */ + float metal_frac; }; #endif /* SWIFT_COOLING_STRUCT_NONE_H */ diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h index a1cc6491115a38d6971a57603fef413c35b52027..5081c7cbe6c4b5168da082ead80687226f9d0c16 100644 --- a/src/cooling/none/cooling.h +++ b/src/cooling/none/cooling.h @@ -31,25 +31,11 @@ /* Local includes. */ #include "error.h" #include "hydro.h" -#include "io_properties.h" #include "parser.h" #include "part.h" #include "physical_constants.h" #include "units.h" -#ifdef HAVE_HDF5 - -/** - * @brief Writes the current model of SPH to the file - * @param h_grpsph The HDF5 group in which to write - */ -__attribute__((always_inline)) INLINE static void cooling_write_flavour( - hid_t h_grpsph) { - - io_write_attribute_s(h_grpsph, "Cooling Model", "None"); -} -#endif - /** * @brief Apply the cooling function to a particle. * @@ -96,13 +82,19 @@ __attribute__((always_inline)) INLINE static float cooling_timestep( * * Nothing to do here. * + * @param phys_const The physical constant in internal units. + * @param us The unit system. + * @param cosmo The current cosmological model. + * @param data The properties of the cooling function. * @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_first_init_part( - const struct part* restrict p, struct xpart* restrict xp, - const struct cooling_function_data* cooling) {} + const struct phys_const* restrict phys_const, + const struct unit_system* restrict us, + const struct cosmology* restrict cosmo, + const struct cooling_function_data* data, const struct part* restrict p, + struct xpart* restrict xp) {} /** * @brief Returns the total radiated energy by this particle. diff --git a/src/cooling/none/cooling_io.h b/src/cooling/none/cooling_io.h new file mode 100644 index 0000000000000000000000000000000000000000..e4c84f506bcd31ff95ededb5be889fbf9a27261b --- /dev/null +++ b/src/cooling/none/cooling_io.h @@ -0,0 +1,56 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + ******************************************************************************/ +#ifndef SWIFT_COOLING_NONE_IO_H +#define SWIFT_COOLING_NONE_IO_H + +/* Config parameters. */ +#include "../config.h" + +/* Local includes */ +#include "io_properties.h" + +#ifdef HAVE_HDF5 + +/** + * @brief Writes the current model of SPH to the file + * @param h_grpsph The HDF5 group in which to write + */ +__attribute__((always_inline)) INLINE static void cooling_write_flavour( + hid_t h_grpsph) { + + io_write_attribute_s(h_grpsph, "Cooling Model", "None"); +} +#endif + +/** + * @brief Specifies which particle fields to write to a dataset + * + * @param xparts The extended particle array. + * @param list The list of i/o properties to write. + * @param cooling The #cooling_function_data + * + * @return Returns the number of fields to write. + */ +__attribute__((always_inline)) INLINE static int cooling_write_particles( + const struct xpart* xparts, struct io_props* list, + const struct cooling_function_data* cooling) { + return 0; +} + +#endif /* SWIFT_COOLING_NONE_IO_H */ diff --git a/src/cooling_io.h b/src/cooling_io.h new file mode 100644 index 0000000000000000000000000000000000000000..88eeae2cabdaa8a0909977b84a7dbcf03145d988 --- /dev/null +++ b/src/cooling_io.h @@ -0,0 +1,40 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published + * by the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + ******************************************************************************/ +#ifndef SWIFT_COOLING_IO_H +#define SWIFT_COOLING_IO_H + +/* Config parameters. */ +#include "../config.h" + +/* Import the i/o routines of the right cooling definition */ +#if defined(COOLING_NONE) +#include "./cooling/none/cooling_io.h" +#elif defined(COOLING_CONST_DU) +#include "./cooling/const_du/cooling_io.h" +#elif defined(COOLING_CONST_LAMBDA) +#include "./cooling/const_lambda/cooling_io.h" +#elif defined(COOLING_GRACKLE) +#include "./cooling/grackle/cooling_io.h" +#elif defined(COOLING_EAGLE) +#include "./cooling/EAGLE/cooling_io.h" +#else +#error "Invalid choice of cooling function." +#endif + +#endif /* SWIFT_COOLING_IO_H */ diff --git a/src/engine.c b/src/engine.c index 6f8c0e9937ee40912964b88a3be69cacf060a0ae..2bdad767aeb198d3c4c9d0dc6113a7580445022e 100644 --- a/src/engine.c +++ b/src/engine.c @@ -5202,15 +5202,18 @@ void engine_unpin() { * @param chemistry The chemistry information. * @param sourceterms The properties of the source terms function. */ -void engine_init( - struct engine *e, struct space *s, const struct swift_params *params, - long long Ngas, long long Ndm, int policy, int verbose, - struct repartition *reparttype, const struct unit_system *internal_units, - const struct phys_const *physical_constants, struct cosmology *cosmo, - const struct hydro_props *hydro, 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_init(struct engine *e, struct space *s, + const struct swift_params *params, long long Ngas, + long long Ndm, int policy, int verbose, + struct repartition *reparttype, + const struct unit_system *internal_units, + const struct phys_const *physical_constants, + struct cosmology *cosmo, const struct hydro_props *hydro, + struct gravity_props *gravity, + const struct external_potential *potential, + const struct cooling_function_data *cooling_func, + const struct chemistry_global_data *chemistry, + struct sourceterms *sourceterms) { /* Clean-up everything */ bzero(e, sizeof(struct engine)); @@ -5973,8 +5976,9 @@ void engine_struct_restore(struct engine *e, FILE *stream) { cooling_struct_restore(cooling_func, stream); e->cooling_func = cooling_func; - struct chemistry_data *chemistry = - (struct chemistry_data *)malloc(sizeof(struct chemistry_data)); + struct chemistry_global_data *chemistry = + (struct chemistry_global_data *)malloc( + sizeof(struct chemistry_global_data)); chemistry_struct_restore(chemistry, stream); e->chemistry = chemistry; diff --git a/src/engine.h b/src/engine.h index 4c37aab593fdffc5f68ea16cb83b29ad782debf0..29a064876317d7e1a7ff1cba5e3cb740b6e171a9 100644 --- a/src/engine.h +++ b/src/engine.h @@ -287,7 +287,7 @@ struct engine { const struct cooling_function_data *cooling_func; /* Properties of the chemistry model */ - const struct chemistry_data *chemistry; + const struct chemistry_global_data *chemistry; /* Properties of source terms */ struct sourceterms *sourceterms; @@ -330,15 +330,18 @@ void engine_drift_top_multipoles(struct engine *e); void engine_reconstruct_multipoles(struct engine *e); void engine_print_stats(struct engine *e); void engine_dump_snapshot(struct engine *e); -void engine_init( - struct engine *e, struct space *s, const struct swift_params *params, - long long Ngas, long long Ndm, int policy, int verbose, - struct repartition *reparttype, const struct unit_system *internal_units, - const struct phys_const *physical_constants, struct cosmology *cosmo, - const struct hydro_props *hydro, 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_init(struct engine *e, struct space *s, + const struct swift_params *params, long long Ngas, + long long Ndm, int policy, int verbose, + struct repartition *reparttype, + const struct unit_system *internal_units, + const struct phys_const *physical_constants, + struct cosmology *cosmo, const struct hydro_props *hydro, + struct gravity_props *gravity, + const struct external_potential *potential, + const struct cooling_function_data *cooling_func, + const struct chemistry_global_data *chemistry, + struct sourceterms *sourceterms); void engine_config(int restart, struct engine *e, const struct swift_params *params, int nr_nodes, int nodeID, int nr_threads, int with_aff, int verbose, diff --git a/src/io_properties.h b/src/io_properties.h index 55e128e890b75db2d01911bd3c3b46388b0c6597..d0005531afa3a0cf4de0302e86dd362e0abc11ab 100644 --- a/src/io_properties.h +++ b/src/io_properties.h @@ -26,6 +26,9 @@ #include "common_io.h" #include "inline.h" +/* Standard includes. */ +#include + /** * @brief The two sorts of data present in the GADGET IC files: compulsory to * start a run or optional. diff --git a/src/parallel_io.c b/src/parallel_io.c index 51c88e95396eb6af1ab7f7eb2d1ed48eada59e0d..b5741fb747a0cd8a2dd9ff1645d4a22ffa780869 100644 --- a/src/parallel_io.c +++ b/src/parallel_io.c @@ -38,7 +38,7 @@ /* Local includes. */ #include "chemistry_io.h" #include "common_io.h" -#include "cooling.h" +#include "cooling_io.h" #include "dimension.h" #include "engine.h" #include "error.h" @@ -860,6 +860,8 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6], int periodic = e->s->periodic; int numFiles = 1; + const struct cooling_function_data* cooling = e->cooling_func; + /* First time, we need to create the XMF file */ if (e->snapshotOutputCount == 0) xmf_create_file(baseName); @@ -1243,6 +1245,8 @@ void write_output_parallel(struct engine* e, const char* baseName, Nparticles = Ngas; hydro_write_particles(parts, xparts, list, &num_fields); num_fields += chemistry_write_particles(parts, list + num_fields); + num_fields += + cooling_write_particles(xparts, list + num_fields, cooling); break; case swift_type_dark_matter: diff --git a/src/runner.c b/src/runner.c index 35996a3ba9a760d5f019ede35abc207bbdcbe1b3..51801a779481a8ce0500d8c5f3f0b20d7ed79faf 100644 --- a/src/runner.c +++ b/src/runner.c @@ -656,7 +656,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) { const struct space *s = e->s; const struct hydro_space *hs = &s->hs; const struct cosmology *cosmo = e->cosmology; - const struct chemistry_data *chemistry = e->chemistry; + const struct chemistry_global_data *chemistry = e->chemistry; const float hydro_h_max = e->hydro_properties->h_max; const float eps = e->hydro_properties->h_tolerance; const float hydro_eta_dim = diff --git a/src/serial_io.c b/src/serial_io.c index dd623f946ce6ec32415586e5048979de3adb58fa..fa48d0ac96e58f9dfc4e9d8c513197b6e91800e4 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -38,7 +38,7 @@ /* Local includes. */ #include "chemistry_io.h" #include "common_io.h" -#include "cooling.h" +#include "cooling_io.h" #include "dimension.h" #include "engine.h" #include "error.h" @@ -744,6 +744,7 @@ void write_output_serial(struct engine* e, const char* baseName, const struct gpart* gparts = e->s->gparts; struct gpart* dmparts = NULL; const struct spart* sparts = e->s->sparts; + const struct cooling_function_data* cooling = e->cooling_func; FILE* xmfFile = 0; /* Number of unassociated gparts */ @@ -992,6 +993,8 @@ void write_output_serial(struct engine* e, const char* baseName, Nparticles = Ngas; hydro_write_particles(parts, xparts, list, &num_fields); num_fields += chemistry_write_particles(parts, list + num_fields); + num_fields += + cooling_write_particles(xparts, list + num_fields, cooling); break; case swift_type_dark_matter: diff --git a/src/single_io.c b/src/single_io.c index 2170bcffc4ce3ab21f1edd168d1dc37b2b4af963..73c7f5d946544ed61ebfae01a6ccd4708a8a5f3c 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -37,7 +37,7 @@ /* Local includes. */ #include "chemistry_io.h" #include "common_io.h" -#include "cooling.h" +#include "cooling_io.h" #include "dimension.h" #include "engine.h" #include "error.h" @@ -612,6 +612,7 @@ void write_output_single(struct engine* e, const char* baseName, const struct gpart* gparts = e->s->gparts; struct gpart* dmparts = NULL; const struct spart* sparts = e->s->sparts; + const struct cooling_function_data* cooling = e->cooling_func; /* Number of unassociated gparts */ const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0; @@ -811,6 +812,8 @@ void write_output_single(struct engine* e, const char* baseName, N = Ngas; hydro_write_particles(parts, xparts, list, &num_fields); num_fields += chemistry_write_particles(parts, list + num_fields); + num_fields += + cooling_write_particles(xparts, list + num_fields, cooling); break; case swift_type_dark_matter: diff --git a/src/space.c b/src/space.c index 4f2213793b270cde1c29281bbf411fcfd97c1b4c..b6aef6ea7db1b7bfe0ecbe6a39157f9be5994b50 100644 --- a/src/space.c +++ b/src/space.c @@ -2317,7 +2317,7 @@ void space_synchronize_particle_positions(struct space *s) { * Calls chemistry_first_init_part() on all the particles */ void space_first_init_parts(struct space *s, - const struct chemistry_data *chemistry, + const struct chemistry_global_data *chemistry, const struct cooling_function_data *cool_func) { const size_t nr_parts = s->nr_parts; @@ -2325,6 +2325,8 @@ void space_first_init_parts(struct space *s, struct xpart *restrict xp = s->xparts; const struct cosmology *cosmo = s->e->cosmology; + const struct phys_const *phys_const = s->e->physical_constants; + const struct unit_system *us = s->e->internal_units; const float a_factor_vel = cosmo->a * cosmo->a; const struct hydro_props *hydro_props = s->e->hydro_properties; @@ -2355,10 +2357,10 @@ void space_first_init_parts(struct space *s, if (u_min > 0.f) hydro_set_init_internal_energy(&p[i], u_min); /* Also initialise the chemistry */ - chemistry_first_init_part(&p[i], &xp[i], chemistry); + chemistry_first_init_part(phys_const, us, cosmo, chemistry, &p[i], &xp[i]); /* And the cooling */ - cooling_first_init_part(&p[i], &xp[i], cool_func); + cooling_first_init_part(phys_const, us, cosmo, cool_func, &p[i], &xp[i]); /* Check part->gpart->part linkeage. */ if (p[i].gpart) p[i].gpart->id_or_neg_offset = -i; diff --git a/src/space.h b/src/space.h index 571b595eda8758c510a014bc784ae4bde85abaea..db9c9edad3d7435ea619e72814e2c97094344090 100644 --- a/src/space.h +++ b/src/space.h @@ -222,7 +222,7 @@ 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, + const struct chemistry_global_data *chemistry, const struct cooling_function_data *cool_func); void space_first_init_gparts(struct space *s, const struct gravity_props *grav_props);