/*******************************************************************************
* 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_CHIMES_IO_H
#define SWIFT_COOLING_CHIMES_IO_H
/* Config parameters. */
#include
/* Local includes */
#include "cooling.h"
#include "engine.h"
#include "io_properties.h"
#ifdef HAVE_HDF5
/**
* @brief Writes the current model of SPH to the file
* @param h_grp The HDF5 group in which to write
* @param h_grp_columns The HDF5 group containing named columns
* @param cooling the parameters of the cooling function.
*/
__attribute__((always_inline)) INLINE static void cooling_write_flavour(
hid_t h_grp, hid_t h_grp_columns,
const struct cooling_function_data* cooling) {
io_write_attribute_s(h_grp, "Cooling Model", "CHIMES");
/* Add the species names to the named columns */
hsize_t dims[1] = {CHIMES_NETWORK_SIZE};
hid_t type = H5Tcopy(H5T_C_S1);
H5Tset_size(type, CHIMES_NAME_STR_LENGTH);
hid_t space = H5Screate_simple(1, dims, NULL);
hid_t dset = H5Dcreate(h_grp_columns, "SpeciesFractions", type, space,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->chimes_species_names_reduced[0]);
H5Dclose(dset);
/* Let's construct the names from the raw full CHIMES list */
char chimes_species_names_snipshots[4][CHIMES_NAME_STR_LENGTH] = {0};
strcpy(chimes_species_names_snipshots[0], chimes_species_names[1]); /* HI */
strcpy(chimes_species_names_snipshots[1], chimes_species_names[4]); /* HeI */
strcpy(chimes_species_names_snipshots[2], chimes_species_names[5]); /* HeII */
strcpy(chimes_species_names_snipshots[3],
chimes_species_names[137]); /* H_2 */
/* Same for the reduced arrays in snipshots */
dims[0] = 4;
space = H5Screate_simple(1, dims, NULL);
dset = H5Dcreate(h_grp_columns, "ReducedSpeciesFractions", type, space,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dset, type, H5S_ALL, H5S_ALL, H5P_DEFAULT,
chimes_species_names_snipshots[0]);
H5Dclose(dset);
H5Tclose(type);
H5Sclose(space);
}
#endif
INLINE static void convert_part_T(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = cooling_get_snapshot_temperature(
e->physical_constants, e->hydro_properties, e->entropy_floor,
e->internal_units, e->cosmology, e->cooling_func, p, xp);
}
INLINE static void convert_part_sub_T(const struct engine* e,
const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = cooling_get_subgrid_temperature(
e->internal_units, e->physical_constants, e->cosmology,
e->hydro_properties, e->entropy_floor, e->cooling_func, p, xp);
}
INLINE static void convert_part_sub_rho(const struct engine* e,
const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = cooling_get_subgrid_density(
e->internal_units, e->physical_constants, e->cosmology,
e->hydro_properties, e->entropy_floor, e->cooling_func, p, xp);
}
INLINE static void convert_part_chimes_abundances(const struct engine* e,
const struct part* p,
const struct xpart* xp,
float* ret) {
/* Create dummy part and xpart structures that we can use to
* temporarily update the CHIMES abundance array for the snapshot. */
struct part dummy_p = *p;
struct xpart dummy_xp = *xp;
/* Update abundance array for particles on the EOS. */
cooling_set_subgrid_properties(e->physical_constants, e->internal_units,
e->cosmology, e->hydro_properties,
e->entropy_floor, e->cooling_func, e->dustevo,
&dummy_p, &dummy_xp);
for (int i = 0; i < CHIMES_NETWORK_SIZE; i++)
ret[i] = (float)dummy_xp.cooling_data.chimes_abundances[i];
}
INLINE static void convert_part_reduced_chimes_abundances(
const struct engine* e, const struct part* p, const struct xpart* xp,
float* ret) {
/* Create dummy part and xpart structures that we can use to
* temporarily update the CHIMES abundance array for the snapshot. */
struct part dummy_p = *p;
struct xpart dummy_xp = *xp;
const struct cooling_function_data* cooling = e->cooling_func;
/* Update abundance array for particles on the EOS. */
cooling_set_subgrid_properties(e->physical_constants, e->internal_units,
e->cosmology, e->hydro_properties,
e->entropy_floor, cooling, e->dustevo,
&dummy_p, &dummy_xp);
const int HI_ind = cooling->ChimesGlobalVars.speciesIndices[sp_HI];
const int HeI_ind = cooling->ChimesGlobalVars.speciesIndices[sp_HeI];
const int HeII_ind = cooling->ChimesGlobalVars.speciesIndices[sp_HeII];
const int H2_ind = cooling->ChimesGlobalVars.speciesIndices[sp_H2];
ret[0] = (float)dummy_xp.cooling_data.chimes_abundances[HI_ind]; /* HI */
ret[1] = (float)dummy_xp.cooling_data.chimes_abundances[HeI_ind]; /* HeI */
ret[2] = (float)dummy_xp.cooling_data.chimes_abundances[HeII_ind]; /* HeII */
ret[3] = (float)dummy_xp.cooling_data.chimes_abundances[H2_ind]; /* H_2 */
}
INLINE static void convert_part_HI_mass(const struct engine* e,
const struct part* p,
const struct xpart* xp, float* ret) {
const float X_H =
chemistry_get_metal_mass_fraction_for_cooling(p)[chemistry_element_H];
const float HI_frac = cooling_get_particle_subgrid_HI_fraction(
e->internal_units, e->physical_constants, e->cosmology,
e->hydro_properties, e->entropy_floor, e->cooling_func, p, xp);
*ret = hydro_get_mass(p) * X_H * HI_frac;
}
INLINE static void convert_part_H2_mass(const struct engine* e,
const struct part* p,
const struct xpart* xp, float* ret) {
const float X_H =
chemistry_get_metal_mass_fraction_for_cooling(p)[chemistry_element_H];
const float H2_frac = cooling_get_particle_subgrid_H2_fraction(
e->internal_units, e->physical_constants, e->cosmology,
e->hydro_properties, e->entropy_floor, e->cooling_func, p, xp);
*ret = hydro_get_mass(p) * X_H * H2_frac * 2.f;
}
INLINE static void convert_part_e_density(const struct engine* e,
const struct part* p,
const struct xpart* xp, double* ret) {
*ret = cooling_get_electron_density(
e->physical_constants, e->hydro_properties, e->entropy_floor,
e->internal_units, e->cosmology, e->cooling_func, e->dustevo, p, xp);
}
INLINE static void convert_part_y_compton(const struct engine* e,
const struct part* p,
const struct xpart* xp, double* ret) {
*ret = cooling_get_ycompton(e->physical_constants, e->hydro_properties,
e->entropy_floor, e->internal_units, e->cosmology,
e->cooling_func, e->dustevo, p, xp);
}
/**
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
* @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 part* parts, const struct xpart* xparts,
struct io_props* list) {
list[0] = io_make_output_field_convert_part(
"SpeciesFractions", FLOAT, CHIMES_NETWORK_SIZE, UNIT_CONV_NO_UNITS, 0.f,
parts, xparts, convert_part_chimes_abundances,
"Species fractions array for all ions and molecules in the CHIMES "
"network. The fraction of species i is defined in terms "
"of its number density relative to hydrogen, i.e. n_i / n_H_tot.");
list[1] = io_make_output_field_convert_part(
"ReducedSpeciesFractions", FLOAT, 4, UNIT_CONV_NO_UNITS, 0.f, parts,
xparts, convert_part_reduced_chimes_abundances,
"Species fractions array for a reduced set of ions and molecules in the "
"CHIMES network. The fraction of species i is defined in terms of its "
"number density relative to hydrogen, i.e. n_i / n_H_tot.");
list[2] = io_make_output_field_convert_part(
"Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
convert_part_T, "Temperature of the particles");
list[3] = io_make_output_field_convert_part(
"SubgridTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts,
xparts, convert_part_sub_T,
"The subgrid temperatures if the particles are within deltaT of the "
"entropy floor the subgrid temperature is calculated assuming a "
"pressure equilibrium on the entropy floor, if the particles are "
"above deltaT of the entropy floor the subgrid temperature is "
"identical to the SPH temperature.");
list[4] = io_make_physical_output_field_convert_part(
"SubgridPhysicalDensities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, parts,
xparts, /*can convert to comoving=*/1, convert_part_sub_rho,
"The subgrid physical density if the particles are within deltaT of the "
"entropy floor the subgrid density is calculated assuming a pressure "
"equilibrium on the entropy floor, if the particles are above deltaT "
"of the entropy floor the subgrid density is identical to the "
"physical SPH density.");
list[5] = io_make_output_field_convert_part(
"AtomicHydrogenMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts, xparts,
convert_part_HI_mass,
"Atomic hydrogen masses containted in the particles. This quantity is "
"obtained from the cooling tables and, if the particle is on the entropy "
"floor, by extrapolating to the equilibrium curve assuming constant "
"pressure.");
list[6] = io_make_output_field_convert_part(
"MolecularHydrogenMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts, xparts,
convert_part_H2_mass,
"Molecular hydrogen masses containted in the particles. This quantity is "
"obtained from the cooling tables and, if the particle is on the entropy "
"floor, by extrapolating to the equilibrium curve assuming constant "
"pressure.");
list[7] = io_make_physical_output_field(
"NetCoolingRates", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS_PER_UNIT_TIME,
0.f, xparts, cooling_data.net_cooling_rate, /*can convert to comoving=*/0,
"Net cooling rate of the particle, i.e. energy per unit mass per unit "
"time. Positive values indicate net cooling, negative values indicate "
"net heating.");
list[8] = io_make_physical_output_field_convert_part(
"TotalElectronNumberDensities", DOUBLE, 1, UNIT_CONV_NUMBER_DENSITY, 0.f,
parts, xparts, /*can convert to comoving=*/0, convert_part_e_density,
"Total electron number densities in the physical frame computed as the "
"sum of the electron number densities from the elements evolved with "
"the non-equilibrium network CHIMES and the electron number densities "
"from the other elements based on the equilibrium cooling tables.");
list[9] = io_make_physical_output_field_convert_part(
"ComptonYParameters", DOUBLE, 1, UNIT_CONV_AREA, 0.f, parts, xparts,
/*can convert to comoving=*/0, convert_part_y_compton,
"Compton y parameters in the physical frame computed based on the "
"cooling tables. This is 0 for star-forming particles.");
return 10;
}
#endif /* SWIFT_COOLING_CHIMES_IO_H */