/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch) * * 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_CHEMISTRY_GEAR_DIFFUSION_H #define SWIFT_CHEMISTRY_GEAR_DIFFUSION_H /** * @file src/chemistry/GEAR_DIFFUSION/chemistry.h * Follows Shen et al. 2010 */ /* Some standard headers. */ #include #include #include /* Local includes. */ #include "chemistry_struct.h" #include "error.h" #include "hydro.h" #include "parser.h" #include "part.h" #include "physical_constants.h" #include "units.h" /** * @brief Copies the chemistry properties of the gas particle over to the * star particle. * * @param p the gas particles. * @param xp the additional properties of the gas particles. * @param sp the new created star particle with its properties. */ INLINE static void chemistry_copy_star_formation_properties( struct part* p, const struct xpart* xp, struct spart* sp) { float mass = hydro_get_mass(p); /* Store the chemistry struct in the star particle */ for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { sp->chemistry_data.metal_mass_fraction[i] = p->chemistry_data.smoothed_metal_mass_fraction[i]; /* Remove the metals taken by the star. */ p->chemistry_data.metal_mass[i] *= mass / (mass + sp->mass); } } /** * @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 with diffusion'."); } /** * @brief Read the solar abundances and scale with them the initial * metallicities. * * @param parameter_file The parsed parameter file. * @param data The properties to initialise. */ static INLINE void chemistry_scale_initial_metallicities( struct swift_params* parameter_file, struct chemistry_global_data* data) { #ifdef HAVE_HDF5 /* Get the yields table */ char filename[DESCRIPTION_BUFFER_SIZE]; parser_get_param_string(parameter_file, "GEARFeedback:yields_table", filename); /* Open file. */ hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); if (file_id < 0) error("unable to open file %s.\n", filename); /* Open group. */ hid_t group_id = H5Gopen(file_id, "Data", H5P_DEFAULT); if (group_id < 0) error("unable to open group Data.\n"); /* Read the data */ float* sol_ab = (float*)malloc(sizeof(float) * GEAR_CHEMISTRY_ELEMENT_COUNT); io_read_array_attribute(group_id, "SolarMassAbundances", FLOAT, sol_ab, GEAR_CHEMISTRY_ELEMENT_COUNT); /* Close group */ hid_t status = H5Gclose(group_id); if (status < 0) error("error closing group."); /* Close file */ status = H5Fclose(file_id); if (status < 0) error("error closing file."); /* Scale the initial metallicities */ char txt[DESCRIPTION_BUFFER_SIZE] = "Scaling initial metallicities by:"; for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { data->initial_metallicities[i] *= sol_ab[i]; char tmp[10]; sprintf(tmp, " %.2g", sol_ab[i]); strcat(txt, tmp); } if (engine_rank == 0) { message("%s", txt); } #else error("Cannot scale the solar abundances without HDF5"); #endif } /** * @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(struct swift_params* parameter_file, const struct unit_system* us, const struct phys_const* phys_const, struct chemistry_global_data* data) { /* read parameters */ const float initial_metallicity = parser_get_param_float( parameter_file, "GEARChemistry:initial_metallicity"); if (initial_metallicity < 0) { message("Setting the initial metallicity from the snapshot."); } else { message("Setting the initial metallicity from the parameter file."); } /* Set the initial metallicities */ for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { data->initial_metallicities[i] = initial_metallicity; } /* Read the diffusion coefficient */ data->C = parser_get_param_float(parameter_file, "GEARChemistry:diffusion_coefficient"); /* Check if need to scale the initial metallicity */ const int scale_metallicity = parser_get_opt_param_int( parameter_file, "GEARChemistry:scale_initial_metallicity", 0); /* Scale the metallicities if required */ if (scale_metallicity) { chemistry_scale_initial_metallicities(parameter_file, data); } } /** * @brief Prepares a particle for the smooth metal calculation. * * Zeroes all the relevant arrays in preparation for the sums taking place in * the various smooth metallicity tasks * * @param p The particle to act upon * @param cd #chemistry_global_data containing chemistry informations. */ __attribute__((always_inline)) INLINE static void chemistry_init_part( struct part* restrict p, const struct chemistry_global_data* cd) { struct chemistry_part_data* cpd = &p->chemistry_data; for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { /* Reset the smoothed metallicity */ cpd->smoothed_metal_mass_fraction[i] = 0.f; /* Reset the diffusion equation */ cpd->metal_mass_dt[i] = 0; } /* Reset the shear tensor */ for (int i = 0; i < 3; i++) { cpd->S[i][0] = 0; cpd->S[i][1] = 0; cpd->S[i][2] = 0; } /* Reset the diffusion. */ cpd->diff_coef = 0; } /** * @brief Finishes the smooth metal calculation. * * Multiplies the smoothed metallicity and number of neighbours by the * appropiate constants and add the self-contribution term. * * This function requires the #hydro_end_density to have been called. * * @param p The particle to act upon. * @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_global_data* cd, const struct cosmology* cosmo) { /* Some smoothing length multiples. */ const float h = p->h; const float h_inv = 1.0f / h; /* 1/h */ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */ const float rho = hydro_get_comoving_density(p); const float factor = h_inv_dim / rho; /* 1 / h^d * rho */ const float rho_inv = 1.0f / rho; /* 1 / rho */ struct chemistry_part_data* cpd = &p->chemistry_data; for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { /* Final operation on the density (add self-contribution). */ cpd->smoothed_metal_mass_fraction[i] += cpd->metal_mass[i] * kernel_root; /* Finish the calculation by inserting the missing h-factors */ cpd->smoothed_metal_mass_fraction[i] *= factor; } /* convert the shear factor into physical */ const float factor_shear = h_inv_dim_plus_one * rho_inv * cosmo->a2_inv; for (int k = 0; k < 3; k++) { cpd->S[k][0] *= factor_shear; cpd->S[k][1] *= factor_shear; cpd->S[k][2] *= factor_shear; } /* Compute the trace over 3 and add the hubble flow. */ float trace_3 = 0; for (int i = 0; i < 3; i++) { cpd->S[i][i] += cosmo->H; trace_3 += cpd->S[i][i]; } trace_3 /= 3.; float S_t[3][3]; for (int i = 0; i < 3; i++) { /* Make the tensor symmetric. */ float avg = 0.5 * (cpd->S[i][0] + cpd->S[0][i]); S_t[i][0] = avg; S_t[0][i] = avg; avg = 0.5 * (cpd->S[i][1] + cpd->S[1][i]); S_t[i][1] = avg; S_t[1][i] = avg; avg = 0.5 * (cpd->S[i][2] + cpd->S[2][i]); S_t[i][2] = avg; S_t[2][i] = avg; /* Remove the trace. */ S_t[i][i] -= trace_3; } /* Compute the norm. */ float norm = 0; for (int i = 0; i < 3; i++) { norm += S_t[i][0] * S_t[i][0]; norm += S_t[i][1] * S_t[i][1]; norm += S_t[i][2] * S_t[i][2]; } norm = sqrt(norm); /* Compute the diffusion coefficient in physical coordinates. * The norm is already in physical coordinates. * We do not include kernel_gamma on purpose. */ const float h_phys = cosmo->a * p->h; cpd->diff_coef = cd->C * norm * h_phys * h_phys; } /** * @brief Updates to the chemistry data after the hydro force loop. * * @param p The particle to act upon. * @param cosmo The current cosmological model. * @param with_cosmology Are we running with the cosmology? * @param time Current time of the simulation. * @param dt Time step (in physical units). */ __attribute__((always_inline)) INLINE static void chemistry_end_force( struct part* restrict p, const struct cosmology* cosmo, const int with_cosmology, const double time, const double dt) { if (dt == 0) { return; } struct chemistry_part_data* ch = &p->chemistry_data; const float h_inv = cosmo->a / p->h; const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */ /* Missing factors in iact. */ const float factor = h_inv_dim * h_inv; const double sum = 0; for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { ch->metal_mass[i] += ch->metal_mass_dt[i] * dt * factor; /* Make sure that the metallicity is 0 <= x <= 1 */ if (ch->metal_mass[i] < 0 || ch->metal_mass[i] > hydro_get_mass(p)) { error("Negative mass or mass fraction larger than 1."); } /* Make sure that we do not have more metals than the sum. */ if (i != GEAR_CHEMISTRY_ELEMENT_COUNT) { sum += ch->metal_mass[i]; } else if (sum > ch->metal_mass[i]) { error("Found more individual elements than the sum of all of them."); } } } /** * @brief Sets all particle fields to sensible values when the #part has 0 ngbs. * * @param p The particle to act upon * @param xp The extended particle data to act upon * @param cd #chemistry_global_data containing chemistry informations. * @param cosmo The current cosmological model. */ __attribute__((always_inline)) INLINE static void chemistry_part_has_no_neighbours(struct part* restrict p, struct xpart* restrict xp, const struct chemistry_global_data* cd, const struct cosmology* cosmo) { /* Set the smoothed fractions with the non smoothed fractions */ for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { p->chemistry_data.smoothed_metal_mass_fraction[i] = p->chemistry_data.metal_mass[i] / hydro_get_mass(p); } } /** * @brief Computes the chemistry-related time-step constraint. * * No constraints in the GEAR model (no diffusion) --> FLT_MAX * * @param phys_const The physical constants in internal units. * @param cosmo The current cosmological model. * @param us The internal system of units. * @param hydro_props The properties of the hydro scheme. * @param cd The global properties of the chemistry scheme. * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float chemistry_timestep( const struct phys_const* restrict phys_const, const struct cosmology* restrict cosmo, const struct unit_system* restrict us, const struct hydro_props* hydro_props, const struct chemistry_global_data* cd, const struct part* restrict p) { return FLT_MAX; } /** * @brief Sets the chemistry properties of the (x-)particles to a valid start * state. * * @param phys_const The #phys_const. * @param us The #unit_system. * @param cosmo The #cosmology. * @param data The global chemistry information. * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data. */ __attribute__((always_inline)) INLINE static void chemistry_first_init_part( 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) { for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { if (data->initial_metallicities[i] < 0) { /* Use the value from the IC. We are reading the metal mass fraction. */ p->chemistry_data.metal_mass[i] *= hydro_get_mass(p); } else { /* Use the value from the parameter file */ p->chemistry_data.metal_mass[i] = data->initial_metallicities[i] * hydro_get_mass(p); } } chemistry_init_part(p, data); } /** * @brief Sets the chemistry properties of the sparticles to a valid start * state. * * @param data The global chemistry information. * @param sp Pointer to the sparticle data. */ __attribute__((always_inline)) INLINE static void chemistry_first_init_spart( const struct chemistry_global_data* data, struct spart* restrict sp) { for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { /* Bug fix (26.07.2024): Check that the initial me metallicities are non negative. */ if (data->initial_metallicities[i] >= 0) { /* Use the value from the parameter file */ sp->chemistry_data.metal_mass_fraction[i] = data->initial_metallicities[i]; } /* else : Use the value from the IC. We are reading the metal mass fraction. So do not overwrite the metallicities */ } } /** * @brief Sets the chemistry properties of the sink particles to a valid start * state. * * @param data The global chemistry information. * @param sink Pointer to the sink particle data. */ __attribute__((always_inline)) INLINE static void chemistry_first_init_sink( const struct chemistry_global_data* data, struct sink* restrict sink) { for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) { /* Use the value from the parameter file */ if (data->initial_metallicities[i] >= 0) { sink->chemistry_data.metal_mass_fraction[i] = data->initial_metallicities[i]; } /* else : read the metallicities from the ICs. */ } } /** * @brief Initialise the chemistry properties of a black hole with * the chemistry properties of the gas it is born from. * * Nothing to do here. * * @param bp_data The black hole data to initialise. * @param p_data The gas data to use. * @param gas_mass The mass of the gas particle. */ __attribute__((always_inline)) INLINE static void chemistry_bpart_from_part( struct chemistry_bpart_data* bp_data, const struct chemistry_part_data* p_data, const double gas_mass) { error("Loic: to be implemented"); } /** * @brief Add the chemistry data of a gas particle to a black hole. * * Nothing to do here. * * @param bp_data The black hole data to add to. * @param p_data The gas data to use. * @param gas_mass The mass of the gas particle. */ __attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart( struct chemistry_bpart_data* bp_data, const struct chemistry_part_data* p_data, const double gas_mass) { error("Loic: to be implemented"); } /** * @brief Add the chemistry data of a black hole to another one. * * Nothing to do here. * * @param bp_data The black hole data to add to. * @param swallowed_data The black hole data to use. */ __attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart( struct chemistry_bpart_data* bp_data, const struct chemistry_bpart_data* swallowed_data) { error("Loic: to be implemented"); } /** * @brief Split the metal content of a particle into n pieces * * @param p The #part. * @param n The number of pieces to split into. */ __attribute__((always_inline)) INLINE static void chemistry_split_part( struct part* p, const double n) { error("Loic: to be implemented"); } /** * @brief Returns the total metallicity (metal mass fraction) of the * star particle to be used in feedback/enrichment related routines. * * @param sp Pointer to the particle data. */ __attribute__((always_inline)) INLINE static double chemistry_get_total_metal_mass_fraction_for_feedback( const struct spart* restrict sp) { return sp->chemistry_data .metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1]; } /** * @brief Returns the abundances (metal mass fraction) of the * star particle to be used in feedback/enrichment related routines. * * @param sp Pointer to the particle data. */ __attribute__((always_inline)) INLINE static double const* chemistry_get_metal_mass_fraction_for_feedback( const struct spart* restrict sp) { return sp->chemistry_data.metal_mass_fraction; } /** * @brief Returns the total metallicity (metal mass fraction) of the * gas particle to be used in cooling related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static double chemistry_get_total_metal_mass_fraction_for_cooling( const struct part* restrict p) { return p->chemistry_data .smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1]; } /** * @brief Returns the abundance array (metal mass fractions) of the * gas particle to be used in cooling related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static double const* chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) { return p->chemistry_data.smoothed_metal_mass_fraction; } /** * @brief Returns the total metallicity (metal mass fraction) of the * gas particle to be used in star formation related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static double chemistry_get_total_metal_mass_fraction_for_star_formation( const struct part* restrict p) { return p->chemistry_data .smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1]; } /** * @brief Returns the abundance array (metal mass fractions) of the * gas particle to be used in star formation related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static double const* chemistry_get_metal_mass_fraction_for_star_formation( const struct part* restrict p) { return p->chemistry_data.smoothed_metal_mass_fraction; } /** * @brief Returns the total metallicity (metal mass fraction) of the * black hole particle to be used in the stats related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float chemistry_get_bh_total_metal_mass_for_stats(const struct bpart* restrict bp) { error("Not implemented"); return 0.f; } /** * @brief Returns the total metallicity (metal mass fraction) of the * gas particle to be used in the stats related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float chemistry_get_total_metal_mass_for_stats(const struct part* restrict p) { return p->chemistry_data.metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT - 1]; } /** * @brief Returns the total metallicity (metal mass fraction) of the * star particle to be used in the stats related routines. * * @param p Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) { return sp->chemistry_data .metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT - 1] * sp->mass; } #endif /* SWIFT_CHEMISTRY_GEAR_DIFFUSION_H */