diff --git a/.clang-format-ignore b/.clang-format-ignore index 36fa238b6..181cf4813 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -30,9 +30,9 @@ src/clib/grackle_units.c src/clib/index_helper.c src/clib/index_helper.h src/clib/init_misc_species_cool_rates.cpp -src/clib/initialize_UVbackground_data.c +src/clib/initialize_UVbackground_data.cpp src/clib/initialize_chemistry_data.cpp -src/clib/initialize_cloudy_data.c +src/clib/initialize_cloudy_data.cpp src/clib/initialize_dust_yields.cpp src/clib/initialize_rates.cpp src/clib/internal_types.hpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 6e3b757c8..c6f52dc03 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -90,8 +90,6 @@ add_library(Grackle_Grackle dynamic_api.c grackle_units.c index_helper.c - initialize_cloudy_data.c initialize_cloudy_data.h - initialize_UVbackground_data.c initialize_UVbackground_data.h rate_functions.c set_default_chemistry_parameters.c solve_chemistry.c @@ -114,8 +112,10 @@ add_library(Grackle_Grackle dust_props.hpp init_misc_species_cool_rates.cpp init_misc_species_cool_rates.hpp initialize_chemistry_data.cpp + initialize_cloudy_data.cpp initialize_cloudy_data.hpp initialize_dust_yields.cpp initialize_dust_yields.hpp initialize_rates.cpp initialize_rates.hpp + initialize_UVbackground_data.cpp initialize_UVbackground_data.hpp internal_types.cpp internal_types.hpp lookup_cool_rates1d.hpp make_consistent.cpp make_consistent.hpp @@ -124,6 +124,7 @@ add_library(Grackle_Grackle scale_fields.cpp scale_fields.hpp solve_rate_cool_g-cpp.cpp solve_rate_cool_g-cpp.h step_rate_newton_raphson.hpp + support/h5io.cpp support/h5io.hpp time_deriv_0d.hpp utils-cpp.cpp utils-cpp.hpp utils-field.hpp diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 2fa8a89ee..a762a8537 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -50,6 +50,7 @@ OBJS_CONFIG_LIB = \ rate_functions.lo \ rate_utils.lo \ gaussj_g.lo \ + support/h5io.lo \ utils.lo \ utils-cpp.lo \ internal_types.lo diff --git a/src/clib/grackle_macros.h b/src/clib/grackle_macros.h index 6daae89f3..dec4527fc 100644 --- a/src/clib/grackle_macros.h +++ b/src/clib/grackle_macros.h @@ -82,20 +82,6 @@ #define GRFLOAT_C(DBL_LITERAL) ( DBL_LITERAL ) #endif -/* HDF5 definitions */ - -#define HDF5_FILE_I4 H5T_STD_I32BE -#define HDF5_FILE_I8 H5T_STD_I64BE -#define HDF5_FILE_R4 H5T_IEEE_F32BE -#define HDF5_FILE_R8 H5T_IEEE_F64BE -#define HDF5_FILE_B8 H5T_STD_B8BE - -#define HDF5_I4 H5T_NATIVE_INT -#define HDF5_I8 H5T_NATIVE_LLONG -#define HDF5_R4 H5T_NATIVE_FLOAT -#define HDF5_R8 H5T_NATIVE_DOUBLE -#define HDF5_R16 H5T_NATIVE_LDOUBLE - /* Precision-dependent definitions */ #ifdef GRACKLE_FLOAT_4 @@ -129,7 +115,6 @@ #define FLOAT_UNDEFINED -99999.0 #define INT_UNDEFINED -99999 -#define MAX_LINE_LENGTH 512 #ifndef tiny #define tiny 1.0e-20 diff --git a/src/clib/initialize_UVbackground_data.c b/src/clib/initialize_UVbackground_data.c deleted file mode 100644 index 8a5bf9250..000000000 --- a/src/clib/initialize_UVbackground_data.c +++ /dev/null @@ -1,364 +0,0 @@ -/*********************************************************************** -/ -/ Initialize UV background data -/ -/ -/ Copyright (c) 2013, Enzo/Grackle Development Team. -/ -/ Distributed under the terms of the Enzo Public Licence. -/ -/ The full license is in the file LICENSE, distributed with this -/ software. -************************************************************************/ - -#include -#include -#include "hdf5.h" -#include "grackle.h" -#include "grackle_macros.h" - -#include "initialize_UVbackground_data.h" - -// function prototypes -int read_dataset(hid_t file_id, char *dset_name, double *buffer); - -/** - * Initializes an empty #UVBtable struct with zeros and NULLs. - */ -void initialize_empty_UVBtable_struct(UVBtable *table) -{ - table->Nz = 0LL; - table->z = NULL; - table->k24 = NULL; - table->k25 = NULL; - table->k26 = NULL; - table->k27 = NULL; - table->k28 = NULL; - table->k29 = NULL; - table->k30 = NULL; - table->k31 = NULL; - table->piHI = NULL; - table->piHeI = NULL; - table->piHeII = NULL; - table->crsHI = NULL; - table->crsHeII = NULL; - table->crsHeI = NULL; -} - - -// Initialize UV Background data -int initialize_UVbackground_data(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates) -{ - long long Nz; - - // Return if no UV background selected or using fully tabulated cooling. - if (my_chemistry->UVbackground == 0 || - my_chemistry->primordial_chemistry == 0) - return SUCCESS; - - - if (grackle_verbose) - fprintf(stdout, "Initializing UV background.\n"); - - - // Read in UV background data from hdf5 file. - - hid_t file_id, dset_id, dspace_id; - herr_t status; - herr_t h5_error = -1; - - if (grackle_verbose) - fprintf(stdout, "Reading UV background data from %s.\n", - my_chemistry->grackle_data_file); - file_id = H5Fopen(my_chemistry->grackle_data_file, - H5F_ACC_RDONLY, H5P_DEFAULT); - - - // Read Info dataset - - dset_id = H5Dopen(file_id, "/UVBRates/Info"); - if (dset_id == h5_error) { - fprintf(stderr, "Can't open 'Info' dataset in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - int strlen = (int)(H5Dget_storage_size(dset_id)); - char info_string[strlen+1]; - - hid_t memtype = H5Tcopy(H5T_C_S1); - H5Tset_size(memtype, strlen+1); - - status = H5Dread(dset_id, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, info_string); - if (status == h5_error) { - fprintf(stderr, "Failed to read dataset 'Info'.\n"); - return FAIL; - } - - H5Tclose(memtype); - H5Dclose(dset_id); - - - - // Open redshift dataset and get number of elements - - dset_id = H5Dopen(file_id, "/UVBRates/z"); - if (dset_id == h5_error) { - fprintf(stderr, "Can't open redshift dataset ('z') in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - dspace_id = H5Dget_space(dset_id); - if (dspace_id == h5_error) { - fprintf(stderr, "Error opening dataspace for dataset 'z' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - Nz = H5Sget_simple_extent_npoints(dspace_id); - if(Nz <= 0) { - fprintf(stderr, "Redshift dataset ('z') has inappropriate size = %lld in %s.\n", - Nz, my_chemistry->grackle_data_file); - return FAIL; - } - - H5Sclose(dspace_id); - H5Dclose(dset_id); - - - - // Now allocate memory for UV background table. - my_rates->UVbackground_table.Nz = Nz; - - my_rates->UVbackground_table.z = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k24 = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k25 = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k26 = malloc(Nz * sizeof(double)); - - if (my_chemistry->primordial_chemistry > 1) { - my_rates->UVbackground_table.k27 = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k28 = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k29 = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k30 = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k31 = malloc(Nz * sizeof(double)); - } - - my_rates->UVbackground_table.piHI = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.piHeII = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.piHeI = malloc(Nz * sizeof(double)); - - if (my_chemistry->self_shielding_method > 0){ - my_rates->UVbackground_table.crsHI = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.crsHeII = malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.crsHeI = malloc(Nz * sizeof(double)); - } - - - // Now read everything. - - - // *** Redshift *** - if(! read_dataset(file_id, "/UVBRates/z", - my_rates->UVbackground_table.z) ) { - fprintf(stderr, "Error reading dataset 'z' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k24 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k24", - my_rates->UVbackground_table.k24) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k24' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k25 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k25", - my_rates->UVbackground_table.k25) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k25' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k26 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k26", - my_rates->UVbackground_table.k26) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k26' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - if (my_chemistry->primordial_chemistry > 1) { - - // *** k27 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k27", - my_rates->UVbackground_table.k27) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k27' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k28 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k28", - my_rates->UVbackground_table.k28) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k28' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k29 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k29", - my_rates->UVbackground_table.k29) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k29' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k30 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k30", - my_rates->UVbackground_table.k30) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k30' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** k31 *** - if(! read_dataset(file_id, "/UVBRates/Chemistry/k31", - my_rates->UVbackground_table.k31) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k31' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - } - - // *** piHI *** - if(! read_dataset(file_id, "/UVBRates/Photoheating/piHI", - my_rates->UVbackground_table.piHI) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHI' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** piHeII *** - if(! read_dataset(file_id, "/UVBRates/Photoheating/piHeII", - my_rates->UVbackground_table.piHeII) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHeII' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // *** piHeI *** - if(! read_dataset(file_id, "/UVBRates/Photoheating/piHeI", - my_rates->UVbackground_table.piHeI) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHeI' in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - - if (my_chemistry->self_shielding_method > 0) { - // *** crsHI *** - if(! read_dataset(file_id, "/UVBRates/CrossSections/hi_avg_crs", - my_rates->UVbackground_table.crsHI) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/hi_avg_crs' in %s.\n", - my_chemistry->grackle_data_file); - fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); - return FAIL; - } - - // *** crsHeII *** - if(! read_dataset(file_id, "/UVBRates/CrossSections/heii_avg_crs", - my_rates->UVbackground_table.crsHeII) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/heii_avg_crs' in %s.\n", - my_chemistry->grackle_data_file); - fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); - return FAIL; - } - - // *** crsHeI *** - if(! read_dataset(file_id, "/UVBRates/CrossSections/hei_avg_crs", - my_rates->UVbackground_table.crsHeI) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/hei_avg_crs' in %s.\n", - my_chemistry->grackle_data_file); - fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); - return FAIL; - } - } - - - H5Fclose(file_id); - - - // Get min/max of redshift vector - my_rates->UVbackground_table.zmin = my_rates->UVbackground_table.z[0]; - my_rates->UVbackground_table.zmax = my_rates->UVbackground_table.z[Nz-1]; - - // Print out some information about the dataset just read in. - if (grackle_verbose) { - fprintf(stdout, "UV background information:\n"); - fprintf(stdout, " %s\n",info_string); - fprintf(stdout, " z_min = %6.3f\n z_max = %6.3f\n", - my_rates->UVbackground_table.zmin, - my_rates->UVbackground_table.zmax); - } - - // Set redshift on/off flags from data if not set. - if (my_chemistry->UVbackground_redshift_on <= FLOAT_UNDEFINED) { - my_chemistry->UVbackground_redshift_on = - my_rates->UVbackground_table.zmax; - if (grackle_verbose) - fprintf(stdout, "Setting UVbackground_redshift_on to %f.\n", - my_chemistry->UVbackground_redshift_on); - } - if (my_chemistry->UVbackground_redshift_fullon <= FLOAT_UNDEFINED) { - my_chemistry->UVbackground_redshift_fullon = - my_rates->UVbackground_table.zmax; - if (grackle_verbose) - fprintf(stdout, "Setting UVbackground_redshift_fullon to %f.\n", - my_chemistry->UVbackground_redshift_fullon); - } - if (my_chemistry->UVbackground_redshift_drop <= FLOAT_UNDEFINED) { - my_chemistry->UVbackground_redshift_drop = - my_rates->UVbackground_table.zmin; - if (grackle_verbose) - fprintf(stdout, "Setting UVbackground_redshift_drop to %f.\n", - my_chemistry->UVbackground_redshift_drop); - } - if (my_chemistry->UVbackground_redshift_off <= FLOAT_UNDEFINED) { - my_chemistry->UVbackground_redshift_off = - my_rates->UVbackground_table.zmin; - if (grackle_verbose) - fprintf(stdout, "Setting UVbackground_redshift_off to %f.\n", - my_chemistry->UVbackground_redshift_off); - } - - return SUCCESS; -} - - - -int read_dataset(hid_t file_id, char *dset_name, double *buffer) { - hid_t dset_id; - herr_t status; - herr_t h5_error = -1; - - dset_id = H5Dopen(file_id, dset_name); - if (dset_id == h5_error) { - fprintf(stderr, "Failed to open dataset 'z'.\n"); - return FAIL; - } - - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); - if (status == h5_error) { - fprintf(stderr, "Failed to read dataset 'z'.\n"); - return FAIL; - } - - H5Dclose(dset_id); - - return SUCCESS; -} diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp new file mode 100644 index 000000000..fafcfb813 --- /dev/null +++ b/src/clib/initialize_UVbackground_data.cpp @@ -0,0 +1,354 @@ +/*********************************************************************** +/ +/ Initialize UV background data +/ +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ + +#include +#include +#include "hdf5.h" +#include "grackle.h" +#include "grackle_macros.h" // FLOAT_UNDEFINED +#include "support/h5io.hpp" + +#include "initialize_UVbackground_data.hpp" + +namespace { // stuff inside an anonymous namespace is only used in this file + +/// Initializes an empty #UVBtable struct with zeros and nullptrs. +void initialize_empty_UVBtable_struct(UVBtable *table) +{ + table->Nz = 0LL; + table->z = nullptr; + table->k24 = nullptr; + table->k25 = nullptr; + table->k26 = nullptr; + table->k27 = nullptr; + table->k28 = nullptr; + table->k29 = nullptr; + table->k30 = nullptr; + table->k31 = nullptr; + table->piHI = nullptr; + table->piHeI = nullptr; + table->piHeII = nullptr; + table->crsHI = nullptr; + table->crsHeII = nullptr; + table->crsHeI = nullptr; +} + +} // anonymous namespace + +// Initialize UV Background data +int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates) +{ + initialize_empty_UVBtable_struct(&my_rates->UVbackground_table); + + // Return if no UV background selected or using fully tabulated cooling. + if (my_chemistry->UVbackground == 0 || + my_chemistry->primordial_chemistry == 0) + return GR_SUCCESS; + + + if (grackle_verbose) + std::fprintf(stdout, "Initializing UV background.\n"); + + + // Read in UV background data from hdf5 file. + if (grackle_verbose) + std::fprintf(stdout, "Reading UV background data from %s.\n", + my_chemistry->grackle_data_file); + hid_t file_id = H5Fopen(my_chemistry->grackle_data_file, + H5F_ACC_RDONLY, H5P_DEFAULT); + + + // Read Info dataset + int buflen = h5io::read_str_dataset(file_id, "/UVBRates/Info", 0, nullptr); + if (buflen < 0) { + std::fprintf(stderr, "Error loading \"/UVBRates/Info\" dataset in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + std::vector info_string(buflen); + if (h5io::read_str_dataset(file_id, "/UVBRates/Info", buflen, + info_string.data()) < 0) { + std::fprintf(stderr, "Error loading \"/UVBRates/Info\" dataset in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // Open redshift dataset and get number of elements + + const h5io::ArrayShape common_shape + = h5io::read_dataset_shape(file_id, "/UVBRates/z"); + if (!h5io::ArrayShape_is_valid(common_shape)) { + return GR_FAIL; // error messages are already printed + } else if (common_shape.ndim != 1 || common_shape.shape[0] < 0) { + std::fprintf( + stderr, + "Redshift dataset (\"/UVBRates/z\") in %s has inappropriate shape\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + long long Nz = static_cast(common_shape.shape[0]); + + // Now allocate memory for UV background table. + my_rates->UVbackground_table.Nz = Nz; + + my_rates->UVbackground_table.z = new double[Nz]; + my_rates->UVbackground_table.k24 = new double[Nz]; + my_rates->UVbackground_table.k25 = new double[Nz]; + my_rates->UVbackground_table.k26 = new double[Nz]; + + if (my_chemistry->primordial_chemistry > 1) { + my_rates->UVbackground_table.k27 = new double[Nz]; + my_rates->UVbackground_table.k28 = new double[Nz]; + my_rates->UVbackground_table.k29 = new double[Nz]; + my_rates->UVbackground_table.k30 = new double[Nz]; + my_rates->UVbackground_table.k31 = new double[Nz]; + } + + my_rates->UVbackground_table.piHI = new double[Nz]; + my_rates->UVbackground_table.piHeII = new double[Nz]; + my_rates->UVbackground_table.piHeI = new double[Nz]; + + if (my_chemistry->self_shielding_method > 0){ + my_rates->UVbackground_table.crsHI = new double[Nz]; + my_rates->UVbackground_table.crsHeII = new double[Nz]; + my_rates->UVbackground_table.crsHeI = new double[Nz]; + } + + + // Now read everything. + + using ::grackle::impl::h5io::read_dataset; + + // *** Redshift *** + if(! read_dataset(file_id, "/UVBRates/z", + my_rates->UVbackground_table.z, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset 'z' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k24 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k24", + my_rates->UVbackground_table.k24, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k24' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k25 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k25", + my_rates->UVbackground_table.k25, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k25' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k26 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k26", + my_rates->UVbackground_table.k26, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k26' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + if (my_chemistry->primordial_chemistry > 1) { + + // *** k27 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k27", + my_rates->UVbackground_table.k27, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k27' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k28 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k28", + my_rates->UVbackground_table.k28, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k28' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k29 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k29", + my_rates->UVbackground_table.k29, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k29' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k30 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k30", + my_rates->UVbackground_table.k30, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k30' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** k31 *** + if(! read_dataset(file_id, "/UVBRates/Chemistry/k31", + my_rates->UVbackground_table.k31, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k31' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + } + + // *** piHI *** + if(! read_dataset(file_id, "/UVBRates/Photoheating/piHI", + my_rates->UVbackground_table.piHI, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHI' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** piHeII *** + if(! read_dataset(file_id, "/UVBRates/Photoheating/piHeII", + my_rates->UVbackground_table.piHeII, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHeII' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + // *** piHeI *** + if(! read_dataset(file_id, "/UVBRates/Photoheating/piHeI", + my_rates->UVbackground_table.piHeI, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHeI' in %s.\n", + my_chemistry->grackle_data_file); + return GR_FAIL; + } + + + if (my_chemistry->self_shielding_method > 0) { + // *** crsHI *** + if(! read_dataset(file_id, "/UVBRates/CrossSections/hi_avg_crs", + my_rates->UVbackground_table.crsHI, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/hi_avg_crs' in %s.\n", + my_chemistry->grackle_data_file); + std::fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); + return GR_FAIL; + } + + // *** crsHeII *** + if(! read_dataset(file_id, "/UVBRates/CrossSections/heii_avg_crs", + my_rates->UVbackground_table.crsHeII, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/heii_avg_crs' in %s.\n", + my_chemistry->grackle_data_file); + std::fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); + return GR_FAIL; + } + + // *** crsHeI *** + if(! read_dataset(file_id, "/UVBRates/CrossSections/hei_avg_crs", + my_rates->UVbackground_table.crsHeI, + /* expected_shape = */ &common_shape) ) { + std::fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/hei_avg_crs' in %s.\n", + my_chemistry->grackle_data_file); + std::fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); + return GR_FAIL; + } + } + + + H5Fclose(file_id); + + + // Get min/max of redshift vector + my_rates->UVbackground_table.zmin = my_rates->UVbackground_table.z[0]; + my_rates->UVbackground_table.zmax = my_rates->UVbackground_table.z[Nz-1]; + + // Print out some information about the dataset just read in. + if (grackle_verbose) { + std::fprintf(stdout, "UV background information:\n"); + std::fprintf(stdout, " %s\n",info_string.data()); + std::fprintf(stdout, " z_min = %6.3f\n z_max = %6.3f\n", + my_rates->UVbackground_table.zmin, + my_rates->UVbackground_table.zmax); + } + + // Set redshift on/off flags from data if not set. + if (my_chemistry->UVbackground_redshift_on <= FLOAT_UNDEFINED) { + my_chemistry->UVbackground_redshift_on = + my_rates->UVbackground_table.zmax; + if (grackle_verbose) + std::fprintf(stdout, "Setting UVbackground_redshift_on to %f.\n", + my_chemistry->UVbackground_redshift_on); + } + if (my_chemistry->UVbackground_redshift_fullon <= FLOAT_UNDEFINED) { + my_chemistry->UVbackground_redshift_fullon = + my_rates->UVbackground_table.zmax; + if (grackle_verbose) + std::fprintf(stdout, "Setting UVbackground_redshift_fullon to %f.\n", + my_chemistry->UVbackground_redshift_fullon); + } + if (my_chemistry->UVbackground_redshift_drop <= FLOAT_UNDEFINED) { + my_chemistry->UVbackground_redshift_drop = + my_rates->UVbackground_table.zmin; + if (grackle_verbose) + std::fprintf(stdout, "Setting UVbackground_redshift_drop to %f.\n", + my_chemistry->UVbackground_redshift_drop); + } + if (my_chemistry->UVbackground_redshift_off <= FLOAT_UNDEFINED) { + my_chemistry->UVbackground_redshift_off = + my_rates->UVbackground_table.zmin; + if (grackle_verbose) + std::fprintf(stdout, "Setting UVbackground_redshift_off to %f.\n", + my_chemistry->UVbackground_redshift_off); + } + + return GR_SUCCESS; +} + +void grackle::impl::free_UVBtable(UVBtable *table) +{ + // define a function that loosely approximates GRACKLE_FREE + auto cleanup_fn = [](double*& ptr) { + if (ptr != nullptr) { + delete[] ptr; + ptr = nullptr; + } + }; + + cleanup_fn(table->z); + cleanup_fn(table->k24); + cleanup_fn(table->k25); + cleanup_fn(table->k26); + cleanup_fn(table->k27); + cleanup_fn(table->k28); + cleanup_fn(table->k29); + cleanup_fn(table->k30); + cleanup_fn(table->k31); + cleanup_fn(table->piHI); + cleanup_fn(table->piHeII); + cleanup_fn(table->piHeI); + cleanup_fn(table->crsHI); + cleanup_fn(table->crsHeII); + cleanup_fn(table->crsHeI); +} diff --git a/src/clib/initialize_UVbackground_data.h b/src/clib/initialize_UVbackground_data.hpp similarity index 76% rename from src/clib/initialize_UVbackground_data.h rename to src/clib/initialize_UVbackground_data.hpp index 5c4580a7b..f24587e2e 100644 --- a/src/clib/initialize_UVbackground_data.h +++ b/src/clib/initialize_UVbackground_data.hpp @@ -15,19 +15,15 @@ #include "grackle.h" -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - -/// Initializes an empty UVBtable struct with zeros and NULLs. -void initialize_empty_UVBtable_struct(UVBtable* table); +namespace grackle::impl { /// Initialize UV Background data int initialize_UVbackground_data(chemistry_data* my_chemistry, chemistry_data_storage* my_rates); -#ifdef __cplusplus -} // extern "C" -#endif /* __cplusplus */ +/// deallocate memory associated with UVBtable struct +void free_UVBtable(UVBtable* table); + +} // namespace grackle::impl #endif /* INITIALIZE_UVBACKGROUND_DATA_H */ diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 294dbf7d6..10dd5118b 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -21,10 +21,10 @@ #include "auto_general.h" #include "interp_table_utils.h" // free_interp_grid_ #include "init_misc_species_cool_rates.hpp" // free_misc_species_cool_rates -#include "initialize_cloudy_data.h" +#include "initialize_cloudy_data.hpp" #include "initialize_dust_yields.hpp" // free_dust_yields #include "initialize_rates.hpp" -#include "initialize_UVbackground_data.h" +#include "initialize_UVbackground_data.hpp" #include "internal_types.hpp" // drop_CollisionalRxnRateCollection #include "opaque_storage.hpp" // gr_opaque_storage #include "phys_constants.h" @@ -429,25 +429,25 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry, /* Primordial tables. */ read_data = my_chemistry->primordial_chemistry == 0; - if (initialize_cloudy_data(my_chemistry, my_rates, - &my_rates->cloudy_primordial, - "Primordial", my_units, read_data) == GR_FAIL) { + if (grackle::impl::initialize_cloudy_data( + my_chemistry, my_rates, &my_rates->cloudy_primordial, + "Primordial", my_units, read_data) != GR_SUCCESS) { fprintf(stderr, "Error in initialize_cloudy_data.\n"); return GR_FAIL; } /* Metal tables. */ read_data = my_chemistry->metal_cooling == TRUE; - if (initialize_cloudy_data(my_chemistry, my_rates, - &my_rates->cloudy_metal, - "Metals", my_units, read_data) == GR_FAIL) { + if (grackle::impl::initialize_cloudy_data( + my_chemistry, my_rates, &my_rates->cloudy_metal, + "Metals", my_units, read_data) != GR_SUCCESS) { fprintf(stderr, "Error in initialize_cloudy_data.\n"); return GR_FAIL; } /* Initialize UV Background data. */ - initialize_empty_UVBtable_struct(&(my_rates->UVbackground_table)); - if (initialize_UVbackground_data(my_chemistry, my_rates) == GR_FAIL) { + if (grackle::impl::initialize_UVbackground_data(my_chemistry, my_rates) + != GR_SUCCESS) { fprintf(stderr, "Error in initialize_UVbackground_data.\n"); return GR_FAIL; } @@ -593,31 +593,11 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry, GRACKLE_FREE(my_rates->grain_growth_rate); } - free_cloudy_data(&my_rates->cloudy_primordial, my_chemistry, /* primordial */ 1); - free_cloudy_data(&my_rates->cloudy_metal, my_chemistry, /* primordial */ 0); - - GRACKLE_FREE(my_rates->UVbackground_table.z); - GRACKLE_FREE(my_rates->UVbackground_table.k24); - GRACKLE_FREE(my_rates->UVbackground_table.k25); - GRACKLE_FREE(my_rates->UVbackground_table.k26); - - if (my_chemistry->primordial_chemistry > 1) { - GRACKLE_FREE(my_rates->UVbackground_table.k27); - GRACKLE_FREE(my_rates->UVbackground_table.k28); - GRACKLE_FREE(my_rates->UVbackground_table.k29); - GRACKLE_FREE(my_rates->UVbackground_table.k30); - GRACKLE_FREE(my_rates->UVbackground_table.k31); - } - - GRACKLE_FREE(my_rates->UVbackground_table.piHI); - GRACKLE_FREE(my_rates->UVbackground_table.piHeII); - GRACKLE_FREE(my_rates->UVbackground_table.piHeI); - - if (my_chemistry->self_shielding_method > 0){ - GRACKLE_FREE(my_rates->UVbackground_table.crsHI); - GRACKLE_FREE(my_rates->UVbackground_table.crsHeII); - GRACKLE_FREE(my_rates->UVbackground_table.crsHeI); - } + grackle::impl::free_cloudy_data(&my_rates->cloudy_primordial, my_chemistry, + /* primordial = */ 1); + grackle::impl::free_cloudy_data(&my_rates->cloudy_metal, my_chemistry, + /* primordial */ 0); + grackle::impl::free_UVBtable(&my_rates->UVbackground_table); if (grackle::impl::free_misc_species_cool_rates(my_chemistry, my_rates) != GR_SUCCESS) { fprintf(stderr, "Error in free_metal_chemistry_rates.\n"); diff --git a/src/clib/initialize_cloudy_data.c b/src/clib/initialize_cloudy_data.c deleted file mode 100644 index 88ab57ada..000000000 --- a/src/clib/initialize_cloudy_data.c +++ /dev/null @@ -1,336 +0,0 @@ -/*********************************************************************** -/ -/ Initialize Cloudy cooling data -/ -/ -/ Copyright (c) 2013, Enzo/Grackle Development Team. -/ -/ Distributed under the terms of the Enzo Public Licence. -/ -/ The full license is in the file LICENSE, distributed with this -/ software. -************************************************************************/ - -#include -#include -#include -#include "hdf5.h" -#include "grackle.h" -#include "grackle_macros.h" -#include "initialize_cloudy_data.h" - -#define SMALL_LOG_VALUE -99.0 - - -/** - * Initializes an empty #cloudy_data struct with zeros and NULLs. - */ -void initialize_empty_cloudy_data_struct(cloudy_data *my_cloudy) -{ - my_cloudy->grid_rank = 0LL; - for (long long i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++){ - my_cloudy->grid_dimension[i] = 0LL; - my_cloudy->grid_parameters[i] = NULL; - } - my_cloudy->heating_data = NULL; - my_cloudy->cooling_data = NULL; - my_cloudy->mmw_data = NULL; - my_cloudy->data_size = 0LL; -} - -// initialize cloudy cooling data -int initialize_cloudy_data(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates, - cloudy_data *my_cloudy, const char *group_name, - code_units *my_units, int read_data) -{ - - long long q, w; - double *temp_data; - long long temp_int; - long long *temp_int_arr; - char parameter_name[MAX_LINE_LENGTH]; - - // Initialize things (to the null-state) even if cloudy cooling is not used. - initialize_empty_cloudy_data_struct(my_cloudy); - - if (read_data == 0) { return SUCCESS; } - - if (grackle_verbose) { - fprintf(stdout,"Initializing Cloudy cooling: %s.\n", group_name); - fprintf(stdout,"cloudy_table_file: %s.\n",my_chemistry->grackle_data_file); - } - - /* Get conversion units. */ - - double co_length_units, co_density_units; - if (my_units->comoving_coordinates == TRUE) { - co_length_units = my_units->length_units; - co_density_units = my_units->density_units; - } - else { - co_length_units = my_units->length_units * - my_units->a_value * my_units->a_units; - co_density_units = my_units->density_units / - POW(my_units->a_value * my_units->a_units, 3); - } - - double tbase1 = my_units->time_units; - double xbase1 = co_length_units / - (my_units->a_value * my_units->a_units); - double dbase1 = co_density_units * - POW(my_units->a_value * my_units->a_units, 3); - double mh = 1.67e-24; - double CoolUnit = (POW(my_units->a_units,5) * POW(xbase1,2) * POW(mh,2)) / - (POW(tbase1,3) * dbase1); - - // Read cooling data in from hdf5 file. - hid_t file_id, dset_id, attr_id; - herr_t status; - herr_t h5_error = -1; - - file_id = H5Fopen(my_chemistry->grackle_data_file, - H5F_ACC_RDONLY, H5P_DEFAULT); - - if (H5Aexists(file_id, "old_style")) { - my_rates->cloudy_data_new = 0; - if (grackle_verbose) - fprintf(stdout, "Loading old-style Cloudy tables.\n"); - } - - // Open cooling dataset and get grid dimensions. - - sprintf(parameter_name, "/CoolingRates/%s/Cooling", group_name); - dset_id = H5Dopen(file_id, parameter_name); - if (dset_id == h5_error) { - fprintf(stderr,"Can't open Cooling in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - // Grid rank. - attr_id = H5Aopen_name(dset_id, "Rank"); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to open Rank attribute in Cooling dataset.\n"); - return FAIL; - } - status = H5Aread(attr_id, HDF5_I8, &temp_int); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to read Rank attribute in Cooling dataset.\n"); - return FAIL; - } - my_cloudy->grid_rank = (long long) temp_int; - if (grackle_verbose) - fprintf(stdout,"Cloudy cooling grid rank: %lld.\n", my_cloudy->grid_rank); - status = H5Aclose(attr_id); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to close Rank attribute in Cooling dataset.\n"); - return FAIL; - } - - // Grid dimension. - temp_int_arr = malloc(my_cloudy->grid_rank * sizeof(long long)); - attr_id = H5Aopen_name(dset_id, "Dimension"); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to open Dimension attribute in Cooling dataset.\n"); - return FAIL; - } - status = H5Aread(attr_id, HDF5_I8,temp_int_arr); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to read Dimension attribute in Cooling dataset.\n"); - return FAIL; - } - if (grackle_verbose) - fprintf(stdout,"Cloudy cooling grid dimensions:"); - for (q = 0;q < my_cloudy->grid_rank;q++) { - my_cloudy->grid_dimension[q] = (long long) temp_int_arr[q]; - if (grackle_verbose) - fprintf(stdout," %lld", my_cloudy->grid_dimension[q]); - } - if (grackle_verbose) - fprintf(stdout,".\n"); - status = H5Aclose(attr_id); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to close Dimension attribute in Cooling dataset.\n"); - return FAIL; - } - free(temp_int_arr); - - // Grid parameters. - for (q = 0;q < my_cloudy->grid_rank;q++) { - - if (q < my_cloudy->grid_rank - 1) { - sprintf(parameter_name,"Parameter%lld",(q+1)); - } - else { - sprintf(parameter_name,"Temperature"); - } - - temp_data = malloc(my_cloudy->grid_dimension[q] * sizeof(double)); - - attr_id = H5Aopen_name(dset_id, parameter_name); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to open %s attribute in Cooling dataset.\n", - parameter_name); - return FAIL; - } - status = H5Aread(attr_id, HDF5_R8, temp_data); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to read %s attribute in Cooling dataset.\n", - parameter_name); - return FAIL; - } - - my_cloudy->grid_parameters[q] = malloc(my_cloudy->grid_dimension[q] * - sizeof(double)); - for (w = 0;w < my_cloudy->grid_dimension[q];w++) { - if (q < my_cloudy->grid_rank - 1) { - my_cloudy->grid_parameters[q][w] = (double) temp_data[w]; - } - else { - // convert temeperature to log - my_cloudy->grid_parameters[q][w] = (double) log10(temp_data[w]); - } - - } - if (grackle_verbose) - fprintf(stdout,"%s: %"GSYM" to %"GSYM" (%lld steps).\n",parameter_name, - my_cloudy->grid_parameters[q][0], - my_cloudy->grid_parameters[q][my_cloudy->grid_dimension[q]-1], - my_cloudy->grid_dimension[q]); - status = H5Aclose(attr_id); - if (attr_id == h5_error) { - fprintf(stderr,"Failed to close %s attribute in Cooling dataset.\n", - parameter_name); - return FAIL; - } - free(temp_data); - - } - - // Read Cooling data. - my_cloudy->data_size = 1; - for (q = 0;q < my_cloudy->grid_rank;q++) { - my_cloudy->data_size *= my_cloudy->grid_dimension[q]; - } - temp_data = malloc(my_cloudy->data_size * sizeof(double)); - - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data); - if (grackle_verbose) - fprintf(stdout,"Reading Cloudy Cooling dataset.\n"); - if (status == h5_error) { - fprintf(stderr,"Failed to read Cooling dataset.\n"); - return FAIL; - } - - my_cloudy->cooling_data = malloc(my_cloudy->data_size * sizeof(double)); - for (q = 0;q < my_cloudy->data_size;q++) { - my_cloudy->cooling_data[q] = temp_data[q] > 0 ? - (double) log10(temp_data[q]) : (double) SMALL_LOG_VALUE; - - // Convert to code units. - my_cloudy->cooling_data[q] -= log10(CoolUnit); - } - free(temp_data); - - status = H5Dclose(dset_id); - if (status == h5_error) { - fprintf(stderr,"Failed to close Cooling dataset.\n"); - return FAIL; - } - - // Read Heating data. - if (my_chemistry->UVbackground == 1) { - - temp_data = malloc(my_cloudy->data_size * sizeof(double)); - - sprintf(parameter_name, "/CoolingRates/%s/Heating", group_name); - dset_id = H5Dopen(file_id, parameter_name); - if (dset_id == h5_error) { - fprintf(stderr,"Can't open Heating in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data); - if (grackle_verbose) - fprintf(stdout,"Reading Cloudy Heating dataset.\n"); - if (status == h5_error) { - fprintf(stderr,"Failed to read Heating dataset.\n"); - return FAIL; - } - - my_cloudy->heating_data = malloc(my_cloudy->data_size * sizeof(double)); - for (q = 0;q < my_cloudy->data_size;q++) { - my_cloudy->heating_data[q] = temp_data[q] > 0 ? - (double) log10(temp_data[q]) : (double) SMALL_LOG_VALUE; - - // Convert to code units. - my_cloudy->heating_data[q] -= log10(CoolUnit); - } - free(temp_data); - - status = H5Dclose(dset_id); - if (status == h5_error) { - fprintf(stderr,"Failed to close Heating dataset.\n"); - return FAIL; - } - } - - // Read MMW data. - if (my_chemistry->primordial_chemistry == 0 && - strcmp(group_name, "Primordial") == 0) { - - my_cloudy->mmw_data = malloc(my_cloudy->data_size * sizeof(double)); - - sprintf(parameter_name, "/CoolingRates/%s/MMW", group_name); - dset_id = H5Dopen(file_id, parameter_name); - if (dset_id == h5_error) { - fprintf(stderr,"Can't open MMW in %s.\n", - my_chemistry->grackle_data_file); - return FAIL; - } - - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, - my_cloudy->mmw_data); - if (grackle_verbose) - fprintf(stdout,"Reading Cloudy MMW dataset.\n"); - if (status == h5_error) { - fprintf(stderr,"Failed to read MMW dataset.\n"); - return FAIL; - } - - status = H5Dclose(dset_id); - if (status == h5_error) { - fprintf(stderr,"Failed to close MMW dataset.\n"); - return FAIL; - } - } - - status = H5Fclose (file_id); - - if (my_cloudy->grid_rank > GRACKLE_CLOUDY_TABLE_MAX_DIMENSION) { - fprintf(stderr,"Error: rank of Cloudy cooling data must be less than or equal to %d.\n", - GRACKLE_CLOUDY_TABLE_MAX_DIMENSION); - return FAIL; - } - - return SUCCESS; -} - -int free_cloudy_data(cloudy_data *my_cloudy, chemistry_data *my_chemistry, int primordial) { - int i; - - for(i = 0; i < my_cloudy->grid_rank; i++) { - GRACKLE_FREE(my_cloudy->grid_parameters[i]); - } - - GRACKLE_FREE(my_cloudy->cooling_data); - if (my_chemistry->UVbackground == 1) { - GRACKLE_FREE(my_cloudy->heating_data); - } - if (my_chemistry->primordial_chemistry == 0 && primordial) { - GRACKLE_FREE(my_cloudy->mmw_data); - } - return GR_SUCCESS; -} diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp new file mode 100644 index 000000000..fa0aab74d --- /dev/null +++ b/src/clib/initialize_cloudy_data.cpp @@ -0,0 +1,290 @@ +/*********************************************************************** +/ +/ Initialize Cloudy cooling data +/ +/ +/ Copyright (c) 2013, Enzo/Grackle Development Team. +/ +/ Distributed under the terms of the Enzo Public Licence. +/ +/ The full license is in the file LICENSE, distributed with this +/ software. +************************************************************************/ + +#include // std::log10, std::pow +#include // std::fprintf, std::snprintf, stderr, stdio +#include // std::strcmp +#include "hdf5.h" +#include "grackle.h" +#include "grackle_macros.h" // TRUE +#include "initialize_cloudy_data.hpp" +#include "support/h5io.hpp" + +#define SMALL_LOG_VALUE (-99.0) +#define MAX_PARAMETER_NAME_LENGTH (512) + + +namespace { // stuff inside of an anonymous namespace is read-only + +/// Initializes an empty #cloudy_data struct with zeros and nullptrs. +void initialize_empty_cloudy_data_struct(cloudy_data *my_cloudy) +{ + my_cloudy->grid_rank = 0LL; + for (long long i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++){ + my_cloudy->grid_dimension[i] = 0LL; + my_cloudy->grid_parameters[i] = nullptr; + } + my_cloudy->heating_data = nullptr; + my_cloudy->cooling_data = nullptr; + my_cloudy->mmw_data = nullptr; + my_cloudy->data_size = 0LL; +} + +double* load_heatcool_data(hid_t file_id, const char* dset_name, + double CoolUnit, + grackle::impl::h5io::ArrayShape expected_shape) { + long long data_size = + grackle::impl::h5io::ArrayShape_elem_count(expected_shape); + double* tmp_data = new double[data_size]; + + if (grackle_verbose) { + std::fprintf(stdout,"Reading from \"%s\" dataset.\n", dset_name); + } + if (grackle::impl::h5io::read_dataset(file_id, dset_name, tmp_data, + &expected_shape) != GR_SUCCESS) { + delete[] tmp_data; + return nullptr; + } + + double* out = new double[data_size]; + for (long long q = 0LL; q < data_size; q++) { + out[q] = tmp_data[q] > 0 ? + (double) std::log10(tmp_data[q]) : (double) SMALL_LOG_VALUE; + + // Convert to code units. + out[q] -= std::log10(CoolUnit); + } + delete[] tmp_data; + + return out; +} + +} // anonymous namespace + + +// initialize cloudy cooling data +int grackle::impl::initialize_cloudy_data( + chemistry_data *my_chemistry, + chemistry_data_storage *my_rates, + cloudy_data *my_cloudy, const char *group_name, + code_units *my_units, int read_data) +{ + + char dset_name[MAX_PARAMETER_NAME_LENGTH]; + const std::size_t name_bufsize = + static_cast(MAX_PARAMETER_NAME_LENGTH); + + // Initialize things (to the null-state) even if cloudy cooling is not used. + initialize_empty_cloudy_data_struct(my_cloudy); + + if (read_data == 0) { return GR_SUCCESS; } + + if (grackle_verbose) { + std::fprintf(stdout,"Initializing Cloudy cooling: %s.\n", group_name); + std::fprintf(stdout,"cloudy_table_file: %s.\n",my_chemistry->grackle_data_file); + } + + /* Get conversion units. */ + + double co_length_units, co_density_units; + if (my_units->comoving_coordinates == TRUE) { + co_length_units = my_units->length_units; + co_density_units = my_units->density_units; + } + else { + co_length_units = my_units->length_units * + my_units->a_value * my_units->a_units; + co_density_units = my_units->density_units / + std::pow(my_units->a_value * my_units->a_units, 3.0); + } + + double tbase1 = my_units->time_units; + double xbase1 = co_length_units / + (my_units->a_value * my_units->a_units); + double dbase1 = co_density_units * + std::pow(my_units->a_value * my_units->a_units, 3.0); + double mh = 1.67e-24; + double CoolUnit = ( + std::pow(my_units->a_units,5.0) * std::pow(xbase1,2.0) * std::pow(mh,2.0) + ) / (std::pow(tbase1, 3.0) * dbase1); + + // Read cooling data in from hdf5 file. + hid_t file_id = H5Fopen(my_chemistry->grackle_data_file, + H5F_ACC_RDONLY, H5P_DEFAULT); + + if (H5Aexists(file_id, "old_style")) { + my_rates->cloudy_data_new = 0; + if (grackle_verbose) + std::fprintf(stdout, "Loading old-style Cloudy tables.\n"); + } + + // Determine the name of the dataset + std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Cooling", + group_name); + + // Parse the grid properties from the dataset's attributes + h5io::GridTableProps grid_props = h5io::parse_GridTableProps(file_id, + dset_name); + if (!h5io::GridTableProps_is_valid(grid_props)) { + h5io::drop_GridTableProps(&grid_props); + H5Fclose (file_id); + // error messages were already printed by h5io::parse_GridTableProps + return GR_FAIL; + } + + // use grid_props to initialize parts of my_cloudy + // TODO: it would be really nice if we explicitly validated that each + // grid_props.axes[i] held the expected values + my_cloudy->grid_rank = static_cast(grid_props.table_shape.ndim); + for (int i = 0; i < grid_props.table_shape.ndim; i++) { + my_cloudy->grid_dimension[i] + = static_cast(grid_props.table_shape.shape[i]); + + bool is_temperature = + std::strcmp("Temperature", grid_props.axes[i].name) == 0; + + double* buf = new double[my_cloudy->grid_dimension[i]]; + for (long long w = 0LL; w < my_cloudy->grid_dimension[i]; w++) { + if (is_temperature) { + buf[w] = std::log10(grid_props.axes[i].values[w]); + } else { + buf[w] = grid_props.axes[i].values[w]; + } + } + my_cloudy->grid_parameters[i] = buf; + } + + if (grackle_verbose) { + std::fprintf(stdout, "Cloudy cooling grid = {\n"); + std::fprintf(stdout, " rank: %lld,\n", my_cloudy->grid_rank); + for (long long i = 0LL; i < my_cloudy->grid_rank; i++) { + std::fprintf(stdout, + " axis %lld (\"%s\"): %g to %g (%lld steps),\n", + i, + grid_props.axes[i].name, + my_cloudy->grid_parameters[i][0], + my_cloudy->grid_parameters[i][my_cloudy->grid_dimension[i]-1], + my_cloudy->grid_dimension[i]); + } + fprintf(stdout, "}\n"); + } + h5io::ArrayShape expected_shape = grid_props.table_shape; + + // Read Cooling data. + my_cloudy->data_size = h5io::ArrayShape_elem_count(expected_shape); + my_cloudy->cooling_data = load_heatcool_data(file_id, dset_name, CoolUnit, + expected_shape); + + if (my_cloudy->cooling_data == nullptr) { + h5io::drop_GridTableProps(&grid_props); + H5Fclose(file_id); + return GR_FAIL; + } + + // Read Heating data. + if (my_chemistry->UVbackground == 1) { + std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Heating", + group_name); + + // Ideally we would uncomment the following block of logic that verifies + // that the Grid Table properties are identical to the cooling table + // -> uncommenting the logic would be useful for enforcing consistency + // in any new data files (all files currently shipped with Grackle + // already satisfy this requirement) + // -> the only downside to uncommenting the block is a little overhead + + //if (h5io::assert_has_consistent_GridTableProps(file_id, dset_name, + // grid_props) != GR_SUCCESS){ + // h5io::drop_GridTableProps(&grid_props); + // H5Fclose(file_id); + // return GR_FAIL; + //} + my_cloudy->heating_data = load_heatcool_data(file_id, dset_name, CoolUnit, + expected_shape); + if (my_cloudy->heating_data == nullptr) { + h5io::drop_GridTableProps(&grid_props); + H5Fclose(file_id); + return GR_FAIL; + } + } + + // Read MMW data. + if (my_chemistry->primordial_chemistry == 0 && + std::strcmp(group_name, "Primordial") == 0) { + const char* mmw_dset_name = "/CoolingRates/Primordial/MMW"; + + // Ideally we would uncomment the following block of logic that verifies + // that the Grid Table properties are identical to the cooling table + // -> uncommenting the logic would be useful for enforcing consistency + // in any new data files (all files currently shipped with Grackle + // already satisfy this requirement) + // -> the only downside to uncommenting the block is a little overhead + + //if (h5io::assert_has_consistent_GridTableProps(file_id, mmw_dset_name, + // grid_props) != GR_SUCCESS){ + // h5io::drop_GridTableProps(&grid_props); + // H5Fclose(file_id); + // return GR_FAIL; + //} + + my_cloudy->mmw_data = new double[my_cloudy->data_size]; + if (h5io::read_dataset(file_id, mmw_dset_name, my_cloudy->mmw_data, + &expected_shape) != GR_SUCCESS) { + h5io::drop_GridTableProps(&grid_props); + H5Fclose(file_id); + return GR_FAIL; + } + + } + + H5Fclose (file_id); + + if (my_cloudy->grid_rank > GRACKLE_CLOUDY_TABLE_MAX_DIMENSION) { + std::fprintf( + stderr, + "Error: rank of Cloudy cooling data must be less than or equal to %d.\n", + GRACKLE_CLOUDY_TABLE_MAX_DIMENSION); + return GR_FAIL; + } + + return GR_SUCCESS; +} + +int grackle::impl::free_cloudy_data(cloudy_data *my_cloudy, + chemistry_data *my_chemistry, + int primordial) { + + for(int i = 0; i < my_cloudy->grid_rank; i++) { + if (my_cloudy->grid_parameters[i] != nullptr) { + delete[] my_cloudy->grid_parameters[i]; + my_cloudy->grid_parameters[i] = nullptr; + } + } + + if (my_cloudy->cooling_data != nullptr) { + delete[] my_cloudy->cooling_data; + my_cloudy->cooling_data = nullptr; + } + + if (my_cloudy->heating_data != nullptr) { + delete[] my_cloudy->heating_data; + my_cloudy->heating_data = nullptr; + } + + if (my_cloudy->mmw_data != nullptr) { + delete[] my_cloudy->mmw_data; + my_cloudy->mmw_data = nullptr; + } + return GR_SUCCESS; +} + diff --git a/src/clib/initialize_cloudy_data.h b/src/clib/initialize_cloudy_data.hpp similarity index 79% rename from src/clib/initialize_cloudy_data.h rename to src/clib/initialize_cloudy_data.hpp index 4159b8ad0..dfc73c977 100644 --- a/src/clib/initialize_cloudy_data.h +++ b/src/clib/initialize_cloudy_data.hpp @@ -15,12 +15,7 @@ #include "grackle.h" -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - -/// Initializes an empty #cloudy_data struct with zeros and NULLs. -void initialize_empty_cloudy_data_struct(cloudy_data* my_cloudy); +namespace grackle::impl { // initialize cloudy cooling data int initialize_cloudy_data(chemistry_data* my_chemistry, @@ -31,8 +26,6 @@ int initialize_cloudy_data(chemistry_data* my_chemistry, int free_cloudy_data(cloudy_data* my_cloudy, chemistry_data* my_chemistry, int primordial); -#ifdef __cplusplus -} // extern "C" -#endif /* __cplusplus */ +} // namespace grackle::impl #endif /* INITIALIZE_CLOUDY_DATA_H */ diff --git a/src/clib/support/h5io.cpp b/src/clib/support/h5io.cpp new file mode 100644 index 000000000..f235031b6 --- /dev/null +++ b/src/clib/support/h5io.cpp @@ -0,0 +1,944 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Implements functions to help read from hdf5 files +/// +/// The error-handling in these functions leaves a lot to be desired... +/// +//===----------------------------------------------------------------------===// +#include +#include +#include + +#include "hdf5.h" +#include "grackle.h" +#include "grackle_macros.h" +#include "h5io.hpp" +#include "../status_reporting.h" + +// HDF5 definitions +// -> these were moved out of the grackle_macros.h header +// -> most of these are unnecessary! + +#define HDF5_FILE_I4 H5T_STD_I32BE +#define HDF5_FILE_I8 H5T_STD_I64BE +#define HDF5_FILE_R4 H5T_IEEE_F32BE +#define HDF5_FILE_R8 H5T_IEEE_F64BE +#define HDF5_FILE_B8 H5T_STD_B8BE + +#define HDF5_I4 H5T_NATIVE_INT +#define HDF5_I8 H5T_NATIVE_LLONG +#define HDF5_R4 H5T_NATIVE_FLOAT +#define HDF5_R8 H5T_NATIVE_DOUBLE + +namespace { // stuff inside an anonymous namespace is local to this file + +/// Determines whether the specified buffer contains a null-terminated byte +/// string of ascii characters. +/// +/// @returns true if buffer only contains ascii characters and contains a null +/// character at buffer[i] for `i < bufsz`. Otherwise it returns true. +/// +/// @par Purpose +/// The function is primarily intended for checking whether a UTF-8 encoded +/// string is also a valid ASCII string (i.e. it only includes UTF-8 code +/// points that are valid ASCII characters). For the uninitiated: +/// - all 128 standard ASCII characters are encoded in 7 bits. UTF-8 was +/// designed for backwards compatibility with these characters; the first +/// 128 UTF-8 code points are all encoded using a single byte and have an +/// **EXACT** one-to-one correspondence with the ascii characters. +/// - all other utf8 codepoints are represented by 2 or more bytes & +/// the leading byte has a value doesn't correspond to an ASCII character +bool is_ascii_string(const char* buffer, int bufsz) { + const char max_val = static_cast(127); + for (int i = 0; i < bufsz; i++) { + if (buffer[i] == '\0') { + return true; + } else if (buffer[i] < '\0' || buffer[i] > max_val) { + // I'm pretty sure we must verfiy that buffer[i] BOTH isn't smaller + // than '\0' AND doesn't exceed max_val since the standard doesn't + // specify whether char is signed or unsigned + return false; + } + } + + return false; // (buffer doesn't contain a null character) +} + +/// does the heavy lifting for read_str_attribute and read_str_dataset +int read_str_data_helper_(hid_t id, bool is_attr, int bufsz, char* buffer) { + if (bufsz < 0 || id == H5I_INVALID_HID) { + return -1; + } + + // I'm a little suspicious of this line of code... + // -> if we ever see Parameter names get clipped, my guess is that this line + // doesn't handle variable length strings properly + // -> with that said, I *think* this is correct + hsize_t storage_size = + (is_attr) ? H5Aget_storage_size(id) : H5Dget_storage_size(id); + if (storage_size > static_cast(INT_MAX - 1)) { + std::fprintf( + stderr, + "the string attribute requires a bigger buffer than expected\n"); + return -1; + } + int required_bufsz = static_cast(storage_size) + 1; + + if (bufsz == 0) { + return required_bufsz; + } else if (required_bufsz > bufsz) { + return -1; + } + + // we inspect a few properties about the datatype on disk + // - its unfortunate that we need to worry about utf8 encoding... This is + // relevant for the _high_density table (I suspect that it may be a default + // within h5py). We'll validate that the characters are compatible with + // ASCII before we return + bool is_variable, uses_utf8_encoding; + { + hid_t disk_typeid = (is_attr) ? H5Aget_type(id) : H5Dget_type(id); + if (disk_typeid == H5I_INVALID_HID) { + std::fprintf(stderr, "Failed to get type of attribute.\n"); + return -1; + } else if (H5Tdetect_class(disk_typeid, H5T_STRING) <= 0) { + std::fprintf(stderr, "Expected to be a string.\n"); + H5Tclose(disk_typeid); + return -1; + } + + htri_t is_variable_rslt = H5Tis_variable_str(disk_typeid); + if (is_variable_rslt < 0) { + std::fprintf(stderr, "Error in H5Tis_variable_str.\n"); + H5Tclose(disk_typeid); + return -1; + } + is_variable = (is_variable_rslt > 0); + + H5T_cset_t cset = H5Tget_cset(disk_typeid); + if (cset == H5T_CSET_ERROR) { + std::fprintf(stderr, "Error in H5Tget_cset.\n"); + H5Tclose(disk_typeid); + return -1; + } + uses_utf8_encoding = (cset == H5T_CSET_UTF8); + + H5Tclose(disk_typeid); + } + + // construct the disk_typeid that we use for describing `buffer` to the hdf5 + // library + hid_t memory_typeid = H5Tcopy(H5T_C_S1); + if (H5Tset_size(memory_typeid, is_variable ? H5T_VARIABLE : bufsz) < 0) { + fprintf(stderr, "error in H5Tset_size (for memory datatype)\n"); + H5Tclose(memory_typeid); + return -1; + } + if (uses_utf8_encoding && (H5Tset_cset(memory_typeid, H5T_CSET_UTF8) < 0)) { + fprintf(stderr, "error in H5Tset_size (for memory datatype)\n"); + H5Tclose(memory_typeid); + return -1; + } + if (H5Tset_strpad(memory_typeid, H5T_STR_NULLTERM) < 0) { + fprintf(stderr, "error in H5Tset_strpad (for memory datatype)\n"); + H5Tclose(memory_typeid); + return -1; + } + + // actually read the data + if (is_variable) { + char* tmp_str = nullptr; + herr_t status = + (is_attr) ? H5Aread(id, memory_typeid, static_cast(&tmp_str)) + : H5Dread(id, memory_typeid, H5S_ALL, H5S_ALL, H5P_DEFAULT, + static_cast(&tmp_str)); + if (status < 0) { + H5Tclose(memory_typeid); + if (tmp_str != nullptr) { + H5free_memory(tmp_str); + } + fprintf(stderr, "error in %s for a variable length string\n", + (is_attr) ? "H5Aread" : "H5Dread"); + return -1; + } + std::strncpy(buffer, tmp_str, static_cast(required_bufsz)); + } else { + herr_t status = (is_attr) ? H5Aread(id, memory_typeid, buffer) + : H5Dread(id, memory_typeid, H5S_ALL, H5S_ALL, + H5P_DEFAULT, buffer); + if (status < 0) { + H5Tclose(memory_typeid); + fprintf(stderr, "error in %s for a fixed length string\n", + (is_attr) ? "H5Aread" : "H5Dread"); + return -1; + } + } + H5Tclose(memory_typeid); + + if (uses_utf8_encoding) { + // if the string was stored on the disk with a utf-8 encoding, we need to + // confirm that the string *only* includes the subset of utf-8 code-points + // that are also valid characters are all compatible with ASCII. + if (!is_ascii_string(buffer, bufsz)) { + std::fprintf(stderr, + "Error: read a UTF-8 encoded string from an attribute that " + "contains non-ASCII code points\n"); + return -1; + } + } + + return required_bufsz; +} + +} // anonymous namespace + +int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, + char* buffer) { + return read_str_data_helper_(attr_id, true, bufsz, buffer); +} + +int grackle::impl::h5io::read_str_dataset(hid_t file_id, const char* dset_name, + int bufsz, char* buffer) { + if (dset_name == nullptr) { + std::fprintf(stderr, "dset_name is a nullptr"); + return GR_FAIL; + } + + hid_t dset_id = H5Dopen(file_id, dset_name); + if (dset_id == H5I_INVALID_HID) { + std::fprintf(stderr, "Failed to open dataset \"%s\".\n", dset_name); + return GR_FAIL; + } + + int out = read_str_data_helper_(dset_id, false, bufsz, buffer); + H5Dclose(dset_id); + return out; +} + +namespace { // stuff inside an anonymous namespace is local to this file + +/// Construct an ArrayShape instance from a HDF5 dataspace +grackle::impl::h5io::ArrayShape shape_from_space(hid_t space_id) { + grackle::impl::h5io::ArrayShape out; + out.ndim = -2; // set up an output value that denotes an error + if (space_id == H5I_INVALID_HID) { + return out; + } + + H5S_class_t space_type = H5Sget_simple_extent_type(space_id); + switch (space_type) { + case H5S_NULL: { + out.ndim = -1; + return out; + } + case H5S_SCALAR: { + out.ndim = 0; + return out; + } + case H5S_SIMPLE: { + int ndim = H5Sget_simple_extent_ndims(space_id); + hsize_t tmp[GRACKLE_CLOUDY_TABLE_MAX_DIMENSION]; + if (ndim <= 0) { // unclear what ndim of 0 would mean + std::fprintf(stderr, "error in call to H5Sget_simple_extent_ndims\n"); + return out; + } else if (ndim > GRACKLE_CLOUDY_TABLE_MAX_DIMENSION) { + std::fprintf( + stderr, + "hdf5 array has %d dimensions, exceeding the hardcoded maximum\n", + ndim); + return out; + } else if (ndim != H5Sget_simple_extent_dims(space_id, tmp, nullptr)) { + std::fprintf(stderr, "error in call to H5Sget_simple_extent_dims\n"); + return out; + } + + for (int i = 0; i < ndim; i++) { + if (tmp[i] > INT64_MAX) { + std::fprintf(stderr, "shape[%d] is too big\n", i); + return out; + } + out.shape[i] = static_cast(tmp[i]); + } + out.ndim = ndim; // make sure to do this last! + return out; + } + default: { + std::fprintf(stderr, "unexpected dataspace issue\n"); + return out; + } + } +} + +/// Construct an ArrayShape object that represents an invalid instance +/// +/// @note +/// This is useful for encoding that a function produced an error +grackle::impl::h5io::ArrayShape mk_invalid_array_shape() { + return shape_from_space(H5I_INVALID_HID); +} + +} // anonymous namespace + +grackle::impl::h5io::ArrayShape grackle::impl::h5io::read_dataset_shape( + hid_t file_id, const char* name) { + if (name == nullptr) { + std::fprintf(stderr, "name is a nullptr"); + return mk_invalid_array_shape(); + } + + hid_t dset_id = H5Dopen(file_id, name); + if (dset_id == H5I_INVALID_HID) { + std::fprintf(stderr, "Failed to open dataset \"%s\".\n", name); + return mk_invalid_array_shape(); + } + + hid_t space_id = H5Dget_space(dset_id); + H5Dclose(dset_id); + + ArrayShape out = shape_from_space(space_id); + H5Sclose(space_id); + + if (ArrayShape_is_null(out)) { + std::fprintf(stderr, "Error reading the dataspace of \"%s\".\n", name); + } + return out; +} + +/// read the dataset named dset_name from file_id into buffer +int grackle::impl::h5io::read_dataset( + hid_t file_id, const char* dset_name, double* buffer, + const grackle::impl::h5io::ArrayShape* expected_shape) { + if (dset_name == nullptr) { + std::fprintf(stderr, "dset_name is a nullptr"); + return GR_FAIL; + } + + hid_t dset_id = H5Dopen(file_id, dset_name); + if (dset_id == H5I_INVALID_HID) { + std::fprintf(stderr, "Failed to open dataset \"%s\".\n", dset_name); + return GR_FAIL; + } + + // we may perform an extra quick check about the dataset's shape + if (expected_shape != nullptr) { + hid_t space_id = H5Dget_space(dset_id); + ArrayShape actual_shape = shape_from_space(space_id); + H5Sclose(space_id); + + if (!ArrayShape_is_equal(*expected_shape, actual_shape)) { + H5Dclose(dset_id); + std::fprintf(stderr, "The \"%s\" dataset has an unexpected shape.\n", + dset_name); + return GR_FAIL; + } + } + + herr_t status = + H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); + H5Dclose(dset_id); + if (status < 0) { + std::fprintf(stderr, "Failed to read dataset \"%s\".\n", dset_name); + return GR_FAIL; + } + + return GR_SUCCESS; +} + +namespace { // stuff inside an anonymous namespace is local to this file + +/// @defgroup AttrNameRecorderGroup Attribute Recording Group +/// This is a collection of crude machinery used when we parse hdf5 dataset +/// attributes that specify Grid Table properties. This machinery is designed +/// to operate in 2 modes: +/// 1. A lightweight counter of attribute names +/// 2. A heavier weight mode that actually tracks alls attribute names. The +/// premise is to only use this mode (when we already know that there is an +/// error) in order to provide a descriptive error message. +///@{ +static constexpr int attr_name_recorder_capacity = + (GRACKLE_CLOUDY_TABLE_MAX_DIMENSION * 3 + 4); + +/// Acts as a crude tool for checking the set of attribute names +/// +/// @note +/// if we ever choose to embrace C++, it would be **MUCH** nicer to use +/// std::set +struct AttrNameRecorder { + bool only_count; + int length; + int capacity; + char* names[attr_name_recorder_capacity]; +}; + +AttrNameRecorder new_AttrNameRecorder(bool only_count) { + AttrNameRecorder out; + out.only_count = only_count; + out.capacity = attr_name_recorder_capacity; + out.length = 0; + for (int i = 0; i < attr_name_recorder_capacity; i++) { + out.names[i] = nullptr; + } + return out; +} + +/// acts as a destructor for all data contained within ptr +void drop_AttrNameRecorder(AttrNameRecorder* ptr) { + if (!ptr->only_count) { + for (int i = 0; i < ptr->length; i++) { + delete[] ptr->names[i]; + } + } + ptr->length = 0; +} + +/// records the attribute name +/// +/// @returns GR_SUCCESS if successful. Any other value denotes an issue +/// +/// @note +/// There are implicit assumptions that ptr and name are not nullptrs. It also +/// makes some sense to assume that name wasn't already recorded. +int AttrNameRecorder_record_name(AttrNameRecorder* ptr, const char* name) { + if (ptr->length == ptr->capacity) { + std::fprintf(stderr, + "encountered more attribute names than hardcoded maximum\n"); + return GR_FAIL; + } + + if (!ptr->only_count) { + std::size_t len_without_nullchr = std::strlen(name); + std::size_t buf_len = len_without_nullchr + 1; + char* buf = new char[buf_len]; + std::memcpy(buf, name, buf_len); + ptr->names[ptr->length] = buf; + } + + ptr->length++; + return GR_SUCCESS; +} + +int AttrNameRecorder_length(AttrNameRecorder* ptr) { return ptr->length; } + +/// writes a representation of all recorded attribute names to a character +/// string @p buffer. +/// +/// Just like ``std::snprintf``, this function writes at most ``buf_size - 1`` +/// characters and the resulting character string is terminated with a null +/// character. As with ``std::snprintf``, when @p buf_size is zero, nothing is +/// written, but the return value is still calculated. +/// +/// @param[in] buffer Pointer to the buffer where characters are written. This +/// can **only** be a nullptr if `buf_size` is 0. +/// @param[in] buf_size Specifies that `buf_size - 1` characters can be written +/// to `buffer`, in addition to the null terminator character. +/// @param[in] obj The object whose contents are being printed +/// +/// \returns If successful, returns number of characters that would be written +/// for a sufficiently large buffer. This value always excludes the terminating +/// null character. Otherwise, a negative value denotes an error. +/// +/// @note +/// Obviously, this could be simplified substantially if we fully embrace C++ +/// (e.g. we could simply return a `std::string` rather than mimic the +/// interface of `std::snprintf`) +int AttrNameRecorder_stringify_attr_names(char* buffer, std::size_t buf_size, + const AttrNameRecorder* obj) { + int length = obj->length; + + if (obj->only_count) { + return std::snprintf(buffer, buf_size, "{<%d unspecified attrs>}", length); + } + + // this is inefficient, but probably okay since its only for error reporting + + const int total_bufsz = (buf_size > static_cast(INT_MAX)) + ? INT_MAX + : static_cast(buf_size); + int total_offset = 0; // accumulated offset from start of buffer + + auto helper_fn = [&total_offset, total_bufsz, buffer](const char* s) -> bool { + bool write = total_bufsz > total_offset; + int rslt = std::snprintf((write) ? buffer + total_offset : nullptr, + (write) ? total_bufsz - total_offset : 0, "%s", s); + bool is_err = rslt < 0; + total_offset = rslt + ((is_err) ? 0 : total_offset); + fprintf(stderr, "component, total_offset: `%s`, %d\n", s, total_offset); + return is_err; + }; + + if (helper_fn("{")) { + return total_offset; + } + + for (int i = 0; i < length; i++) { + if (i > 0 && helper_fn(", ")) { + return total_offset; + } + if (helper_fn("\"") || helper_fn(obj->names[i]) || helper_fn("\"")) { + return total_offset; + } + } + + helper_fn("}"); + return total_offset; +} + +///@} + +/// helper function that calculates ArrayShape from the "Rank" and "Dimension" +/// attributes +/// +/// @param[in] dset_id The dataset identifier +/// @param[in] dset_name Name associated with dset_id (only used for producing +/// nicer error messages) +/// @param[inout] name_recorder Updated to track each attribute involved +/// in parsing a dataset's grid table properties. +/// +/// @returns the parsed array shape for the grid shape properties. The caller +/// should ensure that this function was successful by calling +/// ArrayShape_is_valid on the returned value. +grackle::impl::h5io::ArrayShape shape_from_grid_attrs( + hid_t dset_id, const char* dset_name, AttrNameRecorder* name_recorder) { + hid_t attr_id = H5Aopen_name(dset_id, "Rank"); + if (attr_id == H5I_INVALID_HID) { + std::fprintf(stderr, + "Failed to open \"Rank\" attribute of \"%s\" dataset.\n", + dset_name); + return mk_invalid_array_shape(); + } + long long rank; + if (H5Aread(attr_id, HDF5_I8, &rank) < 0) { + H5Aclose(attr_id); + std::fprintf(stderr, + "Failed to read \"Rank\" attribute of \"%s\" dataset.\n", + dset_name); + return mk_invalid_array_shape(); + } + H5Aclose(attr_id); + if ((rank < 0) || (GRACKLE_CLOUDY_TABLE_MAX_DIMENSION < rank)) { + H5Dclose(dset_id); + std::fprintf( + stderr, + "\"Rank\" attribute of \"%s\" dataset, %lld, is negative or exceeds " + "the hardcoded maximum, %d.\n", + dset_name, rank, GRACKLE_CLOUDY_TABLE_MAX_DIMENSION); + return mk_invalid_array_shape(); + } + if (AttrNameRecorder_record_name(name_recorder, "Rank") != GR_SUCCESS) { + std::fprintf( + stderr, + "issue recording access of \"Rank\" attribute of \"%s\" dataset.", + dset_name); + return mk_invalid_array_shape(); + } + + // read in the Dimension attribute + std::int64_t grid_dimensions[GRACKLE_CLOUDY_TABLE_MAX_DIMENSION]; + for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { + grid_dimensions[i] = -1; + } + attr_id = H5Aopen_name(dset_id, "Dimension"); + // todo: it would be nice to validate that the shape of dimension matches + // our expectation + if (attr_id == H5I_INVALID_HID) { + std::fprintf(stderr, + "Failed to open \"Dimension\" attribute in \"%s\" dataset.\n", + dset_name); + return mk_invalid_array_shape(); + } + if (H5Aread(attr_id, HDF5_I8, grid_dimensions) < 0) { + std::fprintf(stderr, + "Failed to read \"Dimension\" attribute in \"%s\" dataset.\n", + dset_name); + return mk_invalid_array_shape(); + } + H5Aclose(attr_id); + // a quick check + for (int i = 0; i < rank; i++) { + if (grid_dimensions[i] < 1) { + std::fprintf( + stderr, + "\"Dimension\" attribute of \"%s\" dataset has non-positive value\n", + dset_name); + return mk_invalid_array_shape(); + } + } + if (AttrNameRecorder_record_name(name_recorder, "Dimension") != GR_SUCCESS) { + std::fprintf( + stderr, + "issue recording access of \"Dimension\" attribute of \"%s\" dataset.", + dset_name); + return mk_invalid_array_shape(); + } + + // finally, let's format the output + // -> we use mk_invalid_array_shape to suppress compiler warnings that out + // may be uninitialized + grackle::impl::h5io::ArrayShape out = mk_invalid_array_shape(); + out.ndim = static_cast(rank); + for (int i = 0; i < out.ndim; i++) { + out.shape[i] = grid_dimensions[i]; + } + return out; +} + +static constexpr int max_attr_name_length = 64; + +/// updates axes with parsed parameters +/// +/// @param[in] dset_id The dataset identifier +/// @param[in] dset_name Name associated with dset_id (only used for producing +/// nicer error messages) +/// @param[in] grid_shape Specifies the pre-parsed grid shape +/// @param[out] axes Buffers theat will be updated with the grid axes +/// @param[inout] name_recorder Updated to track each attribute involved +/// in parsing a dataset's grid table properties. +/// +/// @returns GR_SUCCESS if successful. Any other value denotes an issue +int set_grid_axes_props(hid_t dset_id, const char* dset_name, + grackle::impl::h5io::ArrayShape grid_shape, + grackle::impl::h5io::GridTableAxis* axes, + AttrNameRecorder* name_recorder) { + using grackle::impl::h5io::read_str_attribute; + + bool legacy_mode = H5Aexists(dset_id, "Temperature") > 0; + + for (int i = 0; i < grid_shape.ndim; i++) { + // Step 1: fill axes[i].name with the name of the quantity that varies + // along axis i and fill val_attr_name with the name of the + // attribute containing the values of the quantitiy + char val_attr_name[max_attr_name_length]; + + if (legacy_mode && ((i + 1) == grid_shape.ndim)) { + axes[i].name = new char[12]; + std::snprintf(axes[i].name, 12, "Temperature"); + + std::snprintf(val_attr_name, max_attr_name_length, "Temperature"); + } else { + char tmp_attr_name[max_attr_name_length]; + std::snprintf(tmp_attr_name, max_attr_name_length, "Parameter%d_Name", + i + 1); + hid_t attr_id = H5Aopen_name(dset_id, tmp_attr_name); + + if (attr_id == H5I_INVALID_HID) { + std::fprintf(stderr, + "Failed to open \"%s\" attribute of \"%s\" dataset.\n", + tmp_attr_name, dset_name); + return GR_FAIL; + } + + int min_buf_length = read_str_attribute(attr_id, 0, nullptr); + if (min_buf_length < 0) { + std::fprintf( + stderr, + "can't determine string size held by \"%s\" attr of \"%s\" dset.\n", + tmp_attr_name, dset_name); + return GR_FAIL; + } + axes[i].name = new char[min_buf_length]; + if (read_str_attribute(attr_id, min_buf_length, axes[i].name) < 0) { + std::fprintf(stderr, + "error loading the \"%s\" attr from \"%s\" dataset\n", + tmp_attr_name, dset_name); + return GR_FAIL; + } + H5Aclose(attr_id); + + if (AttrNameRecorder_record_name(name_recorder, tmp_attr_name) != + GR_SUCCESS) { + std::fprintf(stderr, + "issue recording access of \"%s\" attr of \"%s\" dataset", + tmp_attr_name, dset_name); + return GR_FAIL; + } + + std::snprintf(val_attr_name, max_attr_name_length, "Parameter%d", i + 1); + } + + // step 2: load the values associated with axis i and store them in + // axes[i].values + hid_t attr_id = H5Aopen_name(dset_id, val_attr_name); + if (attr_id == H5I_INVALID_HID) { + std::fprintf(stderr, "Failed to open \"%s\" attr of \"%s\" dataset.\n", + val_attr_name, dset_name); + return GR_FAIL; + } + + hid_t space_id = H5Aget_space(attr_id); + grackle::impl::h5io::ArrayShape axis_shape = shape_from_space(space_id); + H5Sclose(space_id); + if (axis_shape.ndim != 1) { + std::fprintf(stderr, "The \"%s\" dataset's \"%s\" attr isn't 1D\n", + val_attr_name, dset_name); + return GR_FAIL; + } else if (axis_shape.shape[0] != grid_shape.shape[i]) { + std::fprintf( + stderr, + "The \"%s\" dataset's \"Dimension\" attr suggests that axis %d has " + "%lld values. The \"%s\" attribute actually %lld values\n", + dset_name, i, static_cast(grid_shape.shape[i]), + val_attr_name, static_cast(axis_shape.shape[0])); + return GR_FAIL; + } + + axes[i].values = new double[grid_shape.shape[i]]; + if (H5Aread(attr_id, HDF5_R8, axes[i].values) < 0) { + H5Aclose(attr_id); + std::fprintf(stderr, "failed to read \"%s\" attr of \"%s\" dataset.\n", + val_attr_name, dset_name); + return GR_FAIL; + } + H5Aclose(attr_id); + if (AttrNameRecorder_record_name(name_recorder, val_attr_name) != + GR_SUCCESS) { + std::fprintf(stderr, + "issue recording access of \"%s\" attr of \"%s\" dataset", + val_attr_name, dset_name); + return GR_FAIL; + } + } + return GR_SUCCESS; +} + +grackle::impl::h5io::GridTableProps mk_invalid_GridTableProps() { + grackle::impl::h5io::GridTableProps out; + out.table_shape = mk_invalid_array_shape(); + for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { + out.axes[i].name = nullptr; + out.axes[i].values = nullptr; + } + return out; +} + +/// helper function that parses the GridTableProps from dataset attributes +/// +/// @param[in] dset_id The dataset identifier +/// @param[in] dset_name Name associated with dset_id (only used for producing +/// nicer error messages) +/// @param[inout] name_recorder Updated to track each attribute involved +/// in parsing a dataset's grid table properties. +/// +/// @returns Returns the appropriate GridTableProps object. The caller should +/// use the GridTableProps_is_valid function to confirm that the function +/// was successful. +grackle::impl::h5io::GridTableProps parse_GridTableProps_helper( + hid_t dset_id, const char* dset_name, AttrNameRecorder* name_recorder) { + // setup the output object so that we're always prepared to return an object + // - the default object is constructed to denote a failure + grackle::impl::h5io::GridTableProps out = mk_invalid_GridTableProps(); + + // if optional Description attribute is present, record that we accessed it + // (this is done purely for error-handling purposes). + if ((H5Aexists(dset_id, "Description") > 0) && + (AttrNameRecorder_record_name(name_recorder, "Description") != + GR_SUCCESS)) { + return out; + } + + // infer the shape of the table from the attributes + grackle::impl::h5io::ArrayShape inferred_shape = + shape_from_grid_attrs(dset_id, dset_name, name_recorder); + if (!ArrayShape_is_valid(inferred_shape)) { + return out; + } + + // parse the quantities along each axis + if (set_grid_axes_props(dset_id, dset_name, inferred_shape, out.axes, + name_recorder) != GR_SUCCESS) { + drop_GridTableProps(&out); + return out; + } + + out.table_shape = inferred_shape; // we intentionally do this as the very + // last step! + return out; +} + +int get_num_attrs(hid_t dset_id, const char* dset_name) { +#if H5_VERSION_LE(1, 10, 2) + return -2; +#else +#if H5_VERSION_GE(1, 12, 0) + H5O_info2_t info; + herr_t status = H5Oget_info3(dset_id, &info, H5O_INFO_NUM_ATTRS); +#else + H5O_info_t info; + herr_t status = H5Oget_info2(dset_id, &info, H5O_INFO_NUM_ATTRS); +#endif + if (status < 0) { + std::fprintf(stderr, "Can't get num_attrs for \"%s\" dataset.\n", + dset_name); + return -1; + } + if (info.num_attrs > INT_MAX) { + std::fprintf(stderr, "\"%s\" dataset has too many attrs.\n", dset_name); + return -1; + } else { + return static_cast(info.num_attrs); + } +#endif +} + +} // anonymous namespace + +grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( + hid_t file_id, const char* dset_name) { + if (dset_name == nullptr) { // sanity check! + std::fprintf(stderr, "dset_name is a nullptr"); + return mk_invalid_GridTableProps(); + } + + hid_t dset_id = H5Dopen(file_id, dset_name); + if (dset_id == H5I_INVALID_HID) { + std::fprintf(stderr, "Can't open \"%s\" dataset.\n", dset_name); + return mk_invalid_GridTableProps(); + } + + // construct object to count accessed attributes + AttrNameRecorder attr_counter = new_AttrNameRecorder(true); + // actually parse GridTableProps + GridTableProps out = + parse_GridTableProps_helper(dset_id, dset_name, &attr_counter); + // get the attribute count and cleanup the counter + int total_accessed_attrs_count = AttrNameRecorder_length(&attr_counter); + drop_AttrNameRecorder(&attr_counter); + + if (!GridTableProps_is_valid(out)) { + H5Dclose(dset_id); + // parse_GridTableProps_helper already printed appropriate errors messages + return out; + } + + // now, we will perform a check validating that dataset doesn't have + // unexpected attributes (i.e. we want to avoid the situation where the + // creator expects new attributes to introduce some behavior that Grackle + // doesn't know about) + int num_attrs = get_num_attrs(dset_id, dset_name); + if (num_attrs == -2) { + num_attrs = total_accessed_attrs_count; + } else if (num_attrs == -1) { + H5Dclose(dset_id); + drop_GridTableProps(&out); + // get_num_attrs already printed error messages in this case + return out; + } else if (num_attrs != total_accessed_attrs_count) { + drop_GridTableProps(&out); + // to provide a detailed error message that will make it straight-forward + // to debug the underlying problem, we are going to call + // parse_GridTableProps_helper, but this time, we are going to actually + // record every accessed name. + AttrNameRecorder attr_name_recorder = new_AttrNameRecorder(false); + // let's parse GridTableProps (again) + GridTableProps tmp = + parse_GridTableProps_helper(dset_id, dset_name, &attr_name_recorder); + + int bufsz_without_nullchr = + AttrNameRecorder_stringify_attr_names(nullptr, 0, &attr_name_recorder); + int bufsz = bufsz_without_nullchr + 1; + char* stringified_attr_list = nullptr; + if (bufsz > 0) { + stringified_attr_list = new char[bufsz]; + if (AttrNameRecorder_stringify_attr_names(stringified_attr_list, + static_cast(bufsz), + &attr_name_recorder) < 0) { + delete[] stringified_attr_list; + stringified_attr_list = nullptr; + } + } + + if (stringified_attr_list == nullptr) { + std::fprintf( + stderr, + "encountered error while stringifying list of recognized attrs for " + "the \"%s\" dataset\n", + dset_name); + } else { + std::fprintf( + stderr, + "while parsing the %dD grid table properties from attributes of the " + "\"%s\" dataset, encountered (one or more) unrecognized attributes. " + "Recognized attributes include: %s\n", + tmp.table_shape.ndim, dset_name, stringified_attr_list); + delete[] stringified_attr_list; + } + drop_GridTableProps(&tmp); + drop_AttrNameRecorder(&attr_name_recorder); + H5Dclose(dset_id); + return mk_invalid_GridTableProps(); + } + + // check for consistency between the GridTableProps and the dataset shape + hid_t space_id = H5Dget_space(dset_id); + ArrayShape actual_shape = shape_from_space(space_id); + H5Sclose(space_id); + if (!ArrayShape_is_null(actual_shape) && + !ArrayShape_is_equal(out.table_shape, actual_shape)) { + drop_GridTableProps(&out); + H5Dclose(dset_id); + std::fprintf( + stderr, + "The grid table properties parsed from the attributes of the \"%s\" " + "dataset are inconsistent with the dataset's shape.\n", + dset_name); + return mk_invalid_GridTableProps(); + } + + H5Dclose(dset_id); + return out; +} + +bool grackle::impl::h5io::GridTableProps_is_equal( + grackle::impl::h5io::GridTableProps props_a, + grackle::impl::h5io::GridTableProps props_b) { + if (!ArrayShape_is_equal(props_a.table_shape, props_b.table_shape)) { + return false; + } + for (int i = 0; i < props_a.table_shape.ndim; i++) { + if (std::strcmp(props_a.axes[i].name, props_b.axes[i].name) != 0) { + return false; + } + std::int64_t length = props_a.table_shape.shape[i]; + std::int64_t n_equal = 0; + for (std::int64_t j = 0; j < length; j++) { + bool is_equal = (props_a.axes[i].values[j] == props_b.axes[i].values[j]); + n_equal += static_cast(is_equal); + } + if (n_equal != length) { + return false; + } + } + return true; +} + +int grackle::impl::h5io::assert_has_consistent_GridTableProps( + hid_t file_id, const char* dset_name, + grackle::impl::h5io::GridTableProps expected) { + // the current implementation is crude (we may be able to do better) + + GridTableProps actual = parse_GridTableProps(file_id, dset_name); + + if (!GridTableProps_is_valid(actual)) { + drop_GridTableProps(&actual); + fprintf(stderr, + "Error constructing the grid properties for the \"%s\" dataset\n", + dset_name); + return GR_FAIL; + } + + bool is_equal = GridTableProps_is_equal(actual, expected); + drop_GridTableProps(&actual); + if (!is_equal) { + fprintf(stderr, + "the \"%s\" dataset doesn't have the expected grid properties\n", + dset_name); + return GR_FAIL; + } + + return GR_SUCCESS; +} diff --git a/src/clib/support/h5io.hpp b/src/clib/support/h5io.hpp new file mode 100644 index 000000000..80748ff64 --- /dev/null +++ b/src/clib/support/h5io.hpp @@ -0,0 +1,205 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Declares functions to help read from hdf5 files +/// +//===----------------------------------------------------------------------===// + +#ifndef UTILS_H5IO_HPP +#define UTILS_H5IO_HPP + +#include + +#include "hdf5.h" +#include "grackle.h" + +namespace grackle::impl::h5io { + +/// copies the string encoded in the specified hdf5 dataset into ``buffer`` as +/// a null-terminated string, and returns ``min_req_bufsz`` (if successful). +/// +/// @param[in] attr_id Attribute identifier +/// @param[in] bufsz Number of ascii characters (including the null terminating +/// character) that can be written to @p buffer +/// @param[out] buffer Pointer to the buffer where characters are written. This +/// can **only** be a nullptr if @p bufsz is 0. +/// +/// @returns If successful, returns ``min_req_bufsz`` (see below). Otherwise, +/// returns a negative value. +/// +/// ``min_req_bufsz`` is the minimum required @p bufsz that this function must +/// receive for it to attempt to load the string. +/// - this is the maximum length of the string (including the null character). +/// Thus, after succesfully calling this function +/// ``std::strlen(buffer) + 1 <= min_req_bufsz``. +/// - this function's behavoir is described in terms of ``min_req_bufsz``, +/// rather than the exact required buffer length because the exact length +/// can't be determined without loading the buffer. +/// +/// This function fails if @p bufsz is smaller than ``min_req_bufsz``, unless +/// @p bufsz is zero. In that case, nothing is written to @p buffer and the +/// returns ``min_req_bufsz``. The function reports an error if the user tries +/// to reads a utf8-encoded string that contains non-ASCII characters. +/// +/// @note +/// If we are more willing to embrace C++, we could return a std::string or +/// std::vector rather than requiring a pre-allocated buffer +int read_str_attribute(hid_t attr_id, int bufsz, char* buffer); + +/// copies the string encoded in the specified hdf5 dataset into ``buffer`` as +/// a null-terminated string, and returns ``min_req_bufsz`` (if successful). +/// +/// @param[in] file_id File identifier +/// @param[in] dset_name The name of the dataset to read attributes from. +/// @param[in] bufsz Number of ascii characters (including the null terminating +/// character) that can be written to @p buffer +/// @param[out] buffer Pointer to the buffer where characters are written. This +/// can **only** be a nullptr if @p bufsz is 0. +/// +/// @returns If successful, returns ``min_req_bufsz`` (see below). Otherwise, +/// returns a negative value. +/// +/// ``min_req_bufsz`` is the minimum required @p bufsz that this function must +/// receive for it to attempt to load the string. +/// - this is the maximum length of the string (including the null character). +/// Thus, after succesfully calling this function +/// ``std::strlen(buffer) + 1 <= min_req_bufsz``. +/// - this function's behavoir is described in terms of ``min_req_bufsz``, +/// rather than the exact required buffer length because the exact length +/// can't be determined without loading the buffer. +/// +/// This function fails if @p bufsz is smaller than ``min_req_bufsz``, unless +/// @p bufsz is zero. In that case, nothing is written to @p buffer and the +/// returns ``min_req_bufsz``. The function reports an error if the user tries +/// to reads a utf8-encoded string that contains non-ASCII characters. +/// +/// @note +/// If we are more willing to embrace C++, we could return a std::string or +/// std::vector rather than requiring a pre-allocated buffer +/// +/// @note +/// The choice to accept @p file_id and @p dset_name, rather than an already +/// open dataset identifier, was made for consistency with the interfaces of +/// read_dataset and read_dataset_shape. However, the choice is a little +/// "clunky" because this function is usually called twice (once to query +/// ``min_req_bufsz`` and once to load the string). +int read_str_dataset(hid_t file_id, const char* dset_name, int bufsz, + char* buffer); + +/// represents a contiguous array shape +/// +/// @note +/// An ndim of -1 corresponds to a null dataset (i.e. ``H5S_NULL``). Any other +/// negative value denotes an invalid shape. An ndim of 0 denotes a scalar +struct ArrayShape { + int ndim; + std::int64_t shape[GRACKLE_CLOUDY_TABLE_MAX_DIMENSION]; +}; + +/// checks whether shape is valid +inline bool ArrayShape_is_valid(ArrayShape shape) { return shape.ndim >= -1; } + +/// checks whether shape refers to a scalar +inline bool ArrayShape_is_scalar(ArrayShape shape) { return shape.ndim == 0; } + +/// checks whether shape is null +inline bool ArrayShape_is_null(ArrayShape shape) { return shape.ndim == -1; } + +/// calculates the total number of elements in the array +inline std::int64_t ArrayShape_elem_count(ArrayShape shape) { + if (shape.ndim < 0) { + return -1; + } else { // this works even if shape.ndim is 0 + std::int64_t product = 1; + for (int i = 0; i < shape.ndim; i++) { + product *= shape.shape[i]; + } + return product; + } +} + +/// checks whether shape_a and shape_b are the same +/// +/// returns false if the either shape is invalid +inline bool ArrayShape_is_equal(ArrayShape shape_a, ArrayShape shape_b) { + if ((!ArrayShape_is_valid(shape_a)) || (shape_a.ndim != shape_b.ndim)) { + return false; + } + for (int i = 0; i < shape_a.ndim; i++) { + if (shape_a.shape[i] != shape_b.shape[i]) { + return false; + } + } + return true; +} + +/// load the shape of the dataset +ArrayShape read_dataset_shape(hid_t file_id, const char* dset_name); + +/// read the dataset named dset_name from file_id into buffer +/// +/// When expected_shape is provided, the dataset's shape is checked before +/// the data is read +int read_dataset(hid_t file_id, const char* dset_name, double* buffer, + const ArrayShape* expected_shape = nullptr); + +struct GridTableAxis { + char* name; + double* values; +}; + +/// Used to represent properties of an interpolation table +struct GridTableProps { + ArrayShape table_shape; + GridTableAxis axes[GRACKLE_CLOUDY_TABLE_MAX_DIMENSION]; +}; + +inline bool GridTableProps_is_valid(GridTableProps grid_props) { + return ArrayShape_is_valid(grid_props.table_shape); +} + +/// acts as a destructor for the contents within ptr +inline void drop_GridTableProps(GridTableProps* ptr) { + for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { + if (ptr->axes[i].name != nullptr) { + delete[] ptr->axes[i].name; + } + if (ptr->axes[i].values != nullptr) { + delete[] ptr->axes[i].values; + } + ptr->axes[i].name = nullptr; + ptr->axes[i].values = nullptr; + } +} + +/// parses the GridTableProps from dataset attributes +/// +/// @param[in] file_id File identifier +/// @param[in] dset_name The name of the dataset to read attributes from. +/// +/// @returns Returns the appropriate GridTableProps object. The caller should +/// use the GridTableProps_is_valid function to confirm that the function +/// was successful. +GridTableProps parse_GridTableProps(hid_t file_id, const char* dset_name); + +/// checks whether props_a and props_b hold equivalent values +/// +/// @returns true if objects are equal and false in all other cases (including +/// when either object is invalid) +bool GridTableProps_is_equal(GridTableProps props_a, GridTableProps props_b); + +/// checks whether the specified dataset has consistent grid properties +/// +/// @returns GR_SUCCESS if the specified dataset has equivalent properties and +/// a different value in all other cases +int assert_has_consistent_GridTableProps(hid_t file_id, const char* dset_name, + GridTableProps expected); + +} // namespace grackle::impl::h5io + +#endif /* UTILS_H5IO_HPP */