/*******************************************************************************
* 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/CHIMES/colibre_tables.c
* @brief Functions to read COLIBRE tables
*/
/* Config parameters. */
#include
/* Standard includes */
#include
#include
#include
#include
/* Local includes. */
#include "chemistry_struct.h"
#include "cooling_properties.h"
#include "dust.h"
#include "error.h"
#include "exp10.h"
/**
* @brief Reads in COLIBRE table of redshift values
*
* @param cooling #cooling_function_data structure
*/
void get_cooling_redshifts(struct colibre_cooling_tables *table) {
/* Read the list of table redshifts */
char redshift_filename[colibre_table_path_name_length + 16];
sprintf(redshift_filename, "%s/redshifts.dat", table->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 */
if (engine_rank == 0)
message("Found cooling tables at %d redshifts", N_Redshifts);
table->number_of_redshifts = N_Redshifts;
/* Allocate the list of redshifts */
if (swift_memalign("cooling", (void **)&table->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) {
table->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 (table->Redshifts[i + 1] < table->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 table Colibre cooling table structure.
*/
void read_cooling_header(const char *fname,
struct colibre_cooling_tables *table) {
#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", table->cooling_table_path);
/* allocate arrays of bins */
if (posix_memalign((void **)&table->Temp, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_temperature * sizeof(float)) != 0)
error("Failed to allocate temperature table\n");
if (posix_memalign((void **)&table->nH, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_density * sizeof(float)) != 0)
error("Failed to allocate density table\n");
if (posix_memalign((void **)&table->Metallicity, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_metallicity * sizeof(float)) != 0)
error("Failed to allocate metallicity table\n");
if (posix_memalign((void **)&table->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 **)&table->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 **)&table->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 **)&table->atomicmass, SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_elementtypes * sizeof(float)) != 0)
error("Failed to allocate atomic masses array\n");
if (posix_memalign((void **)&table->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 **)&table->Zsol, SWIFT_STRUCT_ALIGNMENT,
1 * sizeof(float)) != 0)
error("Failed to allocate solar metallicity array\n");
if (posix_memalign((void **)&table->Zsol_inv, SWIFT_STRUCT_ALIGNMENT,
1 * sizeof(float)) != 0)
error("Failed to allocate inverse solar metallicity array\n");
if (posix_memalign((void **)&table->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 **)&table->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,
table->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,
table->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,
table->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, "/TotalAbundances", H5P_DEFAULT);
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
table->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,
table->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,
table->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,
table->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);
table->Zsol_inv[0] = 1.f / table->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(table->Metallicity[i]) < tol) {
table->indxZsol = i;
}
}
#if defined(__ICC)
#pragma novector
#endif
for (int i = 0; i < colibre_cooling_N_elementtypes; i++) {
table->atomicmass_inv[i] = 1.f / table->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 = cooling_row_major_index_2d(
i, j, colibre_cooling_N_metallicity, colibre_cooling_N_elementtypes);
table->Abundances[indx1d] = exp10f(table->LogAbundances[indx1d]);
table->Abundances_inv[indx1d] = 1.f / table->Abundances[indx1d];
table->MassFractions[indx1d] = exp10f(table->LogMassFractions[indx1d]);
}
}
#else
error("Need HDF5 to read cooling tables");
#endif
}
/**
* @brief Allocate space for cooling tables
*
* @param table Colibre cooling table structure.
*/
void allocate_cooling_tables(struct colibre_cooling_tables *table) {
/* 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. */
/* Cooling reduced (temperature) */
if (swift_memalign(
"cooling_table.Tcooling_reduced", (void **)&table->Tcooling_reduced,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_temperature *
colibre_cooling_N_metallicity * colibre_cooling_N_density *
colibre_reduced_cooling_N_cooltypes * sizeof(float)) != 0)
error("Failed to allocate Tcooling_reduced array\n");
/* Heating reduced (temperature) */
if (swift_memalign(
"cooling_table.Theating_reduced", (void **)&table->Theating_reduced,
SWIFT_STRUCT_ALIGNMENT,
colibre_cooling_N_loaded_redshifts * colibre_cooling_N_temperature *
colibre_cooling_N_metallicity * colibre_cooling_N_density *
colibre_reduced_cooling_N_heattypes * sizeof(float)) != 0)
error("Failed to allocate Theating_reduced array\n");
/* Electron fraction (temperature) */
if (swift_memalign(
"cooling_table.Tefrac", (void **)&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");
/* Thermal equilibrium temperature */
if (swift_memalign(
"cooling_table.Teq", (void **)&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 **)&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");
/* Pressure at thermal equilibrium temperature */
if (swift_memalign(
"cooling_table.Peq", (void **)&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 for tables");
#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 colibre_cooling_tables *table,
const int low_z_index, const int high_z_index,
float *log_depletion_fractions,
const int do_dust_pairing) {
#ifdef HAVE_HDF5
hid_t dataset;
herr_t status;
/* Temporary tables */
float *cooling_T = NULL;
float *heating_T = NULL;
float *electron_fraction_T = NULL;
float *equilibrium_temperature = NULL;
float *mean_particle_mass_Teq = NULL;
/* Allocate temporary tables */
/* 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");
/* 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");
/* Cooling, all channels, only loaded redshifts. */
if (swift_memalign(
"cooling_table.Tcooling", (void **)&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");
/* Heating, all channels, only loaded redshifts. */
if (swift_memalign(
"cooling_table.Theating", (void **)&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");
/* 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");
/* 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");
/* 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", table->cooling_table_path,
table->cooling_table_filebase, table->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/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 = cooling_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 = cooling_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);
table->Tcooling[internal_index] = cooling_T[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 = cooling_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 = cooling_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);
table->Theating[internal_index] = heating_T[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 = cooling_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 = cooling_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);
table->Telectron_fraction[internal_index] =
electron_fraction_T[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 = cooling_row_major_index_2d(
j, k, colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = cooling_row_major_index_3d(
local_z_index, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
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 = cooling_row_major_index_2d(
j, k, colibre_cooling_N_metallicity, colibre_cooling_N_density);
/* Index in the internal table */
const int internal_index = cooling_row_major_index_3d(
local_z_index, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
table->meanpartmass_Teq[internal_index] =
mean_particle_mass_Teq[hdf5_index];
}
}
/* Compute the pressures at thermal eq. for subgrid properties (if entropy
* floor is used) */
const float log10_kB_cgs = table->log10_kB_cgs;
for (int j = 0; j < colibre_cooling_N_metallicity; j++) {
const int index_XH = cooling_row_major_index_2d(
j, 0, colibre_cooling_N_metallicity, colibre_cooling_N_elementtypes);
const float log10_XH = table->LogMassFractions[index_XH];
for (int k = 0; k < colibre_cooling_N_density; k++) {
const int index_Peq = cooling_row_major_index_3d(
local_z_index, j, k, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
table->logPeq[index_Peq] =
table->nH[k] + table->logTeq[index_Peq] - log10_XH -
log10(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);
}
if (do_dust_pairing) {
/** Modify COLIBRE heating and cooling rates to remove the effects
* of implicit depletion when running with dust evolution **/
depletion_correct_rates_split(
table->Theating, table->Tcooling, log_depletion_fractions,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_elementtypes, colibre_cooling_N_cooltypes,
colibre_cooling_N_heattypes);
}
/** Reduce the cooling and heating tables to include only
* the channels required by CHIMES. **/
int idx_red, idx_T, idx_met, idx_dens, idx_elem, idx_table, idx_reduced_table;
float result;
for (idx_red = 0; idx_red < colibre_cooling_N_loaded_redshifts; idx_red++) {
for (idx_T = 0; idx_T < colibre_cooling_N_temperature; idx_T++) {
for (idx_met = 0; idx_met < colibre_cooling_N_metallicity; idx_met++) {
for (idx_dens = 0; idx_dens < colibre_cooling_N_density; idx_dens++) {
/* Copy the metal elements to the reduced table. */
for (idx_elem = 0; idx_elem < colibre_reduced_cooling_N_elementtypes;
idx_elem++) {
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, idx_elem + element_C,
colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_cooltypes);
idx_reduced_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, idx_elem,
colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_reduced_cooling_N_cooltypes);
table->Tcooling_reduced[idx_reduced_table] =
table->Tcooling[idx_table];
}
/* Copy eeBrems to the reduced table. */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, cooltype_eeBrems,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
idx_reduced_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, reduced_cooltype_eeBrems,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_reduced_cooling_N_cooltypes);
table->Tcooling_reduced[idx_reduced_table] =
table->Tcooling[idx_table];
/* Combine remaining channels if they are
* not already included in CHIMES. */
/* element_OA */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, element_OA,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
result = exp10f(table->Tcooling[idx_table]);
/* molecules */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, cooltype_molecules,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
result += exp10f(table->Tcooling[idx_table]);
/* HD */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, cooltype_HD,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
result += exp10f(table->Tcooling[idx_table]);
/* NetFFM */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, cooltype_NetFFM,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_cooltypes);
result += exp10f(table->Tcooling[idx_table]);
idx_reduced_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, reduced_cooltype_other,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_reduced_cooling_N_cooltypes);
table->Tcooling_reduced[idx_reduced_table] =
log10f(max(result, FLT_MIN));
}
}
}
}
for (idx_red = 0; idx_red < colibre_cooling_N_loaded_redshifts; idx_red++) {
for (idx_T = 0; idx_T < colibre_cooling_N_temperature; idx_T++) {
for (idx_met = 0; idx_met < colibre_cooling_N_metallicity; idx_met++) {
for (idx_dens = 0; idx_dens < colibre_cooling_N_density; idx_dens++) {
/* Copy the metal elements to the reduced table. */
for (idx_elem = 0; idx_elem < colibre_reduced_cooling_N_elementtypes;
idx_elem++) {
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, idx_elem + element_C,
colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_heattypes);
idx_reduced_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, idx_elem,
colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_reduced_cooling_N_heattypes);
table->Theating_reduced[idx_reduced_table] =
table->Theating[idx_table];
}
/* Combine remaining channels if they are
* not already included in CHIMES. */
/* element_OA */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, element_OA,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
result = exp10f(table->Theating[idx_table]);
/* COdiss */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, heattype_COdiss,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
result += exp10f(table->Theating[idx_table]);
/* UTA */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, heattype_UTA,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
result += exp10f(table->Theating[idx_table]);
/* line */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, heattype_line,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
result += exp10f(table->Theating[idx_table]);
/* Hlin */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, heattype_Hlin,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
result += exp10f(table->Theating[idx_table]);
/* ChaT */
idx_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, heattype_ChaT,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_heattypes);
result += exp10f(table->Theating[idx_table]);
idx_reduced_table = cooling_row_major_index_5d(
idx_red, idx_T, idx_met, idx_dens, reduced_heattype_other,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_reduced_cooling_N_heattypes);
table->Theating_reduced[idx_reduced_table] =
log10f(max(result, FLT_MIN));
}
}
}
}
swift_free("cooling_table.Tcooling", table->Tcooling);
swift_free("cooling_table.Theating", table->Theating);
swift_free("cooling-temp", cooling_T);
swift_free("cooling-temp", heating_T);
swift_free("cooling-temp", electron_fraction_T);
swift_free("cooling-temp", equilibrium_temperature);
swift_free("cooling-temp", mean_particle_mass_Teq);
#ifdef SWIFT_DEBUG_CHECKS
message("Done reading in general cooling table");
#endif
#else
error("Need HDF5 to read cooling tables");
#endif
}
/**
* @brief Computes net metal heating rate from Colibre tables.
*
* Computes the net heating rate (heating - cooling) for a given element
* abundance ratio, temperature, redshift, and density. The unit of the net
* cooling rate is Lambda / nH**2 [erg cm^3 s-1] and all input values are in
* cgs.
*
* @param myGasVars The #gasVariables struct.
* @param myGlobalVars The #globalVariables struct.
*/
ChimesFloat colibre_metal_cooling_rate_temperature(
struct gasVariables *myGasVars,
const struct globalVariables *myGlobalVars) {
struct global_hybrid_data_struct *myGlobalData;
myGlobalData = (struct global_hybrid_data_struct *)myGlobalVars->hybrid_data;
struct colibre_cooling_tables *table = myGlobalData->table;
float log_T_cgs = log10f(myGasVars->temperature);
float n_H_cgs = myGasVars->nH_tot;
float Z_absolute = myGasVars->metallicity * myGlobalData->Zsol;
/* Prevent electron fraction from
* becoming negative. */
float noneq_electron_fraction = max(myGasVars->abundances[sp_elec], 0.0f);
const float *abundance_ratio;
struct gas_hybrid_data_struct *myGasData;
myGasData = (struct gas_hybrid_data_struct *)myGasVars->hybrid_data;
abundance_ratio = myGasData->abundance_ratio;
/* Set weights for cooling rates */
float weights_cooling[colibre_reduced_cooling_N_cooltypes];
for (int i = 0; i < colibre_reduced_cooling_N_cooltypes; i++) {
if (i < colibre_reduced_cooling_N_elementtypes) {
/* Only include metals that aren't
* already included in CHIMES.
* Note that the abundance_ratio array
* includes H and He, so add element_C
* to the index. */
if (myGlobalVars->element_included[i] == 0)
weights_cooling[i] = abundance_ratio[i + element_C];
else
weights_cooling[i] = 0.f;
} else {
/* use same abundances as in the tables */
weights_cooling[i] = 1.f;
}
}
/* Set weights for heating rates */
float weights_heating[colibre_reduced_cooling_N_heattypes];
for (int i = 0; i < colibre_reduced_cooling_N_heattypes; i++) {
if (i < colibre_reduced_cooling_N_elementtypes) {
/* Only include metals that aren't
* already included in CHIMES.
* Note that the abundance_ratio array
* includes H and He, so add element_C
* to the index. */
if (myGlobalVars->element_included[i] == 0)
weights_heating[i] = abundance_ratio[i + element_C];
else
weights_heating[i] = 0.f;
} else {
weights_heating[i] = 1.f; /* use same abundances as in the tables */
}
}
// Set weights for electron densities.
float weights_electron[colibre_cooling_N_electrontypes];
for (int i = 0; i < colibre_cooling_N_electrontypes; i++) {
if (i < colibre_cooling_N_elementtypes - 1)
weights_electron[i] = abundance_ratio[i];
else
weights_electron[i] = 1.f;
}
/* Get indices of T, nH, and metallicity*/
int T_index, n_H_index, met_index, red_index;
float d_T, d_n_H, d_met;
float logZZsol = log10f((Z_absolute / table->Zsol[0]) + FLT_MIN);
cooling_get_index_1d(table->Temp, colibre_cooling_N_temperature, log_T_cgs,
&T_index, &d_T);
cooling_get_index_1d(table->Metallicity, colibre_cooling_N_metallicity,
logZZsol, &met_index, &d_met);
cooling_get_index_1d(table->nH, colibre_cooling_N_density, log10(n_H_cgs),
&n_H_index, &d_n_H);
if (table->dz > 0.5)
red_index = 1;
else
red_index = 0;
// n_e / n_H from Colibre
// Electron fraction from H + He + metals, with table metal ratios
float colibre_electron_fraction_table = interpolation3d_fix_x_plus_summation(
table->Telectron_fraction, weights_electron,
colibre_cooling_N_electrontypes - 1, colibre_cooling_N_electrontypes - 1,
red_index, T_index, met_index, n_H_index, d_T, d_met, d_n_H,
colibre_cooling_N_loaded_redshifts, colibre_cooling_N_temperature,
colibre_cooling_N_metallicity, colibre_cooling_N_density,
colibre_cooling_N_electrontypes);
float colibre_electron_fraction_metal_actual;
if (met_index == 0) {
/* Primordial abundances; metal
* contribution is zero. */
colibre_electron_fraction_metal_actual = 0.f;
} else {
/* From metals, with actual metal ratios.
* We also need to exclude any metals
* that are already included in CHIMES. */
for (int i = element_C; i < element_OA; i++) {
if (myGlobalVars->element_included[i - 2] == 1) weights_electron[i] = 0.f;
}
colibre_electron_fraction_metal_actual =
interpolation3d_fix_x_plus_summation(
table->Telectron_fraction, weights_electron, element_C,
colibre_cooling_N_electrontypes - 4, red_index, T_index, met_index,
n_H_index, d_T, d_met, d_n_H, colibre_cooling_N_loaded_redshifts,
colibre_cooling_N_temperature, colibre_cooling_N_metallicity,
colibre_cooling_N_density, colibre_cooling_N_electrontypes);
}
/* Adjust the weights based on the non-eq
* electron fraction from CHIMES, compared
* to the electron fraction that was used
* in the Colibre tables. */
float electron_fraction_ratio =
(noneq_electron_fraction + colibre_electron_fraction_metal_actual) /
(colibre_electron_fraction_table + FLT_MIN);
for (int i = 0; i < colibre_reduced_cooling_N_cooltypes; i++)
weights_cooling[i] *= electron_fraction_ratio;
/* eeBrems needs an extra factor of ne */
weights_cooling[reduced_cooltype_eeBrems] *= electron_fraction_ratio;
int cool_start, heat_start;
if (met_index == 0) {
/* Primordial abundances; metal
* contribution is zero. */
cool_start = reduced_cooltype_eeBrems;
heat_start = reduced_heattype_other;
} else {
cool_start = reduced_cooltype_C;
heat_start = reduced_heattype_C;
}
/* Lambda / n_H**2 */
const float cooling_rate = interpolation3d_fix_x_plus_summation(
table->Tcooling_reduced, weights_cooling, /* */
cool_start, colibre_reduced_cooling_N_cooltypes - 1, /* */
red_index, T_index, met_index, n_H_index, /* */
d_T, d_met, d_n_H, /* */
colibre_cooling_N_loaded_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
colibre_reduced_cooling_N_cooltypes); /* */
/* Gamma / n_H**2 */
const float heating_rate = interpolation3d_fix_x_plus_summation(
table->Theating_reduced, weights_heating, /* */
heat_start, colibre_reduced_cooling_N_heattypes - 1, /* */
red_index, T_index, met_index, n_H_index, /* */
d_T, d_met, d_n_H, /* */
colibre_cooling_N_loaded_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
colibre_reduced_cooling_N_heattypes); /* */
/* Return the net heating rate (Lambda_heat - Lambda_cool) */
return heating_rate - cooling_rate;
}
/**
* @brief Allocate gas_hybrid_data struct in gasVars.
*
* @param myGasVars The #gasVariables struct.
*/
void chimes_allocate_gas_hybrid_data(struct gasVariables *myGasVars) {
myGasVars->hybrid_data =
(void *)malloc(sizeof(struct gas_hybrid_data_struct));
}
/**
* @brief Allocate gas_hybrid_data struct in gasVars.
*
* @param myGasVars The #gasVariables struct.
*/
void chimes_free_gas_hybrid_data(struct gasVariables *myGasVars) {
free(myGasVars->hybrid_data);
}