/*******************************************************************************
* 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 .
*
******************************************************************************/
/**
* @file src/cooling/COLIBRE/cooling_tables.c
* @brief Functions to read COLIBRE tables
*/
/* Config parameters. */
#include
/* This file's header */
#include "cooling_tables.h"
/* Standard includes */
#include
#include
#include
#include
/* Local includes. */
#include "chemistry_struct.h"
#include "cooling_properties.h"
#include "error.h"
#include "exp10.h"
#include "interpolate.h"
/**
* @brief Reads in COLIBRE table of redshift values
*
* @param cooling #cooling_function_data structure
*/
void get_cooling_redshifts(struct cooling_function_data *cooling) {
/* Read the list of table redshifts */
char redshift_filename[colibre_table_path_name_length + 16];
sprintf(redshift_filename, "%s/redshifts.dat", cooling->cooling_table_path);
FILE *infile = fopen(redshift_filename, "r");
if (infile == NULL) {
error("Cannot open the list of cooling table redshifts (%s)",
redshift_filename);
}
int N_Redshifts = -1;
/* Read the file */
if (!feof(infile)) {
char buffer[50];
/* Read the number of redshifts (1st line in the file) */
if (fgets(buffer, 50, infile) != NULL)
N_Redshifts = atoi(buffer);
else
error("Impossible to read the number of redshifts");
/* Be verbose about it (only rank 0) */
if (engine_rank == 0)
message("Found cooling tables at %d redshifts", N_Redshifts);
cooling->number_of_redshifts = N_Redshifts;
/* Allocate the list of redshifts */
if (swift_memalign("cooling", (void **)&cooling->Redshifts,
SWIFT_STRUCT_ALIGNMENT,
N_Redshifts * sizeof(float)) != 0)
error("Failed to allocate redshift table");
/* Read all the redshift values */
int count = 0;
while (!feof(infile)) {
if (fgets(buffer, 50, infile) != NULL) {
cooling->Redshifts[count] = atof(buffer);
count++;
}
}
/* Verify that the file was self-consistent */
if (count != N_Redshifts) {
error(
"Redshift file (%s) does not contain the correct number of redshifts "
"(%d vs. %d)",
redshift_filename, count, N_Redshifts);
}
} else {
error("Redshift file (%s) is empty!", redshift_filename);
}
/* We are done with this file */
fclose(infile);
/* COLIBRE cooling assumes cooling->Redshifts table is in increasing order.
* Test this. */
for (int i = 0; i < N_Redshifts - 1; i++) {
if (cooling->Redshifts[i + 1] < cooling->Redshifts[i]) {
error("table should be in increasing order\n");
}
}
}
/**
* @brief Reads in COLIBRE cooling table header. Consists of tables
* of values for temperature, hydrogen number density, metallicity,
* abundance ratios, and elements used to index the cooling tables.
*
* @param fname Filepath for cooling table from which to read header
* @param cooling Cooling data structure
*/
void read_cooling_header(const char *fname,
struct cooling_function_data *cooling) {
#ifdef HAVE_HDF5
hid_t dataset;
herr_t status;
/* read sizes of array dimensions */
hid_t tempfile_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (tempfile_id < 0) error("unable to open file %s\n", fname);
/* allocate arrays of bins */
if (posix_memalign((void **)&cooling->Temp, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature * sizeof(float)) != 0)
error("Failed to allocate temperature table\n");
if (posix_memalign((void **)&cooling->nH, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate density table\n");
if (posix_memalign((void **)&cooling->Metallicity, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity * sizeof(float)) != 0)
error("Failed to allocate metallicity table\n");
if (posix_memalign((void **)&cooling->Therm, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_internalenergy * sizeof(float)) != 0)
error("Failed to allocate internal energy table\n");
if (posix_memalign((void **)&cooling->LogAbundances, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity *
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate abundance array\n");
if (posix_memalign((void **)&cooling->Abundances, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity *
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate abundance array\n");
if (posix_memalign((void **)&cooling->Abundances_inv, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity *
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate abundance array\n");
if (posix_memalign((void **)&cooling->atomicmass, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate atomic masses array\n");
if (posix_memalign((void **)&cooling->atomicmass_inv, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate inverse atomic masses array\n");
if (posix_memalign((void **)&cooling->Zsol, SWIFT_STRUCT_ALIGNMENT,
1 * sizeof(float)) != 0)
error("Failed to allocate solar metallicity array\n");
if (posix_memalign((void **)&cooling->Zsol_inv, SWIFT_STRUCT_ALIGNMENT,
1 * sizeof(float)) != 0)
error("Failed to allocate inverse solar metallicity array\n");
if (posix_memalign((void **)&cooling->LogMassFractions,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity *
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate log mass fraction array\n");
if (posix_memalign((void **)&cooling->MassFractions, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity *
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate mass fraction array\n");
/* read in bins and misc information */
dataset = H5Dopen(tempfile_id, "/TableBins/TemperatureBins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->Temp);
if (status < 0) error("error reading temperature bins\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/TableBins/DensityBins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->nH);
if (status < 0) error("error reading density bins\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/TableBins/MetallicityBins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->Metallicity);
if (status < 0) error("error reading metallicity bins\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/TableBins/InternalEnergyBins", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->Therm);
if (status < 0) error("error reading internal energy bins\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/TotalAbundances", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->LogAbundances);
if (status < 0) error("error reading total abundances\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/TotalMassFractions", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->LogMassFractions);
if (status < 0) error("error reading total mass fractions\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/ElementMasses", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->atomicmass);
if (status < 0) error("error reading element masses\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
dataset = H5Dopen(tempfile_id, "/SolarMetallicity", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling->Zsol);
if (status < 0) error("error reading solar metallicity \n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
/* Close the file */
H5Fclose(tempfile_id);
cooling->Zsol_inv[0] = 1.f / cooling->Zsol[0];
/* find the metallicity bin that refers to solar metallicity */
const float tol = 1.e-3;
for (int i = 0; i < colibre_cooling_N_metallicity; i++) {
if (fabsf(cooling->Metallicity[i]) < tol) {
cooling->indxZsol = i;
}
}
#if defined(__ICC)
#pragma novector
#endif
for (int i = 0; i < colibre_cooling_N_elementtypes; i++) {
cooling->atomicmass_inv[i] = 1.f / cooling->atomicmass[i];
}
/* set some additional useful abundance arrays */
for (int i = 0; i < colibre_cooling_N_metallicity; i++) {
#if defined(__ICC)
#pragma novector
#endif
for (int j = 0; j < colibre_cooling_N_elementtypes; j++) {
const int indx1d = row_major_index_2d(i, j, colibre_cooling_N_metallicity,
colibre_cooling_N_elementtypes);
cooling->Abundances[indx1d] = exp10f(cooling->LogAbundances[indx1d]);
cooling->Abundances_inv[indx1d] = 1.f / cooling->Abundances[indx1d];
cooling->MassFractions[indx1d] =
exp10f(cooling->LogMassFractions[indx1d]);
}
}
#else
error("Need HDF5 to read cooling tables");
#endif
}
/**
* @brief Allocate space for cooling tables.
*
* @param cooling #cooling_function_data structure
*/
void allocate_cooling_tables(struct cooling_function_data *restrict cooling) {
/* Allocate arrays to store cooling tables. Arrays contain two tables of
* cooling rates with one table being for the redshift above current redshift
* and one below. */
/* Mean particle mass (temperature) */
if (swift_memalign("cooling_table.Tmu", (void **)&cooling->table.Tmu,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_temperature *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate Tmu array\n");
/* Cooling (temperature) */
if (swift_memalign(
"cooling_table.Tcooling", (void **)&cooling->table.Tcooling,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_temperature *
colibre_cooling_N_metallicity * colibre_cooling_N_density *
colibre_cooling_N_cooltypes * sizeof(float)) != 0)
error("Failed to allocate Tcooling array\n");
/* Cooling (internal energy) */
if (swift_memalign("cooling_table.Ucooling",
(void **)&cooling->table.Ucooling, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_internalenergy *
colibre_cooling_N_metallicity *
colibre_cooling_N_density *
colibre_cooling_N_cooltypes * sizeof(float)) != 0)
error("Failed to allocate Ucooling array\n");
/* Heating (temperature) */
if (swift_memalign(
"cooling_table.Theating", (void **)&cooling->table.Theating,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_temperature *
colibre_cooling_N_metallicity * colibre_cooling_N_density *
colibre_cooling_N_heattypes * sizeof(float)) != 0)
error("Failed to allocate Theating array\n");
/* Heating (internal energy) */
if (swift_memalign("cooling_table.Uheating",
(void **)&cooling->table.Uheating, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_internalenergy *
colibre_cooling_N_metallicity *
colibre_cooling_N_density *
colibre_cooling_N_heattypes * sizeof(float)) != 0)
error("Failed to allocate Uheating array\n");
/* Electron fraction (temperature) */
if (swift_memalign(
"cooling_table.Tefrac", (void **)&cooling->table.Telectron_fraction,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_temperature *
colibre_cooling_N_metallicity * colibre_cooling_N_density *
colibre_cooling_N_electrontypes * sizeof(float)) != 0)
error("Failed to allocate Telectron_fraction array\n");
/* Electron fraction (internal energy) */
if (swift_memalign(
"cooling_table.Uefrac", (void **)&cooling->table.Uelectron_fraction,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_internalenergy * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_electrontypes *
sizeof(float)) != 0)
error("Failed to allocate Uelectron_fraction array\n");
/* Internal energy from temperature */
if (swift_memalign("cooling_table.UfromT", (void **)&cooling->table.U_from_T,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_temperature *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate U_from_T array\n");
/* Temperature from interal energy */
if (swift_memalign("cooling_table.TfromU", (void **)&cooling->table.T_from_U,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_internalenergy *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate T_from_U array\n");
/* Thermal equilibrium temperature */
if (swift_memalign("cooling_table.Teq", (void **)&cooling->table.logTeq,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate logTeq array\n");
/* Mean particle mass at thermal equilibrium temperature */
if (swift_memalign(
"cooling_table.mueq", (void **)&cooling->table.meanpartmass_Teq,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate mu array\n");
/* Hydrogen fractions at thermal equilibirum temperature */
if (swift_memalign(
"cooling_table.Hfracs", (void **)&cooling->table.logHfracs_Teq,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_hydrogentypes *
sizeof(float)) != 0)
error("Failed to allocate hydrogen fractions array\n");
/* All hydrogen fractions */
if (swift_memalign(
"cooling_table.Hfracs", (void **)&cooling->table.logHfracs_all,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_temperature *
colibre_cooling_N_metallicity * colibre_cooling_N_density *
colibre_cooling_N_hydrogentypes * sizeof(float)) != 0)
error("Failed to allocate big hydrogen fractions array\n");
/* Pressure at thermal equilibrium temperature */
if (swift_memalign("cooling_table.Peq", (void **)&cooling->table.logPeq,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate logPeq array\n");
#ifdef SWIFT_DEBUG_CHECKS
message("Allocated all cooling arrays");
#endif
}
/**
* @brief Get redshift dependent table of cooling rates.
* Reads in cooling tables (Ploeckinger & Schaye 2020) split in redshift
*
* @param cooling #cooling_function_data structure
* @param low_z_index Index of the lowest redshift table to load.
* @param high_z_index Index of the highest redshift table to load.
*/
void get_cooling_table(struct cooling_function_data *restrict cooling,
const int low_z_index, const int high_z_index,
float *log_depletion_fractions) {
#ifdef HAVE_HDF5
hid_t dataset;
herr_t status;
/* Temporary tables */
float *mean_particle_mass_T = NULL;
float *cooling_T = NULL;
float *cooling_U = NULL;
float *heating_T = NULL;
float *heating_U = NULL;
float *electron_fraction_T = NULL;
float *electron_fraction_U = NULL;
float *energy_from_temperature = NULL;
float *temperature_from_energy = NULL;
float *equilibrium_temperature = NULL;
float *mean_particle_mass_Teq = NULL;
float *hydrogen_fractions_Teq = NULL;
float *hydrogen_fractions_T = NULL;
/* Allocate temporary tables */
/* Mean particle mass (temperature) */
if (swift_memalign("cooling_temp", (void **)&mean_particle_mass_T,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate temporary Tmu array\n");
/* Cooling (temperature) */
if (swift_memalign(
"cooling_temp", (void **)&cooling_T, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_cooltypes *
sizeof(float)) != 0)
error("Failed to allocate temporary Tcooling array\n");
/* Cooling (internal energy) */
if (swift_memalign(
"cooling_temp", (void **)&cooling_U, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_internalenergy * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_cooltypes *
sizeof(float)) != 0)
error("Failed to allocate temporary Ucooling array\n");
/* Heating (temperature) */
if (swift_memalign(
"cooling_temp", (void **)&heating_T, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_heattypes *
sizeof(float)) != 0)
error("Failed to allocate temporary Theating array\n");
/* Heating (internal energy) */
if (swift_memalign(
"cooling_temp", (void **)&heating_U, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_internalenergy * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_heattypes *
sizeof(float)) != 0)
error("Failed to allocate temporary Uheating array\n");
/* Electron fraction (temperature) */
if (swift_memalign(
"cooling_temp", (void **)&electron_fraction_T, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_electrontypes *
sizeof(float)) != 0)
error("Failed to allocate temporary Telectron_fraction array\n");
/* Electron fraction (internal energy) */
if (swift_memalign(
"cooling_temp", (void **)&electron_fraction_U, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_internalenergy * colibre_cooling_N_metallicity *
colibre_cooling_N_density * colibre_cooling_N_electrontypes *
sizeof(float)) != 0)
error("Failed to allocate temporary Uelectron_fraction array\n");
/* Internal energy from temperature */
if (swift_memalign("cooling_temp", (void **)&energy_from_temperature,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate temporary U_from_T array\n");
/* Temperature from interal energy */
if (swift_memalign("cooling_temp", (void **)&temperature_from_energy,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_internalenergy *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate temporary T_from_U array\n");
/* Thermal equilibrium temperature */
if (swift_memalign("cooling_temp", (void **)&equilibrium_temperature,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity * colibre_cooling_N_density *
sizeof(float)) != 0)
error("Failed to allocate temporary logTeq array\n");
/* Mean particle mass at thermal equilibrium temperature */
if (swift_memalign("cooling_temp", (void **)&mean_particle_mass_Teq,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity * colibre_cooling_N_density *
sizeof(float)) != 0)
error("Failed to allocate temporary mu array\n");
/* Hydrogen fractions at thermal equilibirum temperature */
if (swift_memalign("cooling_temp", (void **)&hydrogen_fractions_Teq,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity * colibre_cooling_N_density *
3 * sizeof(float)) != 0)
error("Failed to allocate temporary hydrogen fractions array\n");
/* All hydrogen fractions */
if (swift_memalign("cooling_temp", (void **)&hydrogen_fractions_T,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature *
colibre_cooling_N_metallicity *
colibre_cooling_N_density * 3 * sizeof(float)) != 0)
error("Failed to allocate temporary big hydrogen fractions array\n");
/* Read in tables for redshift above and redshift below current value. */
for (int z_index = low_z_index; z_index <= high_z_index; z_index++) {
/* Index along redhsift dimension for the subset of tables we read */
const int local_z_index = z_index - low_z_index;
#ifdef SWIFT_DEBUG_CHECKS
if (local_z_index >= colibre_cooling_N_loaded_redshifts)
error("Reading invalid number of tables along z axis.");
#endif
/* Open table for this redshift index */
char fname[colibre_table_path_name_length + colibre_table_filebase_length +
12];
sprintf(fname, "%s/%s_z%4.2f.hdf5", cooling->cooling_table_path,
cooling->cooling_table_filebase, cooling->Redshifts[z_index]);
if (engine_rank == 0) message("Reading cooling table %s", fname);
hid_t file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) error("unable to open file %s", fname);
/* Read Tdep/MeanParticleMass into cooling->table.Tmu */
dataset = H5Dopen(file_id, "/Tdep/MeanParticleMass", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
mean_particle_mass_T);
if (status < 0) error("error reading temporary Tmu\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing temporary mean particle mass dataset");
for (int i = 0; i < colibre_cooling_N_temperature; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_3d(
i, j, k, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = row_major_index_4d(
local_z_index, i, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density);
cooling->table.Tmu[internal_index] = mean_particle_mass_T[hdf5_index];
}
}
}
/* Read Tdep/Cooling into cooling->table.Tcooling */
dataset = H5Dopen(file_id, "/Tdep/Cooling", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_T);
if (status < 0) error("error reading Tcooling\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_temperature; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_cooltypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_cooltypes);
cooling->table.Tcooling[internal_index] = cooling_T[hdf5_index];
}
}
}
}
/* Read Udep/Cooling into cooling->table.Ucooling */
dataset = H5Dopen(file_id, "/Udep/Cooling", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_U);
if (status < 0) error("error reading Ucooling\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_internalenergy; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_cooltypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_internalenergy,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_internalenergy, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_cooltypes);
cooling->table.Ucooling[internal_index] = cooling_U[hdf5_index];
}
}
}
}
/* Read Tdep/Heating into cooling->table.Theating */
dataset = H5Dopen(file_id, "/Tdep/Heating", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
heating_T);
if (status < 0) error("error reading Theating\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_temperature; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_heattypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_heattypes);
cooling->table.Theating[internal_index] = heating_T[hdf5_index];
}
}
}
}
/* Read Udep/Heating into cooling->table.Uheating */
dataset = H5Dopen(file_id, "/Udep/Heating", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
heating_U);
if (status < 0) error("error reading Uheating\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_internalenergy; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_heattypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_internalenergy,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_internalenergy, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_heattypes);
cooling->table.Uheating[internal_index] = heating_U[hdf5_index];
}
}
}
}
/* Read Tdep/ElectronFractions into cooling->table.Telectron_fraction */
/* Dataset is named /Tdep/ElectronFractions in the published version of the
* tables and for historical reasons /Tdep/ElectronFractionsVol in the
* version used in the COLIBRE repository. Content is identical but we deal
* here with both names */
if (H5Lexists(file_id, "/Tdep/ElectronFractionsVol", H5P_DEFAULT) > 0) {
dataset = H5Dopen(file_id, "/Tdep/ElectronFractionsVol", H5P_DEFAULT);
} else if (H5Lexists(file_id, "/Tdep/ElectronFractions", H5P_DEFAULT) > 0) {
dataset = H5Dopen(file_id, "/Tdep/ElectronFractions", H5P_DEFAULT);
} else {
error("Could not find the electron_fraction (temperature)!");
}
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
electron_fraction_T);
if (status < 0) error("error reading electron_fraction (temperature)\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_temperature; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_electrontypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_electrontypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_electrontypes);
cooling->table.Telectron_fraction[internal_index] =
electron_fraction_T[hdf5_index];
}
}
}
}
/* Read Udep/ElectronFractions into cooling->table.Uelectron_fraction */
/* Dataset is named /Udep/ElectronFractions in the published version of the
* tables and for historical reasons /Udep/ElectronFractionsVol in the
* version used in the COLIBRE repository. Content is identical but we deal
* here with both names */
if (H5Lexists(file_id, "/Udep/ElectronFractionsVol", H5P_DEFAULT) > 0) {
dataset = H5Dopen(file_id, "/Udep/ElectronFractionsVol", H5P_DEFAULT);
} else if (H5Lexists(file_id, "/Udep/ElectronFractions", H5P_DEFAULT) > 0) {
dataset = H5Dopen(file_id, "/Udep/ElectronFractions", H5P_DEFAULT);
} else {
error("Could not find the electron_fraction (internal energy)!");
}
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
electron_fraction_U);
if (status < 0)
error("error reading electron_fraction (internal energy)\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_internalenergy; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_electrontypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_internalenergy,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_electrontypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_internalenergy, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_electrontypes);
cooling->table.Uelectron_fraction[internal_index] =
electron_fraction_U[hdf5_index];
}
}
}
}
/* Read Tdep/U_from_T into cooling->table.U_from_T */
dataset = H5Dopen(file_id, "/Tdep/U_from_T", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
energy_from_temperature);
if (status < 0) error("error reading U_from_T array\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_temperature; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_3d(
i, j, k, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = row_major_index_4d(
local_z_index, i, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density);
cooling->table.U_from_T[internal_index] =
energy_from_temperature[hdf5_index];
}
}
}
/* Read Udep/T_from_U into cooling->table.T_from_U */
dataset = H5Dopen(file_id, "/Udep/T_from_U", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
temperature_from_energy);
if (status < 0) error("error reading T_from_U array\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing cooling dataset");
for (int i = 0; i < colibre_cooling_N_internalenergy; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_3d(
i, j, k, colibre_cooling_N_internalenergy,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = row_major_index_4d(
local_z_index, i, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_internalenergy, colibre_cooling_N_metallicity,
colibre_cooling_N_density);
cooling->table.T_from_U[internal_index] =
temperature_from_energy[hdf5_index];
}
}
}
/* Read ThermEq/Temperature into cooling->table.logTeq */
dataset = H5Dopen(file_id, "/ThermEq/Temperature", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
equilibrium_temperature);
if (status < 0) error("error reading Teq array\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing logTeq dataset");
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_2d(
j, k, colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = row_major_index_3d(
local_z_index, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
cooling->table.logTeq[internal_index] =
equilibrium_temperature[hdf5_index];
}
}
/* Read ThermEq/MeanParticleMass into cooling->table.meanpartmass_Teq */
dataset = H5Dopen(file_id, "/ThermEq/MeanParticleMass", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
mean_particle_mass_Teq);
if (status < 0) error("error reading mu array\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing mu dataset");
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_2d(
j, k, colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = row_major_index_3d(
local_z_index, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
cooling->table.meanpartmass_Teq[internal_index] =
mean_particle_mass_Teq[hdf5_index];
}
}
/* Read ThermEq/HydrogenFractionsVol into cooling->table.logHfracs_Teq */
dataset = H5Dopen(file_id, "/ThermEq/HydrogenFractionsVol", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
hydrogen_fractions_Teq);
if (status < 0) error("error reading hydrogen fractions array\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing hydrogen fractions dataset");
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_hydrogentypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_3d(
j, k, l, colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_hydrogentypes);
/* Index in the internal table */
const int internal_index = row_major_index_4d(
local_z_index, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_hydrogentypes);
cooling->table.logHfracs_Teq[internal_index] =
hydrogen_fractions_Teq[hdf5_index];
}
}
}
/* Read Tdep/HydrogenFractionsVol into cooling->table.logHfracs_all */
dataset = H5Dopen(file_id, "/Tdep/HydrogenFractionsVol", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
hydrogen_fractions_T);
if (status < 0) error("error reading big hydrogen fractions array\n");
status = H5Dclose(dataset);
if (status < 0) error("error closing big hydrogen fractions dataset");
for (int i = 0; i < colibre_cooling_N_temperature; i++) {
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
for (int k = 0; k < colibre_cooling_N_density; k++) {
for (int l = 0; l < colibre_cooling_N_hydrogentypes; l++) {
/* Index in the HDF5 table */
const int hdf5_index = row_major_index_4d(
i, j, k, l, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_hydrogentypes);
/* Index in the internal table */
const int internal_index = row_major_index_5d(
local_z_index, i, j, k, l, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_hydrogentypes);
cooling->table.logHfracs_all[internal_index] =
hydrogen_fractions_T[hdf5_index];
}
}
}
}
/* Compute the pressures at thermal eq. for subgrid properties (if entropy
* floor is used) */
const float log10_kB_cgs = cooling->log10_kB_cgs;
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
const int index_XH = row_major_index_2d(
j, 0, colibre_cooling_N_metallicity, colibre_cooling_N_elementtypes);
const float log10_XH = cooling->LogMassFractions[index_XH];
for (int k = 0; k < colibre_cooling_N_density; k++) {
const int index_Peq = row_major_index_3d(
local_z_index, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
cooling->table.logPeq[index_Peq] =
cooling->nH[k] + cooling->table.logTeq[index_Peq] - log10_XH -
log10(cooling->table.meanpartmass_Teq[index_Peq]) + log10_kB_cgs;
}
}
/* read depletion table if running with dust evolution modelling */
dust_read_depletion(
file_id, &log_depletion_fractions,
1 /* number of redshifts in hdf5 file */, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_elementtypes);
H5Fclose(file_id);
}
swift_free("cooling-temp", mean_particle_mass_T);
swift_free("cooling-temp", cooling_T);
swift_free("cooling-temp", cooling_U);
swift_free("cooling-temp", heating_T);
swift_free("cooling-temp", heating_U);
swift_free("cooling-temp", electron_fraction_T);
swift_free("cooling-temp", electron_fraction_U);
swift_free("cooling-temp", energy_from_temperature);
swift_free("cooling-temp", temperature_from_energy);
swift_free("cooling-temp", equilibrium_temperature);
swift_free("cooling-temp", mean_particle_mass_Teq);
swift_free("cooling-temp", hydrogen_fractions_Teq);
swift_free("cooling-temp", hydrogen_fractions_T);
#ifdef SWIFT_DEBUG_CHECKS
message("Done reading in general cooling table");
#endif
#else
error("Need HDF5 to read cooling tables");
#endif
}