/******************************************************************************* * 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 }