From d5ff24e84282f79166f60371766451da210fb6d6 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 13:30:49 -0500 Subject: [PATCH 01/38] Initial conversion initialize_cloudy_data.{c->cpp} This commit includes a few more changes than I initially wanted --- .clang-format-ignore | 2 +- src/clib/CMakeLists.txt | 2 +- src/clib/initialize_chemistry_data.cpp | 2 +- ...oudy_data.c => initialize_cloudy_data.cpp} | 191 +++++++++--------- ...oudy_data.h => initialize_cloudy_data.hpp} | 3 - 5 files changed, 103 insertions(+), 97 deletions(-) rename src/clib/{initialize_cloudy_data.c => initialize_cloudy_data.cpp} (55%) rename src/clib/{initialize_cloudy_data.h => initialize_cloudy_data.hpp} (89%) diff --git a/.clang-format-ignore b/.clang-format-ignore index 36fa238b6..900c1bd4b 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -32,7 +32,7 @@ src/clib/index_helper.h src/clib/init_misc_species_cool_rates.cpp src/clib/initialize_UVbackground_data.c 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..68858ad40 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -90,7 +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 @@ -114,6 +113,7 @@ 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 internal_types.cpp internal_types.hpp diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 294dbf7d6..ab18276a3 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -21,7 +21,7 @@ #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" diff --git a/src/clib/initialize_cloudy_data.c b/src/clib/initialize_cloudy_data.cpp similarity index 55% rename from src/clib/initialize_cloudy_data.c rename to src/clib/initialize_cloudy_data.cpp index 88ab57ada..1140d0e67 100644 --- a/src/clib/initialize_cloudy_data.c +++ b/src/clib/initialize_cloudy_data.cpp @@ -12,32 +12,34 @@ ************************************************************************/ #include -#include -#include +#include // std::fprintf, std::sprintf, stderr, stdio +#include // std::log10, std::pow +#include // std::strcmp #include "hdf5.h" #include "grackle.h" #include "grackle_macros.h" -#include "initialize_cloudy_data.h" +#include "initialize_cloudy_data.hpp" -#define SMALL_LOG_VALUE -99.0 +#define SMALL_LOG_VALUE (-99.0) +namespace { // stuff inside of an anonymous namespace is read-only -/** - * Initializes an empty #cloudy_data struct with zeros and NULLs. - */ +/// 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] = NULL; + my_cloudy->grid_parameters[i] = nullptr; } - my_cloudy->heating_data = NULL; - my_cloudy->cooling_data = NULL; - my_cloudy->mmw_data = NULL; + my_cloudy->heating_data = nullptr; + my_cloudy->cooling_data = nullptr; + my_cloudy->mmw_data = nullptr; my_cloudy->data_size = 0LL; } +} // anonymous namespace + // initialize cloudy cooling data int initialize_cloudy_data(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, @@ -54,11 +56,11 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, // 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 (read_data == 0) { return GR_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); + 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. */ @@ -72,17 +74,18 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, 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); + 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 * - POW(my_units->a_value * my_units->a_units, 3); + std::pow(my_units->a_value * my_units->a_units, 3.0); double mh = 1.67e-24; - double CoolUnit = (POW(my_units->a_units,5) * POW(xbase1,2) * POW(mh,2)) / - (POW(tbase1,3) * dbase1); + 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, dset_id, attr_id; @@ -95,64 +98,65 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, if (H5Aexists(file_id, "old_style")) { my_rates->cloudy_data_new = 0; if (grackle_verbose) - fprintf(stdout, "Loading old-style Cloudy tables.\n"); + std::fprintf(stdout, "Loading old-style Cloudy tables.\n"); } // Open cooling dataset and get grid dimensions. - sprintf(parameter_name, "/CoolingRates/%s/Cooling", group_name); + std::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", + std::fprintf(stderr,"Can't open Cooling in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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; + std::fprintf(stderr,"Failed to open Rank attribute in Cooling dataset.\n"); + return GR_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; + std::fprintf(stderr,"Failed to read Rank attribute in Cooling dataset.\n"); + return GR_FAIL; } my_cloudy->grid_rank = (long long) temp_int; if (grackle_verbose) - fprintf(stdout,"Cloudy cooling grid rank: %lld.\n", my_cloudy->grid_rank); + std::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; + std::fprintf(stderr,"Failed to close Rank attribute in Cooling dataset.\n"); + return GR_FAIL; } // Grid dimension. - temp_int_arr = malloc(my_cloudy->grid_rank * sizeof(long long)); + temp_int_arr = (long long*)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; + std::fprintf(stderr,"Failed to open Dimension attribute in Cooling dataset.\n"); + return GR_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; + std::fprintf(stderr,"Failed to read Dimension attribute in Cooling dataset.\n"); + return GR_FAIL; + } + if (grackle_verbose) { + std::fprintf(stdout, "Cloudy cooling grid dimensions:"); } - 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]); + std::fprintf(stdout," %lld", my_cloudy->grid_dimension[q]); } if (grackle_verbose) - fprintf(stdout,".\n"); + std::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; + std::fprintf(stderr,"Failed to close Dimension attribute in Cooling dataset.\n"); + return GR_FAIL; } free(temp_int_arr); @@ -160,28 +164,28 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, for (q = 0;q < my_cloudy->grid_rank;q++) { if (q < my_cloudy->grid_rank - 1) { - sprintf(parameter_name,"Parameter%lld",(q+1)); + std::sprintf(parameter_name,"Parameter%lld",(q+1)); } else { - sprintf(parameter_name,"Temperature"); + std::sprintf(parameter_name,"Temperature"); } - temp_data = malloc(my_cloudy->grid_dimension[q] * sizeof(double)); + temp_data = (double*) 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", + std::fprintf(stderr,"Failed to open %s attribute in Cooling dataset.\n", parameter_name); - return FAIL; + return GR_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; + std::fprintf(stderr,"Failed to read %s attribute in Cooling dataset.\n", + parameter_name); + return GR_FAIL; } - my_cloudy->grid_parameters[q] = malloc(my_cloudy->grid_dimension[q] * + my_cloudy->grid_parameters[q] = (double*) 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) { @@ -189,20 +193,23 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, } else { // convert temeperature to log - my_cloudy->grid_parameters[q][w] = (double) log10(temp_data[w]); + my_cloudy->grid_parameters[q][w] = (double) std::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]); + if (grackle_verbose) { + std::fprintf(stdout, + "%s: %g to %g (%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", + std::fprintf(stderr,"Failed to close %s attribute in Cooling dataset.\n", parameter_name); - return FAIL; + return GR_FAIL; } free(temp_data); @@ -213,109 +220,111 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, 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)); + temp_data = (double*) 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"); + std::fprintf(stdout,"Reading Cloudy Cooling dataset.\n"); if (status == h5_error) { - fprintf(stderr,"Failed to read Cooling dataset.\n"); - return FAIL; + std::fprintf(stderr,"Failed to read Cooling dataset.\n"); + return GR_FAIL; } - my_cloudy->cooling_data = malloc(my_cloudy->data_size * sizeof(double)); + my_cloudy->cooling_data = (double*) 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; + (double) std::log10(temp_data[q]) : (double) SMALL_LOG_VALUE; // Convert to code units. - my_cloudy->cooling_data[q] -= log10(CoolUnit); + my_cloudy->cooling_data[q] -= std::log10(CoolUnit); } free(temp_data); status = H5Dclose(dset_id); if (status == h5_error) { - fprintf(stderr,"Failed to close Cooling dataset.\n"); - return FAIL; + std::fprintf(stderr,"Failed to close Cooling dataset.\n"); + return GR_FAIL; } // Read Heating data. if (my_chemistry->UVbackground == 1) { - temp_data = malloc(my_cloudy->data_size * sizeof(double)); + temp_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); - sprintf(parameter_name, "/CoolingRates/%s/Heating", group_name); + std::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", + std::fprintf(stderr,"Can't open Heating in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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"); + std::fprintf(stdout,"Reading Cloudy Heating dataset.\n"); if (status == h5_error) { - fprintf(stderr,"Failed to read Heating dataset.\n"); - return FAIL; + std::fprintf(stderr,"Failed to read Heating dataset.\n"); + return GR_FAIL; } - my_cloudy->heating_data = malloc(my_cloudy->data_size * sizeof(double)); + my_cloudy->heating_data = (double*) 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; + (double) std::log10(temp_data[q]) : (double) SMALL_LOG_VALUE; // Convert to code units. - my_cloudy->heating_data[q] -= log10(CoolUnit); + my_cloudy->heating_data[q] -= std::log10(CoolUnit); } free(temp_data); status = H5Dclose(dset_id); if (status == h5_error) { - fprintf(stderr,"Failed to close Heating dataset.\n"); - return FAIL; + std::fprintf(stderr,"Failed to close Heating dataset.\n"); + return GR_FAIL; } } // Read MMW data. if (my_chemistry->primordial_chemistry == 0 && - strcmp(group_name, "Primordial") == 0) { + std::strcmp(group_name, "Primordial") == 0) { - my_cloudy->mmw_data = malloc(my_cloudy->data_size * sizeof(double)); + my_cloudy->mmw_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); - sprintf(parameter_name, "/CoolingRates/%s/MMW", group_name); + std::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", + std::fprintf(stderr,"Can't open MMW in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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"); + std::fprintf(stdout,"Reading Cloudy MMW dataset.\n"); if (status == h5_error) { - fprintf(stderr,"Failed to read MMW dataset.\n"); - return FAIL; + std::fprintf(stderr,"Failed to read MMW dataset.\n"); + return GR_FAIL; } status = H5Dclose(dset_id); if (status == h5_error) { - fprintf(stderr,"Failed to close MMW dataset.\n"); - return FAIL; + std::fprintf(stderr,"Failed to close MMW dataset.\n"); + return GR_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", + std::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 GR_FAIL; } - return SUCCESS; + return GR_SUCCESS; } int free_cloudy_data(cloudy_data *my_cloudy, chemistry_data *my_chemistry, int primordial) { diff --git a/src/clib/initialize_cloudy_data.h b/src/clib/initialize_cloudy_data.hpp similarity index 89% rename from src/clib/initialize_cloudy_data.h rename to src/clib/initialize_cloudy_data.hpp index 4159b8ad0..fc6301493 100644 --- a/src/clib/initialize_cloudy_data.h +++ b/src/clib/initialize_cloudy_data.hpp @@ -19,9 +19,6 @@ extern "C" { #endif /* __cplusplus */ -/// Initializes an empty #cloudy_data struct with zeros and NULLs. -void initialize_empty_cloudy_data_struct(cloudy_data* my_cloudy); - // initialize cloudy cooling data int initialize_cloudy_data(chemistry_data* my_chemistry, chemistry_data_storage* my_rates, From ec0b5c6d6e2ba25f1112d60e0064252915b3608b Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 13:34:51 -0500 Subject: [PATCH 02/38] a small tweak --- src/clib/initialize_cloudy_data.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 1140d0e67..d1ed0a7aa 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -47,7 +47,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, code_units *my_units, int read_data) { - long long q, w; + long long w; double *temp_data; long long temp_int; long long *temp_int_arr; @@ -146,7 +146,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, if (grackle_verbose) { std::fprintf(stdout, "Cloudy cooling grid dimensions:"); } - for (q = 0;q < my_cloudy->grid_rank;q++) { + for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { my_cloudy->grid_dimension[q] = (long long) temp_int_arr[q]; if (grackle_verbose) std::fprintf(stdout," %lld", my_cloudy->grid_dimension[q]); @@ -161,7 +161,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, free(temp_int_arr); // Grid parameters. - for (q = 0;q < my_cloudy->grid_rank;q++) { + for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { if (q < my_cloudy->grid_rank - 1) { std::sprintf(parameter_name,"Parameter%lld",(q+1)); @@ -187,7 +187,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, my_cloudy->grid_parameters[q] = (double*) malloc(my_cloudy->grid_dimension[q] * sizeof(double)); - for (w = 0;w < my_cloudy->grid_dimension[q];w++) { + for (long long w = 0LL; w < my_cloudy->grid_dimension[q]; w++) { if (q < my_cloudy->grid_rank - 1) { my_cloudy->grid_parameters[q][w] = (double) temp_data[w]; } @@ -217,7 +217,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, // Read Cooling data. my_cloudy->data_size = 1; - for (q = 0;q < my_cloudy->grid_rank;q++) { + for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { my_cloudy->data_size *= my_cloudy->grid_dimension[q]; } temp_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); @@ -231,7 +231,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, } my_cloudy->cooling_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); - for (q = 0;q < my_cloudy->data_size;q++) { + for (long long q = 0LL; q < my_cloudy->data_size; q++) { my_cloudy->cooling_data[q] = temp_data[q] > 0 ? (double) std::log10(temp_data[q]) : (double) SMALL_LOG_VALUE; @@ -268,7 +268,7 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, } my_cloudy->heating_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); - for (q = 0;q < my_cloudy->data_size;q++) { + for (long long q = 0LL; q < my_cloudy->data_size; q++) { my_cloudy->heating_data[q] = temp_data[q] > 0 ? (double) std::log10(temp_data[q]) : (double) SMALL_LOG_VALUE; From 5435b76302e00b995c55b324d36af50412b9c660 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 13:46:44 -0500 Subject: [PATCH 03/38] Replace sprintf with snprintf This was spurred on by deprecation warnings included on macOS's implementation of the standard library --- src/clib/initialize_cloudy_data.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index d1ed0a7aa..cf6ff6734 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -12,7 +12,7 @@ ************************************************************************/ #include -#include // std::fprintf, std::sprintf, stderr, stdio +#include // std::fprintf, std::snprintf, stderr, stdio #include // std::log10, std::pow #include // std::strcmp #include "hdf5.h" @@ -47,11 +47,11 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, code_units *my_units, int read_data) { - long long w; double *temp_data; long long temp_int; long long *temp_int_arr; char parameter_name[MAX_LINE_LENGTH]; + const std::size_t pname_bufsize = static_cast(MAX_LINE_LENGTH); // Initialize things (to the null-state) even if cloudy cooling is not used. initialize_empty_cloudy_data_struct(my_cloudy); @@ -103,7 +103,8 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, // Open cooling dataset and get grid dimensions. - std::sprintf(parameter_name, "/CoolingRates/%s/Cooling", group_name); + std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/Cooling", + group_name); dset_id = H5Dopen(file_id, parameter_name); if (dset_id == h5_error) { std::fprintf(stderr,"Can't open Cooling in %s.\n", @@ -164,10 +165,10 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { if (q < my_cloudy->grid_rank - 1) { - std::sprintf(parameter_name,"Parameter%lld",(q+1)); + std::snprintf(parameter_name, pname_bufsize, "Parameter%lld",(q+1)); } else { - std::sprintf(parameter_name,"Temperature"); + std::snprintf(parameter_name, pname_bufsize, "Temperature"); } temp_data = (double*) malloc(my_cloudy->grid_dimension[q] * sizeof(double)); @@ -251,7 +252,8 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, temp_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); - std::sprintf(parameter_name, "/CoolingRates/%s/Heating", group_name); + std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/Heating", + group_name); dset_id = H5Dopen(file_id, parameter_name); if (dset_id == h5_error) { std::fprintf(stderr,"Can't open Heating in %s.\n", @@ -290,7 +292,8 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, my_cloudy->mmw_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); - std::sprintf(parameter_name, "/CoolingRates/%s/MMW", group_name); + std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/MMW", + group_name); dset_id = H5Dopen(file_id, parameter_name); if (dset_id == h5_error) { std::fprintf(stderr,"Can't open MMW in %s.\n", From 978187c38fc6d2beefdc6d8bc646d8c10c733c2c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 13:57:30 -0500 Subject: [PATCH 04/38] mv `(initialize|free)_cloudy_data` into namespace --- src/clib/initialize_chemistry_data.cpp | 18 ++++++++++-------- src/clib/initialize_cloudy_data.cpp | 15 ++++++++++----- src/clib/initialize_cloudy_data.hpp | 8 ++------ 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index ab18276a3..37fa313f9 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -429,18 +429,18 @@ 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; } @@ -593,8 +593,10 @@ 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::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_FREE(my_rates->UVbackground_table.z); GRACKLE_FREE(my_rates->UVbackground_table.k24); diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index cf6ff6734..9c9d0f9a2 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -40,11 +40,13 @@ void initialize_empty_cloudy_data_struct(cloudy_data *my_cloudy) } // anonymous namespace + // 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) +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) { double *temp_data; @@ -330,7 +332,9 @@ int initialize_cloudy_data(chemistry_data *my_chemistry, return GR_SUCCESS; } -int free_cloudy_data(cloudy_data *my_cloudy, chemistry_data *my_chemistry, int primordial) { +int grackle::impl::free_cloudy_data(cloudy_data *my_cloudy, + chemistry_data *my_chemistry, + int primordial) { int i; for(i = 0; i < my_cloudy->grid_rank; i++) { @@ -346,3 +350,4 @@ int free_cloudy_data(cloudy_data *my_cloudy, chemistry_data *my_chemistry, int p } return GR_SUCCESS; } + diff --git a/src/clib/initialize_cloudy_data.hpp b/src/clib/initialize_cloudy_data.hpp index fc6301493..dfc73c977 100644 --- a/src/clib/initialize_cloudy_data.hpp +++ b/src/clib/initialize_cloudy_data.hpp @@ -15,9 +15,7 @@ #include "grackle.h" -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ +namespace grackle::impl { // initialize cloudy cooling data int initialize_cloudy_data(chemistry_data* my_chemistry, @@ -28,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 */ From 41c3f94d1053dc72da77efe9e7a87526b3c656be Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:02:43 -0500 Subject: [PATCH 05/38] rename and mv declaration of MAX_LINE_LENGTH_MACRO --- src/clib/grackle_macros.h | 1 - src/clib/initialize_cloudy_data.cpp | 7 +++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/clib/grackle_macros.h b/src/clib/grackle_macros.h index 6daae89f3..995965cef 100644 --- a/src/clib/grackle_macros.h +++ b/src/clib/grackle_macros.h @@ -129,7 +129,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_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 9c9d0f9a2..5427cdfec 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -21,6 +21,8 @@ #include "initialize_cloudy_data.hpp" #define SMALL_LOG_VALUE (-99.0) +#define MAX_PARAMETER_NAME_LENGTH (512) + namespace { // stuff inside of an anonymous namespace is read-only @@ -52,8 +54,9 @@ int grackle::impl::initialize_cloudy_data( double *temp_data; long long temp_int; long long *temp_int_arr; - char parameter_name[MAX_LINE_LENGTH]; - const std::size_t pname_bufsize = static_cast(MAX_LINE_LENGTH); + char parameter_name[MAX_PARAMETER_NAME_LENGTH]; + const std::size_t pname_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); From d9708c4c529aa78fe996b6dd46d06c2664e121ad Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:06:58 -0500 Subject: [PATCH 06/38] change handling of temp_int_arr --- src/clib/initialize_cloudy_data.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 5427cdfec..f5b13dda7 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -53,7 +53,6 @@ int grackle::impl::initialize_cloudy_data( double *temp_data; long long temp_int; - long long *temp_int_arr; char parameter_name[MAX_PARAMETER_NAME_LENGTH]; const std::size_t pname_bufsize = static_cast(MAX_PARAMETER_NAME_LENGTH); @@ -138,13 +137,13 @@ int grackle::impl::initialize_cloudy_data( } // Grid dimension. - temp_int_arr = (long long*)malloc(my_cloudy->grid_rank * sizeof(long long)); + long long* temp_int_arr = new long long[my_cloudy->grid_rank]; attr_id = H5Aopen_name(dset_id, "Dimension"); if (attr_id == h5_error) { std::fprintf(stderr,"Failed to open Dimension attribute in Cooling dataset.\n"); return GR_FAIL; } - status = H5Aread(attr_id, HDF5_I8,temp_int_arr); + status = H5Aread(attr_id, HDF5_I8, temp_int_arr); if (attr_id == h5_error) { std::fprintf(stderr,"Failed to read Dimension attribute in Cooling dataset.\n"); return GR_FAIL; @@ -164,7 +163,7 @@ int grackle::impl::initialize_cloudy_data( std::fprintf(stderr,"Failed to close Dimension attribute in Cooling dataset.\n"); return GR_FAIL; } - free(temp_int_arr); + delete[] temp_int_arr; // Grid parameters. for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { From 8cbbe48b96d3e353d204693d934e6bbfaa331fa3 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:13:49 -0500 Subject: [PATCH 07/38] Remove a use of malloc & free while reading in a parameter grid --- src/clib/initialize_cloudy_data.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index f5b13dda7..6855cbc03 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -175,15 +175,15 @@ int grackle::impl::initialize_cloudy_data( std::snprintf(parameter_name, pname_bufsize, "Temperature"); } - temp_data = (double*) malloc(my_cloudy->grid_dimension[q] * sizeof(double)); - attr_id = H5Aopen_name(dset_id, parameter_name); if (attr_id == h5_error) { std::fprintf(stderr,"Failed to open %s attribute in Cooling dataset.\n", parameter_name); return GR_FAIL; } - status = H5Aread(attr_id, HDF5_R8, temp_data); + + double* tmp_param_buf = new double[my_cloudy->grid_dimension[q]]; + status = H5Aread(attr_id, HDF5_R8, tmp_param_buf); if (attr_id == h5_error) { std::fprintf(stderr,"Failed to read %s attribute in Cooling dataset.\n", parameter_name); @@ -194,14 +194,15 @@ int grackle::impl::initialize_cloudy_data( sizeof(double)); for (long long w = 0LL; w < my_cloudy->grid_dimension[q]; w++) { if (q < my_cloudy->grid_rank - 1) { - my_cloudy->grid_parameters[q][w] = (double) temp_data[w]; + my_cloudy->grid_parameters[q][w] = (double) tmp_param_buf[w]; } else { // convert temeperature to log - my_cloudy->grid_parameters[q][w] = (double) std::log10(temp_data[w]); + my_cloudy->grid_parameters[q][w] = (double) std::log10(tmp_param_buf[w]); } } + delete[] tmp_param_buf; if (grackle_verbose) { std::fprintf(stdout, "%s: %g to %g (%lld steps).\n", @@ -216,8 +217,6 @@ int grackle::impl::initialize_cloudy_data( parameter_name); return GR_FAIL; } - free(temp_data); - } // Read Cooling data. From abf852bb9163be46c51162d414170f083ecbbbf5 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:25:09 -0500 Subject: [PATCH 08/38] replace another use of malloc & free with new & delete --- src/clib/initialize_cloudy_data.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 6855cbc03..2ac8baa0c 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -224,25 +224,27 @@ int grackle::impl::initialize_cloudy_data( for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { my_cloudy->data_size *= my_cloudy->grid_dimension[q]; } - temp_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); + double* tmp_cool_data = new double[my_cloudy->data_size]; - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data); + status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, + tmp_cool_data); if (grackle_verbose) std::fprintf(stdout,"Reading Cloudy Cooling dataset.\n"); if (status == h5_error) { std::fprintf(stderr,"Failed to read Cooling dataset.\n"); + delete[] tmp_cool_data; return GR_FAIL; } my_cloudy->cooling_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); for (long long q = 0LL; q < my_cloudy->data_size; q++) { - my_cloudy->cooling_data[q] = temp_data[q] > 0 ? - (double) std::log10(temp_data[q]) : (double) SMALL_LOG_VALUE; + my_cloudy->cooling_data[q] = tmp_cool_data[q] > 0 ? + (double) std::log10(tmp_cool_data[q]) : (double) SMALL_LOG_VALUE; // Convert to code units. my_cloudy->cooling_data[q] -= std::log10(CoolUnit); } - free(temp_data); + delete[] tmp_cool_data; status = H5Dclose(dset_id); if (status == h5_error) { From 2ac76b8c9f899a964646d4a8e29df7868d982f5e Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:33:41 -0500 Subject: [PATCH 09/38] replace a 3rd use of malloc & free with new & delete --- src/clib/initialize_cloudy_data.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 2ac8baa0c..1cf27839f 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -51,7 +51,6 @@ int grackle::impl::initialize_cloudy_data( code_units *my_units, int read_data) { - double *temp_data; long long temp_int; char parameter_name[MAX_PARAMETER_NAME_LENGTH]; const std::size_t pname_bufsize = @@ -255,7 +254,7 @@ int grackle::impl::initialize_cloudy_data( // Read Heating data. if (my_chemistry->UVbackground == 1) { - temp_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); + double* tmp_heat_data = new double[my_cloudy->data_size]; std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/Heating", group_name); @@ -266,7 +265,8 @@ int grackle::impl::initialize_cloudy_data( return GR_FAIL; } - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data); + status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, + tmp_heat_data); if (grackle_verbose) std::fprintf(stdout,"Reading Cloudy Heating dataset.\n"); if (status == h5_error) { @@ -276,13 +276,13 @@ int grackle::impl::initialize_cloudy_data( my_cloudy->heating_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); for (long long q = 0LL; q < my_cloudy->data_size; q++) { - my_cloudy->heating_data[q] = temp_data[q] > 0 ? - (double) std::log10(temp_data[q]) : (double) SMALL_LOG_VALUE; + my_cloudy->heating_data[q] = tmp_heat_data[q] > 0 ? + (double) std::log10(tmp_heat_data[q]) : (double) SMALL_LOG_VALUE; // Convert to code units. my_cloudy->heating_data[q] -= std::log10(CoolUnit); } - free(temp_data); + delete[] tmp_heat_data; status = H5Dclose(dset_id); if (status == h5_error) { From 15cd75daf62db253f59a180e9b9d31d028586f74 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:43:06 -0500 Subject: [PATCH 10/38] remove remaining occurrences of malloc from initialize_cloudy_data --- src/clib/initialize_cloudy_data.cpp | 34 ++++++++++++++++++----------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 1cf27839f..2e04ff1c2 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -189,8 +189,7 @@ int grackle::impl::initialize_cloudy_data( return GR_FAIL; } - my_cloudy->grid_parameters[q] = (double*) malloc(my_cloudy->grid_dimension[q] * - sizeof(double)); + my_cloudy->grid_parameters[q] = new double[my_cloudy->grid_dimension[q]]; for (long long w = 0LL; w < my_cloudy->grid_dimension[q]; w++) { if (q < my_cloudy->grid_rank - 1) { my_cloudy->grid_parameters[q][w] = (double) tmp_param_buf[w]; @@ -235,7 +234,7 @@ int grackle::impl::initialize_cloudy_data( return GR_FAIL; } - my_cloudy->cooling_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); + my_cloudy->cooling_data = new double[my_cloudy->data_size]; for (long long q = 0LL; q < my_cloudy->data_size; q++) { my_cloudy->cooling_data[q] = tmp_cool_data[q] > 0 ? (double) std::log10(tmp_cool_data[q]) : (double) SMALL_LOG_VALUE; @@ -274,7 +273,7 @@ int grackle::impl::initialize_cloudy_data( return GR_FAIL; } - my_cloudy->heating_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); + my_cloudy->heating_data = new double[my_cloudy->data_size]; for (long long q = 0LL; q < my_cloudy->data_size; q++) { my_cloudy->heating_data[q] = tmp_heat_data[q] > 0 ? (double) std::log10(tmp_heat_data[q]) : (double) SMALL_LOG_VALUE; @@ -295,7 +294,7 @@ int grackle::impl::initialize_cloudy_data( if (my_chemistry->primordial_chemistry == 0 && std::strcmp(group_name, "Primordial") == 0) { - my_cloudy->mmw_data = (double*) malloc(my_cloudy->data_size * sizeof(double)); + my_cloudy->mmw_data = new double[my_cloudy->data_size]; std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/MMW", group_name); @@ -338,18 +337,27 @@ int grackle::impl::initialize_cloudy_data( int grackle::impl::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]); + 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; + } } - GRACKLE_FREE(my_cloudy->cooling_data); - if (my_chemistry->UVbackground == 1) { - GRACKLE_FREE(my_cloudy->heating_data); + if (my_cloudy->cooling_data != nullptr) { + delete[] my_cloudy->cooling_data; + my_cloudy->cooling_data = nullptr; } - if (my_chemistry->primordial_chemistry == 0 && primordial) { - GRACKLE_FREE(my_cloudy->mmw_data); + + 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; } From 3bc86498ea996306e1705b56f3c74963d1edb80e Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 14:45:38 -0500 Subject: [PATCH 11/38] cleanup include directives of initialize_cloudy_data.cpp --- src/clib/initialize_cloudy_data.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 2e04ff1c2..013e449ef 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -11,13 +11,12 @@ / software. ************************************************************************/ -#include -#include // std::fprintf, std::snprintf, stderr, stdio #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" +#include "grackle_macros.h" // HDF5_I8, HDF5_R8, TRUE #include "initialize_cloudy_data.hpp" #define SMALL_LOG_VALUE (-99.0) From 11964ff172fdb7c13ba9a23a83461628c4c4f773 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 15:04:54 -0500 Subject: [PATCH 12/38] initial conversion of `initialize_UVbackground_data.{c->cpp}` --- .clang-format-ignore | 2 +- src/clib/CMakeLists.txt | 2 +- ...ata.c => initialize_UVbackground_data.cpp} | 36 +++++++++---------- ...ata.h => initialize_UVbackground_data.hpp} | 0 src/clib/initialize_chemistry_data.cpp | 2 +- 5 files changed, 21 insertions(+), 21 deletions(-) rename src/clib/{initialize_UVbackground_data.c => initialize_UVbackground_data.cpp} (88%) rename src/clib/{initialize_UVbackground_data.h => initialize_UVbackground_data.hpp} (100%) diff --git a/.clang-format-ignore b/.clang-format-ignore index 900c1bd4b..181cf4813 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -30,7 +30,7 @@ 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.cpp src/clib/initialize_dust_yields.cpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 68858ad40..d80784907 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -90,7 +90,6 @@ add_library(Grackle_Grackle dynamic_api.c grackle_units.c index_helper.c - initialize_UVbackground_data.c initialize_UVbackground_data.h rate_functions.c set_default_chemistry_parameters.c solve_chemistry.c @@ -116,6 +115,7 @@ add_library(Grackle_Grackle 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 diff --git a/src/clib/initialize_UVbackground_data.c b/src/clib/initialize_UVbackground_data.cpp similarity index 88% rename from src/clib/initialize_UVbackground_data.c rename to src/clib/initialize_UVbackground_data.cpp index 8a5bf9250..e0b8632ad 100644 --- a/src/clib/initialize_UVbackground_data.c +++ b/src/clib/initialize_UVbackground_data.cpp @@ -17,10 +17,10 @@ #include "grackle.h" #include "grackle_macros.h" -#include "initialize_UVbackground_data.h" +#include "initialize_UVbackground_data.hpp" // function prototypes -int read_dataset(hid_t file_id, char *dset_name, double *buffer); +int read_dataset(hid_t file_id, const char *dset_name, double *buffer); /** * Initializes an empty #UVBtable struct with zeros and NULLs. @@ -132,27 +132,27 @@ int initialize_UVbackground_data(chemistry_data *my_chemistry, // 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)); + my_rates->UVbackground_table.z = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k24 = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k25 = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k26 = (double*)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.k27 = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k28 = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k29 = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k30 = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.k31 = (double*)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)); + my_rates->UVbackground_table.piHI = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.piHeII = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.piHeI = (double*)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)); + my_rates->UVbackground_table.crsHI = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.crsHeII = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.crsHeI = (double*)malloc(Nz * sizeof(double)); } @@ -341,7 +341,7 @@ int initialize_UVbackground_data(chemistry_data *my_chemistry, -int read_dataset(hid_t file_id, char *dset_name, double *buffer) { +int read_dataset(hid_t file_id, const char *dset_name, double *buffer) { hid_t dset_id; herr_t status; herr_t h5_error = -1; diff --git a/src/clib/initialize_UVbackground_data.h b/src/clib/initialize_UVbackground_data.hpp similarity index 100% rename from src/clib/initialize_UVbackground_data.h rename to src/clib/initialize_UVbackground_data.hpp diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 37fa313f9..4790fd8f6 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -24,7 +24,7 @@ #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" From c1319739698aed48b552303d8b61055993d64146 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 15:23:48 -0500 Subject: [PATCH 13/38] mv UVBtable methods into grackle::impl and factor out cleanup logic --- src/clib/initialize_UVbackground_data.cpp | 25 +++++++++++++++---- src/clib/initialize_UVbackground_data.hpp | 13 +++++----- src/clib/initialize_chemistry_data.cpp | 29 ++++------------------- 3 files changed, 31 insertions(+), 36 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index e0b8632ad..c2477fe57 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -25,7 +25,7 @@ int read_dataset(hid_t file_id, const char *dset_name, double *buffer); /** * Initializes an empty #UVBtable struct with zeros and NULLs. */ -void initialize_empty_UVBtable_struct(UVBtable *table) +void grackle::impl::initialize_empty_UVBtable_struct(UVBtable *table) { table->Nz = 0LL; table->z = NULL; @@ -45,10 +45,9 @@ void initialize_empty_UVBtable_struct(UVBtable *table) table->crsHeI = NULL; } - // Initialize UV Background data -int initialize_UVbackground_data(chemistry_data *my_chemistry, - chemistry_data_storage *my_rates) +int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, + chemistry_data_storage *my_rates) { long long Nz; @@ -339,6 +338,24 @@ int initialize_UVbackground_data(chemistry_data *my_chemistry, return SUCCESS; } +void grackle::impl::free_UVBtable(UVBtable *table) +{ + GRACKLE_FREE(table->z); + GRACKLE_FREE(table->k24); + GRACKLE_FREE(table->k25); + GRACKLE_FREE(table->k26); + GRACKLE_FREE(table->k27); + GRACKLE_FREE(table->k28); + GRACKLE_FREE(table->k29); + GRACKLE_FREE(table->k30); + GRACKLE_FREE(table->k31); + GRACKLE_FREE(table->piHI); + GRACKLE_FREE(table->piHeII); + GRACKLE_FREE(table->piHeI); + GRACKLE_FREE(table->crsHI); + GRACKLE_FREE(table->crsHeII); + GRACKLE_FREE(table->crsHeI); +} int read_dataset(hid_t file_id, const char *dset_name, double *buffer) { diff --git a/src/clib/initialize_UVbackground_data.hpp b/src/clib/initialize_UVbackground_data.hpp index 5c4580a7b..f2028291c 100644 --- a/src/clib/initialize_UVbackground_data.hpp +++ b/src/clib/initialize_UVbackground_data.hpp @@ -15,19 +15,18 @@ #include "grackle.h" -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ +namespace grackle::impl { -/// Initializes an empty UVBtable struct with zeros and NULLs. +/// Initializes an empty UVBtable struct with zeros and nullptrs. void initialize_empty_UVBtable_struct(UVBtable* table); /// 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 4790fd8f6..27b6e4cb5 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -446,8 +446,9 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry, } /* Initialize UV Background data. */ - initialize_empty_UVBtable_struct(&(my_rates->UVbackground_table)); - if (initialize_UVbackground_data(my_chemistry, my_rates) == GR_FAIL) { + grackle::impl::initialize_empty_UVBtable_struct(&my_rates->UVbackground_table); + if (grackle::impl::initialize_UVbackground_data(my_chemistry, my_rates) + != GR_SUCCESS) { fprintf(stderr, "Error in initialize_UVbackground_data.\n"); return GR_FAIL; } @@ -597,29 +598,7 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry, /* primordial = */ 1); grackle::impl::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_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"); From ad8388d044bef66b2505b2cce152e6b01bcde26d Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 15:30:52 -0500 Subject: [PATCH 14/38] initialize_UVbackground_data: replace malloc with new --- src/clib/initialize_UVbackground_data.cpp | 100 ++++++++++++---------- 1 file changed, 53 insertions(+), 47 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index c2477fe57..285628c2f 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -11,7 +11,6 @@ / software. ************************************************************************/ -#include #include #include "hdf5.h" #include "grackle.h" @@ -28,21 +27,21 @@ int read_dataset(hid_t file_id, const char *dset_name, double *buffer); void grackle::impl::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; + 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; } // Initialize UV Background data @@ -127,31 +126,30 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, H5Dclose(dset_id); - // Now allocate memory for UV background table. my_rates->UVbackground_table.Nz = Nz; - my_rates->UVbackground_table.z = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k24 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k25 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k26 = (double*)malloc(Nz * sizeof(double)); + 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 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k28 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k29 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k30 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.k31 = (double*)malloc(Nz * sizeof(double)); + 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 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.piHeII = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.piHeI = (double*)malloc(Nz * sizeof(double)); + 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 = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.crsHeII = (double*)malloc(Nz * sizeof(double)); - my_rates->UVbackground_table.crsHeI = (double*)malloc(Nz * sizeof(double)); + my_rates->UVbackground_table.crsHI = new double[Nz]; + my_rates->UVbackground_table.crsHeII = new double[Nz]; + my_rates->UVbackground_table.crsHeI = new double[Nz]; } @@ -340,21 +338,29 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, void grackle::impl::free_UVBtable(UVBtable *table) { - GRACKLE_FREE(table->z); - GRACKLE_FREE(table->k24); - GRACKLE_FREE(table->k25); - GRACKLE_FREE(table->k26); - GRACKLE_FREE(table->k27); - GRACKLE_FREE(table->k28); - GRACKLE_FREE(table->k29); - GRACKLE_FREE(table->k30); - GRACKLE_FREE(table->k31); - GRACKLE_FREE(table->piHI); - GRACKLE_FREE(table->piHeII); - GRACKLE_FREE(table->piHeI); - GRACKLE_FREE(table->crsHI); - GRACKLE_FREE(table->crsHeII); - GRACKLE_FREE(table->crsHeI); + // 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); } From dfd1a8779923c6c9ae45e74fa7c72a5b0c5f4b49 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 15:35:20 -0500 Subject: [PATCH 15/38] initialize_UVbackground_data: more cppification --- src/clib/initialize_UVbackground_data.cpp | 120 +++++++++++----------- 1 file changed, 60 insertions(+), 60 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index 285628c2f..935c1275c 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -11,7 +11,7 @@ / software. ************************************************************************/ -#include +#include #include "hdf5.h" #include "grackle.h" #include "grackle_macros.h" @@ -53,11 +53,11 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // Return if no UV background selected or using fully tabulated cooling. if (my_chemistry->UVbackground == 0 || my_chemistry->primordial_chemistry == 0) - return SUCCESS; + return GR_SUCCESS; if (grackle_verbose) - fprintf(stdout, "Initializing UV background.\n"); + std::fprintf(stdout, "Initializing UV background.\n"); // Read in UV background data from hdf5 file. @@ -67,7 +67,7 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, herr_t h5_error = -1; if (grackle_verbose) - fprintf(stdout, "Reading UV background data from %s.\n", + std::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); @@ -77,9 +77,9 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, dset_id = H5Dopen(file_id, "/UVBRates/Info"); if (dset_id == h5_error) { - fprintf(stderr, "Can't open 'Info' dataset in %s.\n", + std::fprintf(stderr, "Can't open 'Info' dataset in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } int strlen = (int)(H5Dget_storage_size(dset_id)); @@ -90,8 +90,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, 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; + std::fprintf(stderr, "Failed to read dataset 'Info'.\n"); + return GR_FAIL; } H5Tclose(memtype); @@ -103,23 +103,23 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, dset_id = H5Dopen(file_id, "/UVBRates/z"); if (dset_id == h5_error) { - fprintf(stderr, "Can't open redshift dataset ('z') in %s.\n", + std::fprintf(stderr, "Can't open redshift dataset ('z') in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } dspace_id = H5Dget_space(dset_id); if (dspace_id == h5_error) { - fprintf(stderr, "Error opening dataspace for dataset 'z' in %s.\n", + std::fprintf(stderr, "Error opening dataspace for dataset 'z' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } Nz = H5Sget_simple_extent_npoints(dspace_id); if(Nz <= 0) { - fprintf(stderr, "Redshift dataset ('z') has inappropriate size = %lld in %s.\n", + std::fprintf(stderr, "Redshift dataset ('z') has inappropriate size = %lld in %s.\n", Nz, my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } H5Sclose(dspace_id); @@ -159,33 +159,33 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** Redshift *** if(! read_dataset(file_id, "/UVBRates/z", my_rates->UVbackground_table.z) ) { - fprintf(stderr, "Error reading dataset 'z' in %s.\n", + std::fprintf(stderr, "Error reading dataset 'z' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k24' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k25' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k26' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } if (my_chemistry->primordial_chemistry > 1) { @@ -193,41 +193,41 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** 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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k27' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k28' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k29' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k30' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Chemistry/k31' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } } @@ -235,25 +235,25 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** 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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHI' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHeII' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_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", + std::fprintf(stderr, "Error reading dataset '/UVBRates/Photoheating/piHeI' in %s.\n", my_chemistry->grackle_data_file); - return FAIL; + return GR_FAIL; } @@ -261,28 +261,28 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** 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", + std::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; + 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) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/heii_avg_crs' in %s.\n", + std::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; + 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) ) { - fprintf(stderr, "Error reading dataset '/UVBRates/CrossSections/hei_avg_crs' in %s.\n", + std::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; + std::fprintf(stderr, "In order to use self-shielding, you must use the shielding datasets\n"); + return GR_FAIL; } } @@ -296,9 +296,9 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // 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", + std::fprintf(stdout, "UV background information:\n"); + std::fprintf(stdout, " %s\n",info_string); + std::fprintf(stdout, " z_min = %6.3f\n z_max = %6.3f\n", my_rates->UVbackground_table.zmin, my_rates->UVbackground_table.zmax); } @@ -308,32 +308,32 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, my_chemistry->UVbackground_redshift_on = my_rates->UVbackground_table.zmax; if (grackle_verbose) - fprintf(stdout, "Setting UVbackground_redshift_on to %f.\n", + 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) - fprintf(stdout, "Setting UVbackground_redshift_fullon to %f.\n", + 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) - fprintf(stdout, "Setting UVbackground_redshift_drop to %f.\n", + 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) - fprintf(stdout, "Setting UVbackground_redshift_off to %f.\n", + std::fprintf(stdout, "Setting UVbackground_redshift_off to %f.\n", my_chemistry->UVbackground_redshift_off); } - return SUCCESS; + return GR_SUCCESS; } void grackle::impl::free_UVBtable(UVBtable *table) @@ -371,17 +371,17 @@ int read_dataset(hid_t file_id, const char *dset_name, double *buffer) { dset_id = H5Dopen(file_id, dset_name); if (dset_id == h5_error) { - fprintf(stderr, "Failed to open dataset 'z'.\n"); - return FAIL; + std::fprintf(stderr, "Failed to open dataset 'z'.\n"); + return GR_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; + std::fprintf(stderr, "Failed to read dataset 'z'.\n"); + return GR_FAIL; } H5Dclose(dset_id); - return SUCCESS; + return GR_SUCCESS; } From 597f9cf64cb39bc2bcb7409d70c611d41ba18c9a Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 15:37:52 -0500 Subject: [PATCH 16/38] initialize_UVbackground_data: use std::vector instead of a C++ compiler extension --- src/clib/initialize_UVbackground_data.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index 935c1275c..f399acb6f 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -12,6 +12,7 @@ ************************************************************************/ #include +#include #include "hdf5.h" #include "grackle.h" #include "grackle_macros.h" @@ -83,12 +84,13 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, } int strlen = (int)(H5Dget_storage_size(dset_id)); - char info_string[strlen+1]; + std::vector 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); + status = H5Dread(dset_id, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, + info_string.data()); if (status == h5_error) { std::fprintf(stderr, "Failed to read dataset 'Info'.\n"); return GR_FAIL; @@ -297,7 +299,7 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // 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); + 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); From 303a3adc3a801e33014cbf0469ff9932a420e994 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 15:42:15 -0500 Subject: [PATCH 17/38] call initialize_empty_UVBtable_struct within initialize_UVbackground_data --- src/clib/initialize_UVbackground_data.cpp | 11 +++++++---- src/clib/initialize_UVbackground_data.hpp | 3 --- src/clib/initialize_chemistry_data.cpp | 1 - 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index f399acb6f..8f6bddfc0 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -22,10 +22,10 @@ // function prototypes int read_dataset(hid_t file_id, const char *dset_name, double *buffer); -/** - * Initializes an empty #UVBtable struct with zeros and NULLs. - */ -void grackle::impl::initialize_empty_UVBtable_struct(UVBtable *table) +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; @@ -45,10 +45,13 @@ void grackle::impl::initialize_empty_UVBtable_struct(UVBtable *table) 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); long long Nz; // Return if no UV background selected or using fully tabulated cooling. diff --git a/src/clib/initialize_UVbackground_data.hpp b/src/clib/initialize_UVbackground_data.hpp index f2028291c..f24587e2e 100644 --- a/src/clib/initialize_UVbackground_data.hpp +++ b/src/clib/initialize_UVbackground_data.hpp @@ -17,9 +17,6 @@ namespace grackle::impl { -/// Initializes an empty UVBtable struct with zeros and nullptrs. -void initialize_empty_UVBtable_struct(UVBtable* table); - /// Initialize UV Background data int initialize_UVbackground_data(chemistry_data* my_chemistry, chemistry_data_storage* my_rates); diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 27b6e4cb5..10dd5118b 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -446,7 +446,6 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry, } /* Initialize UV Background data. */ - grackle::impl::initialize_empty_UVBtable_struct(&my_rates->UVbackground_table); if (grackle::impl::initialize_UVbackground_data(my_chemistry, my_rates) != GR_SUCCESS) { fprintf(stderr, "Error in initialize_UVbackground_data.\n"); From e3e915a4f8ba2f4682962c2e52035567aed74ed4 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 16:45:45 -0500 Subject: [PATCH 18/38] introduce h5io utility file --- src/clib/CMakeLists.txt | 1 + src/clib/initialize_UVbackground_data.cpp | 29 ++------------- src/clib/utils/h5io.cpp | 43 +++++++++++++++++++++++ src/clib/utils/h5io.hpp | 25 +++++++++++++ 4 files changed, 72 insertions(+), 26 deletions(-) create mode 100644 src/clib/utils/h5io.cpp create mode 100644 src/clib/utils/h5io.hpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index d80784907..7bcb5ecaa 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -125,6 +125,7 @@ add_library(Grackle_Grackle solve_rate_cool_g-cpp.cpp solve_rate_cool_g-cpp.h step_rate_newton_raphson.hpp time_deriv_0d.hpp + utils/h5io.cpp utils/h5io.hpp utils-cpp.cpp utils-cpp.hpp utils-field.hpp fortran_func_wrappers.hpp diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index 8f6bddfc0..62d28ef38 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -16,12 +16,10 @@ #include "hdf5.h" #include "grackle.h" #include "grackle_macros.h" +#include "utils/h5io.hpp" #include "initialize_UVbackground_data.hpp" -// function prototypes -int read_dataset(hid_t file_id, const char *dset_name, double *buffer); - namespace { // stuff inside an anonymous namespace is only used in this file /// Initializes an empty #UVBtable struct with zeros and nullptrs. @@ -160,6 +158,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // Now read everything. + using ::grackle::impl::h5io::read_dataset; + // *** Redshift *** if(! read_dataset(file_id, "/UVBRates/z", @@ -367,26 +367,3 @@ void grackle::impl::free_UVBtable(UVBtable *table) cleanup_fn(table->crsHeII); cleanup_fn(table->crsHeI); } - - -int read_dataset(hid_t file_id, const 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) { - std::fprintf(stderr, "Failed to open dataset 'z'.\n"); - return GR_FAIL; - } - - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); - if (status == h5_error) { - std::fprintf(stderr, "Failed to read dataset 'z'.\n"); - return GR_FAIL; - } - - H5Dclose(dset_id); - - return GR_SUCCESS; -} diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp new file mode 100644 index 000000000..4c7a7ed38 --- /dev/null +++ b/src/clib/utils/h5io.cpp @@ -0,0 +1,43 @@ +//===----------------------------------------------------------------------===// +// +// 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 +/// +//===----------------------------------------------------------------------===// +#include + +#include "hdf5.h" +#include "grackle.h" +#include "grackle_macros.h" + +#include "h5io.hpp" +#include "../status_reporting.h" + +/// 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) { + hid_t dset_id; + herr_t status; + herr_t h5_error = -1; + + dset_id = H5Dopen(file_id, dset_name); + if (dset_id == h5_error) { + std::fprintf(stderr, "Failed to open dataset 'z'.\n"); + return GR_FAIL; + } + + status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); + if (status == h5_error) { + std::fprintf(stderr, "Failed to read dataset 'z'.\n"); + return GR_FAIL; + } + + H5Dclose(dset_id); + + return GR_SUCCESS; +} diff --git a/src/clib/utils/h5io.hpp b/src/clib/utils/h5io.hpp new file mode 100644 index 000000000..3e73d793c --- /dev/null +++ b/src/clib/utils/h5io.hpp @@ -0,0 +1,25 @@ +//===----------------------------------------------------------------------===// +// +// 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 "hdf5.h" + +namespace grackle::impl::h5io { + +/// read the dataset named dset_name from file_id into buffer +int read_dataset(hid_t file_id, const char* dset_name, double* buffer); + +} // namespace grackle::impl::h5io + +#endif /* UTILS_H5IO_HPP */ From bcd3074891c71da1aa96692e17f621e03ad627a4 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 17:09:22 -0500 Subject: [PATCH 19/38] start using read_dataset in initialize_cloudy_data --- src/clib/initialize_cloudy_data.cpp | 46 ++++++++++++++++------------- src/clib/utils/h5io.cpp | 19 +++++++----- 2 files changed, 36 insertions(+), 29 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 013e449ef..309793db2 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -18,6 +18,7 @@ #include "grackle.h" #include "grackle_macros.h" // HDF5_I8, HDF5_R8, TRUE #include "initialize_cloudy_data.hpp" +#include "utils/h5io.hpp" #define SMALL_LOG_VALUE (-99.0) #define MAX_PARAMETER_NAME_LENGTH (512) @@ -51,8 +52,9 @@ int grackle::impl::initialize_cloudy_data( { long long temp_int; + char dset_name[MAX_PARAMETER_NAME_LENGTH]; char parameter_name[MAX_PARAMETER_NAME_LENGTH]; - const std::size_t pname_bufsize = + 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. @@ -105,12 +107,12 @@ int grackle::impl::initialize_cloudy_data( // Open cooling dataset and get grid dimensions. - std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/Cooling", + std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Cooling", group_name); - dset_id = H5Dopen(file_id, parameter_name); + dset_id = H5Dopen(file_id, dset_name); if (dset_id == h5_error) { - std::fprintf(stderr,"Can't open Cooling in %s.\n", - my_chemistry->grackle_data_file); + std::fprintf(stderr,"Can't open \"%s\" dataset in %s.\n", + dset_name, my_chemistry->grackle_data_file); return GR_FAIL; } @@ -167,10 +169,10 @@ int grackle::impl::initialize_cloudy_data( for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { if (q < my_cloudy->grid_rank - 1) { - std::snprintf(parameter_name, pname_bufsize, "Parameter%lld",(q+1)); + std::snprintf(parameter_name, name_bufsize, "Parameter%lld",(q+1)); } else { - std::snprintf(parameter_name, pname_bufsize, "Temperature"); + std::snprintf(parameter_name, name_bufsize, "Temperature"); } attr_id = H5Aopen_name(dset_id, parameter_name); @@ -216,20 +218,27 @@ int grackle::impl::initialize_cloudy_data( } } + // to make the logic more concise, we do something a *little* silly + // -> we effectively close the dataset and then call a function that reopens + // and closes the dataset + status = H5Dclose(dset_id); + if (status == h5_error) { + std::fprintf(stderr,"Error while to closing %s dataset.\n", dset_name); + return GR_FAIL; + } + // Read Cooling data. my_cloudy->data_size = 1; for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { my_cloudy->data_size *= my_cloudy->grid_dimension[q]; } double* tmp_cool_data = new double[my_cloudy->data_size]; - - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, - tmp_cool_data); - if (grackle_verbose) - std::fprintf(stdout,"Reading Cloudy Cooling dataset.\n"); - if (status == h5_error) { - std::fprintf(stderr,"Failed to read Cooling dataset.\n"); + if (grackle_verbose) { + std::fprintf(stdout,"Reading from \"%s\" dataset.\n", dset_name); + } + if (h5io::read_dataset(file_id, dset_name, tmp_cool_data) != GR_SUCCESS) { delete[] tmp_cool_data; + std::fprintf(stderr,"Failed to read Cooling dataset.\n"); return GR_FAIL; } @@ -243,18 +252,13 @@ int grackle::impl::initialize_cloudy_data( } delete[] tmp_cool_data; - status = H5Dclose(dset_id); - if (status == h5_error) { - std::fprintf(stderr,"Failed to close Cooling dataset.\n"); - return GR_FAIL; - } // Read Heating data. if (my_chemistry->UVbackground == 1) { double* tmp_heat_data = new double[my_cloudy->data_size]; - std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/Heating", + std::snprintf(parameter_name, name_bufsize, "/CoolingRates/%s/Heating", group_name); dset_id = H5Dopen(file_id, parameter_name); if (dset_id == h5_error) { @@ -295,7 +299,7 @@ int grackle::impl::initialize_cloudy_data( my_cloudy->mmw_data = new double[my_cloudy->data_size]; - std::snprintf(parameter_name, pname_bufsize, "/CoolingRates/%s/MMW", + std::snprintf(parameter_name, name_bufsize, "/CoolingRates/%s/MMW", group_name); dset_id = H5Dopen(file_id, parameter_name); if (dset_id == h5_error) { diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 4c7a7ed38..a1672fee3 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -14,26 +14,29 @@ #include "hdf5.h" #include "grackle.h" #include "grackle_macros.h" - #include "h5io.hpp" #include "../status_reporting.h" /// 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) { - hid_t dset_id; - herr_t status; - herr_t h5_error = -1; + if (dset_name == nullptr) { + std::fprintf(stderr, "dset_name is a nullptr"); + return GR_FAIL; + } + + const herr_t h5_error = -1; - dset_id = H5Dopen(file_id, dset_name); + hid_t dset_id = H5Dopen(file_id, dset_name); if (dset_id == h5_error) { - std::fprintf(stderr, "Failed to open dataset 'z'.\n"); + std::fprintf(stderr, "Failed to open dataset \"%s\".\n", dset_name); return GR_FAIL; } - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); + herr_t status = + H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer); if (status == h5_error) { - std::fprintf(stderr, "Failed to read dataset 'z'.\n"); + std::fprintf(stderr, "Failed to read dataset \"%s\".\n", dset_name); return GR_FAIL; } From 6854adfe0008c086b9e83ee7594adfd951ea7e28 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 17:18:42 -0500 Subject: [PATCH 20/38] another use of read_dataset in initialize_cloudy_data --- src/clib/initialize_cloudy_data.cpp | 25 ++++--------------------- 1 file changed, 4 insertions(+), 21 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 309793db2..812dff555 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -238,7 +238,6 @@ int grackle::impl::initialize_cloudy_data( } if (h5io::read_dataset(file_id, dset_name, tmp_cool_data) != GR_SUCCESS) { delete[] tmp_cool_data; - std::fprintf(stderr,"Failed to read Cooling dataset.\n"); return GR_FAIL; } @@ -256,23 +255,13 @@ int grackle::impl::initialize_cloudy_data( // Read Heating data. if (my_chemistry->UVbackground == 1) { - double* tmp_heat_data = new double[my_cloudy->data_size]; - std::snprintf(parameter_name, name_bufsize, "/CoolingRates/%s/Heating", + std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Heating", group_name); - dset_id = H5Dopen(file_id, parameter_name); - if (dset_id == h5_error) { - std::fprintf(stderr,"Can't open Heating in %s.\n", - my_chemistry->grackle_data_file); - return GR_FAIL; - } - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, - tmp_heat_data); - if (grackle_verbose) - std::fprintf(stdout,"Reading Cloudy Heating dataset.\n"); - if (status == h5_error) { - std::fprintf(stderr,"Failed to read Heating dataset.\n"); + double* tmp_heat_data = new double[my_cloudy->data_size]; + if (h5io::read_dataset(file_id, dset_name, tmp_heat_data) != GR_SUCCESS) { + delete[] tmp_heat_data; return GR_FAIL; } @@ -285,12 +274,6 @@ int grackle::impl::initialize_cloudy_data( my_cloudy->heating_data[q] -= std::log10(CoolUnit); } delete[] tmp_heat_data; - - status = H5Dclose(dset_id); - if (status == h5_error) { - std::fprintf(stderr,"Failed to close Heating dataset.\n"); - return GR_FAIL; - } } // Read MMW data. From 76e5f3a91a151cbf7fd226140e519dbb79e5d83a Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 17:29:48 -0500 Subject: [PATCH 21/38] the last commit in the sequence --- src/clib/initialize_cloudy_data.cpp | 27 ++++----------------------- 1 file changed, 4 insertions(+), 23 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 812dff555..686093349 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -53,7 +53,6 @@ int grackle::impl::initialize_cloudy_data( long long temp_int; char dset_name[MAX_PARAMETER_NAME_LENGTH]; - char parameter_name[MAX_PARAMETER_NAME_LENGTH]; const std::size_t name_bufsize = static_cast(MAX_PARAMETER_NAME_LENGTH); @@ -167,6 +166,7 @@ int grackle::impl::initialize_cloudy_data( // Grid parameters. for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { + char parameter_name[MAX_PARAMETER_NAME_LENGTH]; if (q < my_cloudy->grid_rank - 1) { std::snprintf(parameter_name, name_bufsize, "Parameter%lld",(q+1)); @@ -255,7 +255,6 @@ int grackle::impl::initialize_cloudy_data( // Read Heating data. if (my_chemistry->UVbackground == 1) { - std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Heating", group_name); @@ -281,30 +280,12 @@ int grackle::impl::initialize_cloudy_data( std::strcmp(group_name, "Primordial") == 0) { my_cloudy->mmw_data = new double[my_cloudy->data_size]; - - std::snprintf(parameter_name, name_bufsize, "/CoolingRates/%s/MMW", - group_name); - dset_id = H5Dopen(file_id, parameter_name); - if (dset_id == h5_error) { - std::fprintf(stderr,"Can't open MMW in %s.\n", - my_chemistry->grackle_data_file); + if (h5io::read_dataset(file_id, "/CoolingRates/Primordial/MMW", + my_cloudy->mmw_data) != GR_SUCCESS) { + // nothing to cleanup right now return GR_FAIL; } - status = H5Dread(dset_id, HDF5_R8, H5S_ALL, H5S_ALL, H5P_DEFAULT, - my_cloudy->mmw_data); - if (grackle_verbose) - std::fprintf(stdout,"Reading Cloudy MMW dataset.\n"); - if (status == h5_error) { - std::fprintf(stderr,"Failed to read MMW dataset.\n"); - return GR_FAIL; - } - - status = H5Dclose(dset_id); - if (status == h5_error) { - std::fprintf(stderr,"Failed to close MMW dataset.\n"); - return GR_FAIL; - } } status = H5Fclose (file_id); From 70fe3331c459d04f43ea35321cec0c75b21d5640 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 24 Nov 2025 19:40:28 -0500 Subject: [PATCH 22/38] fix the failing build --- src/clib/Make.config.objects | 1 + 1 file changed, 1 insertion(+) diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 2fa8a89ee..858d7d1ee 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 \ + utils/h5io.lo \ utils.lo \ utils-cpp.lo \ internal_types.lo From 9498310a1bb19920b37ebb208950c993902aa174 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 25 Nov 2025 09:43:29 -0500 Subject: [PATCH 23/38] introduce abstractions to load the shape of hdf5 datasets --- src/clib/initialize_UVbackground_data.cpp | 34 ++----- src/clib/utils/h5io.cpp | 108 ++++++++++++++++++++-- src/clib/utils/h5io.hpp | 46 ++++++++- 3 files changed, 155 insertions(+), 33 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index 62d28ef38..e0e128fff 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -50,7 +50,6 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, chemistry_data_storage *my_rates) { initialize_empty_UVBtable_struct(&my_rates->UVbackground_table); - long long Nz; // Return if no UV background selected or using fully tabulated cooling. if (my_chemistry->UVbackground == 0 || @@ -100,34 +99,21 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, 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) { - std::fprintf(stderr, "Can't open redshift dataset ('z') in %s.\n", - my_chemistry->grackle_data_file); + 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; } - dspace_id = H5Dget_space(dset_id); - if (dspace_id == h5_error) { - std::fprintf(stderr, "Error opening dataspace for dataset 'z' in %s.\n", - my_chemistry->grackle_data_file); - return GR_FAIL; - } - - Nz = H5Sget_simple_extent_npoints(dspace_id); - if(Nz <= 0) { - std::fprintf(stderr, "Redshift dataset ('z') has inappropriate size = %lld in %s.\n", - Nz, my_chemistry->grackle_data_file); - return GR_FAIL; - } - - H5Sclose(dspace_id); - H5Dclose(dset_id); - + long long Nz = static_cast(common_shape.shape[0]); // Now allocate memory for UV background table. my_rates->UVbackground_table.Nz = Nz; diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index a1672fee3..5c681191d 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -17,30 +17,122 @@ #include "h5io.hpp" #include "../status_reporting.h" +namespace { // stuff inside an anonymous namespace is local to this file + +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; + } + } +} + +} // 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 shape_from_space(H5I_INVALID_HID); + } + + 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 shape_from_space(H5I_INVALID_HID); + } + + 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) { +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; } - const herr_t h5_error = -1; - hid_t dset_id = H5Dopen(file_id, dset_name); - if (dset_id == h5_error) { + 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_same(*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); - if (status == h5_error) { + H5Dclose(dset_id); + if (status < 0) { std::fprintf(stderr, "Failed to read dataset \"%s\".\n", dset_name); return GR_FAIL; } - H5Dclose(dset_id); - return GR_SUCCESS; } diff --git a/src/clib/utils/h5io.hpp b/src/clib/utils/h5io.hpp index 3e73d793c..767884818 100644 --- a/src/clib/utils/h5io.hpp +++ b/src/clib/utils/h5io.hpp @@ -13,12 +13,56 @@ #ifndef UTILS_H5IO_HPP #define UTILS_H5IO_HPP +#include + #include "hdf5.h" +#include "grackle.h" namespace grackle::impl::h5io { +/// 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 is valid +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; } + +/// checks whether shape_a and shape_b are the same +/// +/// returns false if the either shape is invalid +inline bool ArrayShape_is_same(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 -int read_dataset(hid_t file_id, const char* dset_name, double* 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); } // namespace grackle::impl::h5io From 0ed1ef3446e7cce33b5bd99408868e5f21e8defe Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 25 Nov 2025 10:04:53 -0500 Subject: [PATCH 24/38] more explicitly verify shapes when loading UVbackground data --- src/clib/initialize_UVbackground_data.cpp | 46 +++++++++++++++-------- 1 file changed, 30 insertions(+), 16 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index e0e128fff..b34fb6aa2 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -146,10 +146,10 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, using ::grackle::impl::h5io::read_dataset; - // *** Redshift *** if(! read_dataset(file_id, "/UVBRates/z", - my_rates->UVbackground_table.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; @@ -157,7 +157,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k24 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k24", - my_rates->UVbackground_table.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; @@ -165,7 +166,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k25 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k25", - my_rates->UVbackground_table.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; @@ -173,7 +175,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k26 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k26", - my_rates->UVbackground_table.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; @@ -183,7 +186,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k27 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k27", - my_rates->UVbackground_table.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; @@ -191,7 +195,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k28 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k28", - my_rates->UVbackground_table.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; @@ -199,7 +204,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k29 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k29", - my_rates->UVbackground_table.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; @@ -207,7 +213,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k30 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k30", - my_rates->UVbackground_table.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; @@ -215,7 +222,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** k31 *** if(! read_dataset(file_id, "/UVBRates/Chemistry/k31", - my_rates->UVbackground_table.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; @@ -225,7 +233,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** piHI *** if(! read_dataset(file_id, "/UVBRates/Photoheating/piHI", - my_rates->UVbackground_table.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; @@ -233,7 +242,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** piHeII *** if(! read_dataset(file_id, "/UVBRates/Photoheating/piHeII", - my_rates->UVbackground_table.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; @@ -241,7 +251,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** piHeI *** if(! read_dataset(file_id, "/UVBRates/Photoheating/piHeI", - my_rates->UVbackground_table.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; @@ -251,7 +262,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, if (my_chemistry->self_shielding_method > 0) { // *** crsHI *** if(! read_dataset(file_id, "/UVBRates/CrossSections/hi_avg_crs", - my_rates->UVbackground_table.crsHI) ) { + 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"); @@ -260,7 +272,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** crsHeII *** if(! read_dataset(file_id, "/UVBRates/CrossSections/heii_avg_crs", - my_rates->UVbackground_table.crsHeII) ) { + 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"); @@ -269,7 +282,8 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // *** crsHeI *** if(! read_dataset(file_id, "/UVBRates/CrossSections/hei_avg_crs", - my_rates->UVbackground_table.crsHeI) ) { + 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"); From 619da6c3dfd1efbfdf761f91956363a82c358af6 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 26 Nov 2025 12:33:54 -0500 Subject: [PATCH 25/38] Overhauled parsing of the table-axes for interpolation --- src/clib/initialize_cloudy_data.cpp | 172 ++++-------- src/clib/utils/h5io.cpp | 418 +++++++++++++++++++++++++++- src/clib/utils/h5io.hpp | 68 +++++ 3 files changed, 538 insertions(+), 120 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 686093349..1c6e4a53c 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -104,139 +104,67 @@ int grackle::impl::initialize_cloudy_data( std::fprintf(stdout, "Loading old-style Cloudy tables.\n"); } - // Open cooling dataset and get grid dimensions. - + // Determine the name of the dataset std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Cooling", group_name); - dset_id = H5Dopen(file_id, dset_name); - if (dset_id == h5_error) { - std::fprintf(stderr,"Can't open \"%s\" dataset in %s.\n", - dset_name, my_chemistry->grackle_data_file); - return GR_FAIL; - } - - // Grid rank. - attr_id = H5Aopen_name(dset_id, "Rank"); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to open Rank attribute in Cooling dataset.\n"); - return GR_FAIL; - } - status = H5Aread(attr_id, HDF5_I8, &temp_int); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to read Rank attribute in Cooling dataset.\n"); - return GR_FAIL; - } - my_cloudy->grid_rank = (long long) temp_int; - if (grackle_verbose) - std::fprintf(stdout,"Cloudy cooling grid rank: %lld.\n", my_cloudy->grid_rank); - status = H5Aclose(attr_id); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to close Rank attribute in Cooling dataset.\n"); - return GR_FAIL; - } - // Grid dimension. - long long* temp_int_arr = new long long[my_cloudy->grid_rank]; - attr_id = H5Aopen_name(dset_id, "Dimension"); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to open Dimension attribute in Cooling dataset.\n"); + // 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; } - status = H5Aread(attr_id, HDF5_I8, temp_int_arr); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to read Dimension attribute in Cooling dataset.\n"); - return GR_FAIL; - } - if (grackle_verbose) { - std::fprintf(stdout, "Cloudy cooling grid dimensions:"); - } - for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { - my_cloudy->grid_dimension[q] = (long long) temp_int_arr[q]; - if (grackle_verbose) - std::fprintf(stdout," %lld", my_cloudy->grid_dimension[q]); - } - if (grackle_verbose) - std::fprintf(stdout,".\n"); - status = H5Aclose(attr_id); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to close Dimension attribute in Cooling dataset.\n"); - return GR_FAIL; - } - delete[] temp_int_arr; - // Grid parameters. - for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { - char parameter_name[MAX_PARAMETER_NAME_LENGTH]; - - if (q < my_cloudy->grid_rank - 1) { - std::snprintf(parameter_name, name_bufsize, "Parameter%lld",(q+1)); - } - else { - std::snprintf(parameter_name, name_bufsize, "Temperature"); - } - - attr_id = H5Aopen_name(dset_id, parameter_name); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to open %s attribute in Cooling dataset.\n", - parameter_name); - return GR_FAIL; - } - - double* tmp_param_buf = new double[my_cloudy->grid_dimension[q]]; - status = H5Aread(attr_id, HDF5_R8, tmp_param_buf); - if (attr_id == h5_error) { - std::fprintf(stderr,"Failed to read %s attribute in Cooling dataset.\n", - parameter_name); - return GR_FAIL; - } - - my_cloudy->grid_parameters[q] = new double[my_cloudy->grid_dimension[q]]; - for (long long w = 0LL; w < my_cloudy->grid_dimension[q]; w++) { - if (q < my_cloudy->grid_rank - 1) { - my_cloudy->grid_parameters[q][w] = (double) tmp_param_buf[w]; + // 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]; } - else { - // convert temeperature to log - my_cloudy->grid_parameters[q][w] = (double) std::log10(tmp_param_buf[w]); - } - - } - delete[] tmp_param_buf; - if (grackle_verbose) { - std::fprintf(stdout, - "%s: %g to %g (%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) { - std::fprintf(stderr,"Failed to close %s attribute in Cooling dataset.\n", - parameter_name); - return GR_FAIL; } + my_cloudy->grid_parameters[i] = buf; } - // to make the logic more concise, we do something a *little* silly - // -> we effectively close the dataset and then call a function that reopens - // and closes the dataset - status = H5Dclose(dset_id); - if (status == h5_error) { - std::fprintf(stderr,"Error while to closing %s dataset.\n", dset_name); - return GR_FAIL; + 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; + h5io::drop_GridTableProps(&grid_props); // Read Cooling data. - my_cloudy->data_size = 1; - for (long long q = 0LL; q < my_cloudy->grid_rank; q++) { - my_cloudy->data_size *= my_cloudy->grid_dimension[q]; - } + my_cloudy->data_size = h5io::ArrayShape_elem_count(expected_shape); double* tmp_cool_data = new double[my_cloudy->data_size]; if (grackle_verbose) { std::fprintf(stdout,"Reading from \"%s\" dataset.\n", dset_name); } - if (h5io::read_dataset(file_id, dset_name, tmp_cool_data) != GR_SUCCESS) { + if (h5io::read_dataset(file_id, dset_name, tmp_cool_data, &expected_shape) + != GR_SUCCESS) { delete[] tmp_cool_data; return GR_FAIL; } @@ -254,12 +182,16 @@ int grackle::impl::initialize_cloudy_data( // Read Heating data. if (my_chemistry->UVbackground == 1) { + // Ideally we would parse the attributes describing how the table varies + // and confirm it's consistent with the Cooling table + // -> at this point, this is easy to do, but it introduces overhead std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Heating", group_name); double* tmp_heat_data = new double[my_cloudy->data_size]; - if (h5io::read_dataset(file_id, dset_name, tmp_heat_data) != GR_SUCCESS) { + if (h5io::read_dataset(file_id, dset_name, tmp_heat_data, &expected_shape) + != GR_SUCCESS) { delete[] tmp_heat_data; return GR_FAIL; } @@ -278,10 +210,14 @@ int grackle::impl::initialize_cloudy_data( // Read MMW data. if (my_chemistry->primordial_chemistry == 0 && std::strcmp(group_name, "Primordial") == 0) { + // Ideally we would parse the attributes describing how the table varies + // and confirm it's consistent with the Cooling table + // -> at this point, this is easy to do, but it introduces overhead my_cloudy->mmw_data = new double[my_cloudy->data_size]; if (h5io::read_dataset(file_id, "/CoolingRates/Primordial/MMW", - my_cloudy->mmw_data) != GR_SUCCESS) { + my_cloudy->mmw_data, &expected_shape) + != GR_SUCCESS) { // nothing to cleanup right now return GR_FAIL; } diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 5c681191d..5be2c4f65 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -9,7 +9,9 @@ /// Implements functions to help read from hdf5 files /// //===----------------------------------------------------------------------===// +#include #include +#include #include "hdf5.h" #include "grackle.h" @@ -17,6 +19,135 @@ #include "h5io.hpp" #include "../status_reporting.h" +int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, + char* buffer) { + if (bufsz < 0 || attr_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 = H5Aget_storage_size(attr_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 = H5Aget_type(attr_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 + // - it might be nice if we actually cached this type and carried it around + // so that we can avoid remaking the type every time we load a string + 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; + } + if (is_variable) { + char* tmp_str = nullptr; + if (H5Aread(attr_id, memory_typeid, static_cast(&tmp_str)) < 0) { + // is there a memory leak? + fprintf(stderr, "error in H5Aread for a variable length string\n"); + H5Tclose(memory_typeid); + return -1; + } + std::strncpy(buffer, tmp_str, static_cast(required_bufsz)); + H5free_memory(tmp_str); + } else { + if (H5Aread(attr_id, memory_typeid, buffer) < 0) { + H5Tclose(memory_typeid); + fprintf(stderr, "error in H5Aread\n"); + return -1; + } + } + H5Tclose(memory_typeid); + + if (uses_utf8_encoding) { + // if the string used utf8 encoding on disk, we need to confirm that the + // characters are all compatible with ASCII characterset. + // + // For the uninitiated + // - for the uninitiated, all 128 standard ASCII characters are encoded in + // 7bits and have an identical representation in utf8-encoding. + // - all other utf8 codepoints are represented by 2 or more bytes & + // the leading byte has a value doesn't correspond to an ASCII character + const char max_val = static_cast(127); + for (int i = 0; i < bufsz; i++) { + char chr = buffer[i]; + if (buffer[i] == '\0') { + break; + } 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 + std::fprintf(stderr, + "Error: read a utf8-encoded string from an attribute that " + "contains non-ASCII characters\n"); + return -1; + } + } + } + + return required_bufsz; +} + namespace { // stuff inside an anonymous namespace is local to this file grackle::impl::h5io::ArrayShape shape_from_space(hid_t space_id) { @@ -70,19 +201,23 @@ grackle::impl::h5io::ArrayShape shape_from_space(hid_t space_id) { } } +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 shape_from_space(H5I_INVALID_HID); + 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 shape_from_space(H5I_INVALID_HID); + return mk_invalid_array_shape(); } hid_t space_id = H5Dget_space(dset_id); @@ -136,3 +271,282 @@ int grackle::impl::h5io::read_dataset( return GR_SUCCESS; } + +namespace { // stuff inside an anonymous namespace is local to this file + +static constexpr int attr_name_set_capacity = 30; + +/// Acts as a crude set of attribute names +/// +/// @note +/// if we ever choose to embrace C++, it would be **MUCH** nicer to use +/// std::set +struct AttrNameSet { + char* vals[attr_name_set_capacity]; + int capacity; + int length; +}; + +/// helper function that calculates ArrayShape from the "Rank" and "Dimension" +/// attributes +/// +/// @note +/// we only carry around dset_name argument for nicer error messages +grackle::impl::h5io::ArrayShape shape_from_grid_attrs(hid_t dset_id, + const char* dset_name) { + 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(); + } + + // 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(); + } + } + + // finally, let's format the output + grackle::impl::h5io::ArrayShape out; + 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 +/// +/// returns the number of uniqued attributes accessed by this function +int set_grid_axes_props(hid_t dset_id, const char* dset_name, + grackle::impl::h5io::ArrayShape grid_shape, + grackle::impl::h5io::GridTableAxis* axes) { + using grackle::impl::h5io::read_str_attribute; + int accessed_attr_count = 0; + + for (int i = 0; i < grid_shape.ndim; i++) { + bool is_last_axis = (i + 1) == grid_shape.ndim; + + // Step 1: get the name of the attribute holding the values along axis i + // AND fill axes[i].name with the name of the varying quantity + char val_attr_name[max_attr_name_length]; + std::snprintf(val_attr_name, max_attr_name_length, "Parameter%d", i + 1); + + if (H5Aexists(dset_id, val_attr_name) > 0) { + 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 -1; + } + + 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 -1; + } + + 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 -1; + } + H5Aclose(attr_id); + accessed_attr_count++; + + } else if (is_last_axis && H5Aexists(dset_id, "Temperature") > 0) { + axes[i].name = new char[12]; + std::snprintf(val_attr_name, max_attr_name_length, "Temperature"); + std::snprintf(axes[i].name, 12, "Temperature"); + + } else { + const char* extra_detail = + (is_last_axis) ? "" : " Neither is \"Temperature\"."; + std::fprintf( + stderr, + "Failed to determine attribute associated with axis %d of the " + "\"%s\" dataset. \"%s\" is not a known attribute.%s\n", + i, val_attr_name, dset_name, extra_detail); + return -1; + } + + // step 2: load the values associated with axes[i].name 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 -1; + } + + 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 || axis_shape.shape[0] != grid_shape.shape[i]) { + std::fprintf( + stderr, + "The \"%s\" attr of the \"%s\" dataset isn't an array with a shape " + "that is consistent with the dataset's \"Dimension\" attr.\n", + val_attr_name, dset_name); + return -1; + } + + 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 -1; + } + H5Aclose(attr_id); + accessed_attr_count++; + } + return accessed_attr_count; +} + +int get_num_attrs(hid_t dset_id, const char* dset_name) { + H5O_info1_t info; + if (H5Oget_info2(dset_id, &info, H5O_INFO_NUM_ATTRS) < 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); + } +} + +} // anonymous namespace + +grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( + hid_t file_id, const char* dset_name) { + // 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; + out.table_shape = shape_from_space(H5I_INVALID_HID); + for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { + out.axes[i].name = nullptr; + out.axes[i].values = nullptr; + } + + if (dset_name == nullptr) { // sanity check! + std::fprintf(stderr, "dset_name is a nullptr"); + return out; + } + + 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 out; + } + + // infer the shape of the table from the attributes + grackle::impl::h5io::ArrayShape inferred_shape = + shape_from_grid_attrs(dset_id, dset_name); + if (!ArrayShape_is_valid(inferred_shape)) { + H5Dclose(dset_id); + // shape_from_grid_attrs already printed error messages in this case + return out; + } + int shape_attr_count = 2; + + // parse the quantities along each axis + int axes_prop_attr_count = + set_grid_axes_props(dset_id, dset_name, inferred_shape, out.axes); + if (axes_prop_attr_count < 0) { + H5Dclose(dset_id); + drop_GridTableProps(&out); + // shape_from_grid_attrs already printed error messages in this case + return out; + } + + // perform a check to confirm that we have accessed every available attribute + // -> we may want to disable this check... + int num_attrs = get_num_attrs(dset_id, dset_name); + if (num_attrs < 0) { + H5Dclose(dset_id); + drop_GridTableProps(&out); + // get_num_attrs already printed error messages in this case + return out; + } + + H5Dclose(dset_id); + + int total_accessed_attrs_count = shape_attr_count + axes_prop_attr_count; + if (num_attrs != total_accessed_attrs_count) { + drop_GridTableProps(&out); + std::fprintf( + stderr, + "the number of accessed attributes, %d, for the \"%s\" dataset isn't " + "equal to the number of available attributes, %d. This may be a sign " + "of a problem\n", + total_accessed_attrs_count, dset_name, num_attrs); + return out; + } + + out.table_shape = inferred_shape; // we intentionally do this as the very + // last step! + return out; +} diff --git a/src/clib/utils/h5io.hpp b/src/clib/utils/h5io.hpp index 767884818..7bc64f899 100644 --- a/src/clib/utils/h5io.hpp +++ b/src/clib/utils/h5io.hpp @@ -20,6 +20,30 @@ namespace grackle::impl::h5io { +/// copies the string encoded in the specified hdf5 attribute and into +/// ``buffer`` as a null-terminated string, and returns ``min_req_bufsz``. +/// +/// ``min_req_bufsz`` is the minimum required ``bufsz`` (i.e. ``bufsz`` is the +/// the length of ``buffer``) where this function will try to load the +/// attribute. +/// - this is the maximum length of the string (including the null character). +/// To put it another way, after this function succeeds, std::strlen(buffer), +/// will return a ``x`` that is bounded by ``0 <= x <= 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 ``bufsz`` is smaller than ``min_req_bufsz``, unless +/// unless ``bufsz`` is zero. In the event that ``bufsz`` is zero, nothing is +/// written and ``buffer`` may be a ``nullptr``. In this case, the function +/// returns ``min_req_bufsz``. If this function fails, a negative value is +/// returned. +/// +/// @note +/// The function reports an error if it reads a utf8-encoded string that +/// contains non-ASCII characters +int read_str_attribute(hid_t attr_id, int bufsz, char* buffer); + /// represents a contiguous array shape /// /// @note @@ -39,6 +63,19 @@ 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 @@ -64,6 +101,37 @@ ArrayShape read_dataset_shape(hid_t file_id, const char* dset_name); 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); +} + +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 +GridTableProps parse_GridTableProps(hid_t file_id, const char* dset_name); + } // namespace grackle::impl::h5io #endif /* UTILS_H5IO_HPP */ From 62ac9c55cd11f7a387b59d02b6d2c16bdc07c2f0 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 26 Nov 2025 17:47:12 -0500 Subject: [PATCH 26/38] remove unused variables --- src/clib/initialize_cloudy_data.cpp | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 1c6e4a53c..b56062578 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -51,7 +51,6 @@ int grackle::impl::initialize_cloudy_data( code_units *my_units, int read_data) { - long long temp_int; char dset_name[MAX_PARAMETER_NAME_LENGTH]; const std::size_t name_bufsize = static_cast(MAX_PARAMETER_NAME_LENGTH); @@ -91,12 +90,8 @@ int grackle::impl::initialize_cloudy_data( ) / (std::pow(tbase1, 3.0) * 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); + 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; @@ -224,7 +219,7 @@ int grackle::impl::initialize_cloudy_data( } - status = H5Fclose (file_id); + H5Fclose (file_id); if (my_cloudy->grid_rank > GRACKLE_CLOUDY_TABLE_MAX_DIMENSION) { std::fprintf( From 753eb9c8503522ead0b9f695a0094216c079be4d Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 26 Nov 2025 18:07:01 -0500 Subject: [PATCH 27/38] fixup get_num_attrs --- src/clib/utils/h5io.cpp | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 5be2c4f65..9522c5cb1 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -130,7 +130,6 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, // the leading byte has a value doesn't correspond to an ASCII character const char max_val = static_cast(127); for (int i = 0; i < bufsz; i++) { - char chr = buffer[i]; if (buffer[i] == '\0') { break; } else if (buffer[i] < '\0' || buffer[i] > max_val) { @@ -464,8 +463,17 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, } int get_num_attrs(hid_t dset_id, const char* dset_name) { - H5O_info1_t info; - if (H5Oget_info2(dset_id, &info, H5O_INFO_NUM_ATTRS) < 0) { +#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; @@ -476,6 +484,7 @@ int get_num_attrs(hid_t dset_id, const char* dset_name) { } else { return static_cast(info.num_attrs); } +#endif } } // anonymous namespace @@ -522,10 +531,14 @@ grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( return out; } + int total_accessed_attrs_count = shape_attr_count + axes_prop_attr_count; + // perform a check to confirm that we have accessed every available attribute // -> we may want to disable this check... int num_attrs = get_num_attrs(dset_id, dset_name); - if (num_attrs < 0) { + 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 @@ -534,7 +547,6 @@ grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( H5Dclose(dset_id); - int total_accessed_attrs_count = shape_attr_count + axes_prop_attr_count; if (num_attrs != total_accessed_attrs_count) { drop_GridTableProps(&out); std::fprintf( From fb25a6b28392549e32e0eeeb62cd7cbfff718626 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 30 Nov 2025 16:47:33 -0500 Subject: [PATCH 28/38] more meaningful error message --- src/clib/utils/h5io.cpp | 371 +++++++++++++++++++++++++++++++++------- src/clib/utils/h5io.hpp | 11 +- 2 files changed, 315 insertions(+), 67 deletions(-) diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 9522c5cb1..c33baca7d 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -273,26 +273,159 @@ int grackle::impl::h5io::read_dataset( namespace { // stuff inside an anonymous namespace is local to this file -static constexpr int attr_name_set_capacity = 30; - -/// Acts as a crude set of attribute names +/// @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 AttrNameSet { - char* vals[attr_name_set_capacity]; - int capacity; +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 /// -/// @note -/// we only carry around dset_name argument for nicer error messages -grackle::impl::h5io::ArrayShape shape_from_grid_attrs(hid_t dset_id, - const char* dset_name) { +/// @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, @@ -318,6 +451,13 @@ grackle::impl::h5io::ArrayShape shape_from_grid_attrs(hid_t dset_id, 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]; @@ -350,6 +490,13 @@ grackle::impl::h5io::ArrayShape shape_from_grid_attrs(hid_t dset_id, 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 grackle::impl::h5io::ArrayShape out; @@ -364,10 +511,19 @@ static constexpr int max_attr_name_length = 64; /// updates axes with parsed parameters /// -/// returns the number of uniqued attributes accessed by this function +/// @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) { + grackle::impl::h5io::GridTableAxis* axes, + AttrNameRecorder* name_recorder) { using grackle::impl::h5io::read_str_attribute; int accessed_attr_count = 0; @@ -389,7 +545,7 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, std::fprintf(stderr, "Failed to open \"%s\" attribute of \"%s\" dataset.\n", tmp_attr_name, dset_name); - return -1; + return GR_FAIL; } int min_buf_length = read_str_attribute(attr_id, 0, nullptr); @@ -398,7 +554,7 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, stderr, "can't determine string size held by \"%s\" attr of \"%s\" dset.\n", tmp_attr_name, dset_name); - return -1; + return GR_FAIL; } axes[i].name = new char[min_buf_length]; @@ -406,10 +562,17 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, std::fprintf(stderr, "error loading the \"%s\" attr from \"%s\" dataset\n", tmp_attr_name, dset_name); - return -1; + return GR_FAIL; } H5Aclose(attr_id); - accessed_attr_count++; + + 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; + } } else if (is_last_axis && H5Aexists(dset_id, "Temperature") > 0) { axes[i].name = new char[12]; @@ -424,7 +587,7 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, "Failed to determine attribute associated with axis %d of the " "\"%s\" dataset. \"%s\" is not a known attribute.%s\n", i, val_attr_name, dset_name, extra_detail); - return -1; + return GR_FAIL; } // step 2: load the values associated with axes[i].name and store them in @@ -434,7 +597,7 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_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 -1; + return GR_FAIL; } hid_t space_id = H5Aget_space(attr_id); @@ -446,7 +609,7 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, "The \"%s\" attr of the \"%s\" dataset isn't an array with a shape " "that is consistent with the dataset's \"Dimension\" attr.\n", val_attr_name, dset_name); - return -1; + return GR_FAIL; } axes[i].values = new double[grid_shape.shape[i]]; @@ -454,12 +617,72 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, H5Aclose(attr_id); std::fprintf(stderr, "failed to read \"%s\" attr of \"%s\" dataset.\n", val_attr_name, dset_name); - return -1; + return GR_FAIL; } H5Aclose(attr_id); - accessed_attr_count++; + 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 accessed_attr_count; + 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) { @@ -491,50 +714,36 @@ int get_num_attrs(hid_t dset_id, const char* dset_name) { grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( hid_t file_id, const char* dset_name) { - // 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; - out.table_shape = shape_from_space(H5I_INVALID_HID); - for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { - out.axes[i].name = nullptr; - out.axes[i].values = nullptr; - } - if (dset_name == nullptr) { // sanity check! std::fprintf(stderr, "dset_name is a nullptr"); - return out; + 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 out; + return mk_invalid_GridTableProps(); } - // infer the shape of the table from the attributes - grackle::impl::h5io::ArrayShape inferred_shape = - shape_from_grid_attrs(dset_id, dset_name); - if (!ArrayShape_is_valid(inferred_shape)) { - H5Dclose(dset_id); - // shape_from_grid_attrs already printed error messages in this case - return out; - } - int shape_attr_count = 2; + // 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); - // parse the quantities along each axis - int axes_prop_attr_count = - set_grid_axes_props(dset_id, dset_name, inferred_shape, out.axes); - if (axes_prop_attr_count < 0) { + if (!GridTableProps_is_valid(out)) { H5Dclose(dset_id); - drop_GridTableProps(&out); - // shape_from_grid_attrs already printed error messages in this case + // parse_GridTableProps_helper already printed appropriate errors messages return out; } - int total_accessed_attrs_count = shape_attr_count + axes_prop_attr_count; - - // perform a check to confirm that we have accessed every available attribute - // -> we may want to disable this check... + // 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; @@ -543,22 +752,52 @@ grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( drop_GridTableProps(&out); // get_num_attrs already printed error messages in this case return out; - } - - H5Dclose(dset_id); - - if (num_attrs != total_accessed_attrs_count) { + } else if (num_attrs != total_accessed_attrs_count) { drop_GridTableProps(&out); - std::fprintf( - stderr, - "the number of accessed attributes, %d, for the \"%s\" dataset isn't " - "equal to the number of available attributes, %d. This may be a sign " - "of a problem\n", - total_accessed_attrs_count, dset_name, num_attrs); - return 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(); } - out.table_shape = inferred_shape; // we intentionally do this as the very - // last step! + H5Dclose(dset_id); return out; } diff --git a/src/clib/utils/h5io.hpp b/src/clib/utils/h5io.hpp index 7bc64f899..528e928c9 100644 --- a/src/clib/utils/h5io.hpp +++ b/src/clib/utils/h5io.hpp @@ -57,7 +57,7 @@ struct ArrayShape { /// checks whether shape is valid inline bool ArrayShape_is_valid(ArrayShape shape) { return shape.ndim >= -1; } -/// checks whether shape is valid +/// checks whether shape refers to a scalar inline bool ArrayShape_is_scalar(ArrayShape shape) { return shape.ndim == 0; } /// checks whether shape is null @@ -116,6 +116,7 @@ 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) { @@ -130,6 +131,14 @@ inline void drop_GridTableProps(GridTableProps* ptr) { } /// parses the GridTableProps from dataset attributes +/// +/// @param[in] file_id File identifier +/// @param[in] dset_name The name of the dataset that the parameters are +/// attached to. +/// +/// @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); } // namespace grackle::impl::h5io From 0859fb86c9d644ad292d33bcde994ceebbd0e4bf Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 1 Dec 2025 08:04:49 -0500 Subject: [PATCH 29/38] address compiler warnings --- src/clib/initialize_UVbackground_data.cpp | 2 +- src/clib/utils/h5io.cpp | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index b34fb6aa2..d9af11420 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -63,7 +63,7 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // Read in UV background data from hdf5 file. - hid_t file_id, dset_id, dspace_id; + hid_t file_id, dset_id; herr_t status; herr_t h5_error = -1; diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index c33baca7d..337f085e1 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -499,7 +499,9 @@ grackle::impl::h5io::ArrayShape shape_from_grid_attrs( } // finally, let's format the output - grackle::impl::h5io::ArrayShape out; + // -> 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]; @@ -525,7 +527,6 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, grackle::impl::h5io::GridTableAxis* axes, AttrNameRecorder* name_recorder) { using grackle::impl::h5io::read_str_attribute; - int accessed_attr_count = 0; for (int i = 0; i < grid_shape.ndim; i++) { bool is_last_axis = (i + 1) == grid_shape.ndim; From 34f5a33dcf73496a7d85d5992f6a70570e2093b5 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 1 Dec 2025 08:12:28 -0500 Subject: [PATCH 30/38] move a bunch of hdf5-specific macros out of grackle_macros.h --- src/clib/grackle_macros.h | 14 -------------- src/clib/initialize_UVbackground_data.cpp | 2 +- src/clib/initialize_cloudy_data.cpp | 2 +- src/clib/utils/h5io.cpp | 15 +++++++++++++++ 4 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/clib/grackle_macros.h b/src/clib/grackle_macros.h index 995965cef..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 diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index d9af11420..84b8a1917 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -15,7 +15,7 @@ #include #include "hdf5.h" #include "grackle.h" -#include "grackle_macros.h" +#include "grackle_macros.h" // FLOAT_UNDEFINED #include "utils/h5io.hpp" #include "initialize_UVbackground_data.hpp" diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index b56062578..4f4901e55 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -16,7 +16,7 @@ #include // std::strcmp #include "hdf5.h" #include "grackle.h" -#include "grackle_macros.h" // HDF5_I8, HDF5_R8, TRUE +#include "grackle_macros.h" // TRUE #include "initialize_cloudy_data.hpp" #include "utils/h5io.hpp" diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 337f085e1..1558dd571 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -19,6 +19,21 @@ #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 + int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, char* buffer) { if (bufsz < 0 || attr_id == H5I_INVALID_HID) { From 75b687315be0dda948616390eeb3a6c2b6fbad7d Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 1 Dec 2025 08:15:32 -0500 Subject: [PATCH 31/38] add a few docstrings --- src/clib/utils/h5io.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 1558dd571..45c92cdee 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -8,6 +8,8 @@ /// @file /// Implements functions to help read from hdf5 files /// +/// The error-handling in these functions leaves a lot to be desired... +/// //===----------------------------------------------------------------------===// #include #include @@ -164,6 +166,7 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, 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 @@ -215,6 +218,10 @@ grackle::impl::h5io::ArrayShape shape_from_space(hid_t space_id) { } } +/// 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); } From c11ddea33c0ede6e7021ab826e0424663798cc1b Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 1 Dec 2025 08:46:14 -0500 Subject: [PATCH 32/38] light refactor of read_str_attribute --- src/clib/utils/h5io.cpp | 65 ++++++++++++++++++++++++++++------------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 45c92cdee..ea8209cc7 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -36,6 +36,42 @@ #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) +} + +} // anonymous namespace + int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, char* buffer) { if (bufsz < 0 || attr_id == H5I_INVALID_HID) { @@ -137,27 +173,14 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, H5Tclose(memory_typeid); if (uses_utf8_encoding) { - // if the string used utf8 encoding on disk, we need to confirm that the - // characters are all compatible with ASCII characterset. - // - // For the uninitiated - // - for the uninitiated, all 128 standard ASCII characters are encoded in - // 7bits and have an identical representation in utf8-encoding. - // - all other utf8 codepoints are represented by 2 or more bytes & - // the leading byte has a value doesn't correspond to an ASCII character - const char max_val = static_cast(127); - for (int i = 0; i < bufsz; i++) { - if (buffer[i] == '\0') { - break; - } 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 - std::fprintf(stderr, - "Error: read a utf8-encoded string from an attribute that " - "contains non-ASCII characters\n"); - return -1; - } + // 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; } } From fdc3c72508b1a4fa99d38a6c96148b5a2c0e2013 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 1 Dec 2025 10:09:49 -0500 Subject: [PATCH 33/38] wrap the logic for loading a string dataset --- src/clib/initialize_UVbackground_data.cpp | 37 ++++------- src/clib/utils/h5io.cpp | 64 ++++++++++++++----- src/clib/utils/h5io.hpp | 78 ++++++++++++++++++----- 3 files changed, 122 insertions(+), 57 deletions(-) diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index 84b8a1917..0ce07c74c 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -62,43 +62,28 @@ int grackle::impl::initialize_UVbackground_data(chemistry_data *my_chemistry, // Read in UV background data from hdf5 file. - - hid_t file_id, dset_id; - herr_t status; - herr_t h5_error = -1; - if (grackle_verbose) std::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); + hid_t 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) { - std::fprintf(stderr, "Can't open 'Info' dataset in %s.\n", - my_chemistry->grackle_data_file); + 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; } - - int strlen = (int)(H5Dget_storage_size(dset_id)); - std::vector 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.data()); - if (status == h5_error) { - std::fprintf(stderr, "Failed to read dataset 'Info'.\n"); + 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; } - H5Tclose(memtype); - H5Dclose(dset_id); - // Open redshift dataset and get number of elements const h5io::ArrayShape common_shape diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index ea8209cc7..1d56f439a 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -70,11 +70,9 @@ bool is_ascii_string(const char* buffer, int bufsz) { return false; // (buffer doesn't contain a null character) } -} // anonymous namespace - -int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, - char* buffer) { - if (bufsz < 0 || attr_id == H5I_INVALID_HID) { +/// 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; } @@ -82,7 +80,8 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, // -> 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 = H5Aget_storage_size(attr_id); + 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, @@ -104,7 +103,7 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, // ASCII before we return bool is_variable, uses_utf8_encoding; { - hid_t disk_typeid = H5Aget_type(attr_id); + 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; @@ -135,8 +134,6 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, // construct the disk_typeid that we use for describing `buffer` to the hdf5 // library - // - it might be nice if we actually cached this type and carried it around - // so that we can avoid remaking the type every time we load a string 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"); @@ -153,20 +150,32 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, H5Tclose(memory_typeid); return -1; } + + // actually read the data if (is_variable) { char* tmp_str = nullptr; - if (H5Aread(attr_id, memory_typeid, static_cast(&tmp_str)) < 0) { - // is there a memory leak? - fprintf(stderr, "error in H5Aread for a variable length string\n"); + 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)); - H5free_memory(tmp_str); } else { - if (H5Aread(attr_id, memory_typeid, buffer) < 0) { + 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 H5Aread\n"); + fprintf(stderr, "error in %s for a fixed length string\n", + (is_attr) ? "H5Aread" : "H5Dread"); return -1; } } @@ -187,6 +196,31 @@ int grackle::impl::h5io::read_str_attribute(hid_t attr_id, int bufsz, 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 diff --git a/src/clib/utils/h5io.hpp b/src/clib/utils/h5io.hpp index 528e928c9..9832aa6c9 100644 --- a/src/clib/utils/h5io.hpp +++ b/src/clib/utils/h5io.hpp @@ -20,30 +20,77 @@ namespace grackle::impl::h5io { -/// copies the string encoded in the specified hdf5 attribute and into -/// ``buffer`` as a null-terminated string, and returns ``min_req_bufsz``. +/// copies the string encoded in the specified hdf5 dataset into ``buffer`` as +/// a null-terminated string, and returns ``min_req_bufsz`` (if successful). /// -/// ``min_req_bufsz`` is the minimum required ``bufsz`` (i.e. ``bufsz`` is the -/// the length of ``buffer``) where this function will try to load the -/// attribute. +/// @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). -/// To put it another way, after this function succeeds, std::strlen(buffer), -/// will return a ``x`` that is bounded by ``0 <= x <= min_req_bufsz``. +/// 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 ``bufsz`` is smaller than ``min_req_bufsz``, unless -/// unless ``bufsz`` is zero. In the event that ``bufsz`` is zero, nothing is -/// written and ``buffer`` may be a ``nullptr``. In this case, the function -/// returns ``min_req_bufsz``. If this function fails, a negative value is -/// returned. +/// 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 -/// The function reports an error if it reads a utf8-encoded string that -/// contains non-ASCII characters +/// 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 @@ -133,8 +180,7 @@ inline void drop_GridTableProps(GridTableProps* ptr) { /// parses the GridTableProps from dataset attributes /// /// @param[in] file_id File identifier -/// @param[in] dset_name The name of the dataset that the parameters are -/// attached to. +/// @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 From fa4f0e8cb48ca007999b280c0305d801eb7dfb58 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 2 Dec 2025 13:28:23 -0500 Subject: [PATCH 34/38] add a consistency check --- src/clib/utils/h5io.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 1d56f439a..d968bdfd2 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -878,6 +878,22 @@ grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( 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_same(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; } From 4454a61ef1bfbaa072871cef354d79d4ce3d58ba Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 2 Dec 2025 15:53:29 -0500 Subject: [PATCH 35/38] light refactor --- src/clib/utils/h5io.cpp | 49 ++++++++++++++++++----------------------- 1 file changed, 22 insertions(+), 27 deletions(-) diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index d968bdfd2..6bc61ba8c 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -607,15 +607,20 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, AttrNameRecorder* name_recorder) { using grackle::impl::h5io::read_str_attribute; - for (int i = 0; i < grid_shape.ndim; i++) { - bool is_last_axis = (i + 1) == grid_shape.ndim; + bool legacy_mode = H5Aexists(dset_id, "Temperature") > 0; - // Step 1: get the name of the attribute holding the values along axis i - // AND fill axes[i].name with the name of the varying quantity + 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]; - std::snprintf(val_attr_name, max_attr_name_length, "Parameter%d", i + 1); - if (H5Aexists(dset_id, val_attr_name) > 0) { + 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); @@ -636,7 +641,6 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, 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, @@ -654,25 +658,11 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, return GR_FAIL; } - } else if (is_last_axis && H5Aexists(dset_id, "Temperature") > 0) { - axes[i].name = new char[12]; - std::snprintf(val_attr_name, max_attr_name_length, "Temperature"); - std::snprintf(axes[i].name, 12, "Temperature"); - - } else { - const char* extra_detail = - (is_last_axis) ? "" : " Neither is \"Temperature\"."; - std::fprintf( - stderr, - "Failed to determine attribute associated with axis %d of the " - "\"%s\" dataset. \"%s\" is not a known attribute.%s\n", - i, val_attr_name, dset_name, extra_detail); - return GR_FAIL; + std::snprintf(val_attr_name, max_attr_name_length, "Parameter%d", i + 1); } - // step 2: load the values associated with axes[i].name and store them in + // 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", @@ -683,12 +673,17 @@ int set_grid_axes_props(hid_t dset_id, const char* dset_name, 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 || axis_shape.shape[0] != grid_shape.shape[i]) { + 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\" attr of the \"%s\" dataset isn't an array with a shape " - "that is consistent with the dataset's \"Dimension\" attr.\n", - val_attr_name, dset_name); + "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; } From 88e56f853f665c6c608eb6b31973b03f2c71993f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 3 Dec 2025 01:14:08 -0500 Subject: [PATCH 36/38] a last change --- src/clib/initialize_cloudy_data.cpp | 96 ++++++++++++++++++----------- src/clib/utils/h5io.cpp | 54 +++++++++++++++- src/clib/utils/h5io.hpp | 15 ++++- 3 files changed, 125 insertions(+), 40 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 4f4901e55..867845642 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -40,6 +40,35 @@ void initialize_empty_cloudy_data_struct(cloudy_data *my_cloudy) 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 @@ -150,30 +179,17 @@ int grackle::impl::initialize_cloudy_data( fprintf(stdout, "}\n"); } h5io::ArrayShape expected_shape = grid_props.table_shape; - h5io::drop_GridTableProps(&grid_props); // Read Cooling data. my_cloudy->data_size = h5io::ArrayShape_elem_count(expected_shape); - double* tmp_cool_data = new double[my_cloudy->data_size]; - if (grackle_verbose) { - std::fprintf(stdout,"Reading from \"%s\" dataset.\n", dset_name); - } - if (h5io::read_dataset(file_id, dset_name, tmp_cool_data, &expected_shape) - != GR_SUCCESS) { - delete[] tmp_cool_data; - return GR_FAIL; - } - - my_cloudy->cooling_data = new double[my_cloudy->data_size]; - for (long long q = 0LL; q < my_cloudy->data_size; q++) { - my_cloudy->cooling_data[q] = tmp_cool_data[q] > 0 ? - (double) std::log10(tmp_cool_data[q]) : (double) SMALL_LOG_VALUE; + my_cloudy->cooling_data = load_heatcool_data(file_id, dset_name, CoolUnit, + expected_shape); - // Convert to code units. - my_cloudy->cooling_data[q] -= std::log10(CoolUnit); + if (my_cloudy->cooling_data == nullptr) { + h5io::drop_GridTableProps(&grid_props); + H5Fclose(file_id); + return GR_FAIL; } - delete[] tmp_cool_data; - // Read Heating data. if (my_chemistry->UVbackground == 1) { @@ -184,22 +200,19 @@ int grackle::impl::initialize_cloudy_data( std::snprintf(dset_name, name_bufsize, "/CoolingRates/%s/Heating", group_name); - double* tmp_heat_data = new double[my_cloudy->data_size]; - if (h5io::read_dataset(file_id, dset_name, tmp_heat_data, &expected_shape) - != GR_SUCCESS) { - delete[] tmp_heat_data; + //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; } - - my_cloudy->heating_data = new double[my_cloudy->data_size]; - for (long long q = 0LL; q < my_cloudy->data_size; q++) { - my_cloudy->heating_data[q] = tmp_heat_data[q] > 0 ? - (double) std::log10(tmp_heat_data[q]) : (double) SMALL_LOG_VALUE; - - // Convert to code units. - my_cloudy->heating_data[q] -= std::log10(CoolUnit); - } - delete[] tmp_heat_data; } // Read MMW data. @@ -209,11 +222,20 @@ int grackle::impl::initialize_cloudy_data( // and confirm it's consistent with the Cooling table // -> at this point, this is easy to do, but it introduces overhead + const char* mmw_dset_name = "/CoolingRates/Primordial/MMW"; + + //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, "/CoolingRates/Primordial/MMW", - my_cloudy->mmw_data, &expected_shape) - != GR_SUCCESS) { - // nothing to cleanup right now + 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; } diff --git a/src/clib/utils/h5io.cpp b/src/clib/utils/h5io.cpp index 6bc61ba8c..f235031b6 100644 --- a/src/clib/utils/h5io.cpp +++ b/src/clib/utils/h5io.cpp @@ -331,7 +331,7 @@ int grackle::impl::h5io::read_dataset( ArrayShape actual_shape = shape_from_space(space_id); H5Sclose(space_id); - if (!ArrayShape_is_same(*expected_shape, actual_shape)) { + if (!ArrayShape_is_equal(*expected_shape, actual_shape)) { H5Dclose(dset_id); std::fprintf(stderr, "The \"%s\" dataset has an unexpected shape.\n", dset_name); @@ -878,7 +878,7 @@ grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_GridTableProps( ArrayShape actual_shape = shape_from_space(space_id); H5Sclose(space_id); if (!ArrayShape_is_null(actual_shape) && - !ArrayShape_is_same(out.table_shape, actual_shape)) { + !ArrayShape_is_equal(out.table_shape, actual_shape)) { drop_GridTableProps(&out); H5Dclose(dset_id); std::fprintf( @@ -892,3 +892,53 @@ grackle::impl::h5io::GridTableProps grackle::impl::h5io::parse_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/utils/h5io.hpp b/src/clib/utils/h5io.hpp index 9832aa6c9..80748ff64 100644 --- a/src/clib/utils/h5io.hpp +++ b/src/clib/utils/h5io.hpp @@ -126,7 +126,7 @@ inline std::int64_t ArrayShape_elem_count(ArrayShape shape) { /// checks whether shape_a and shape_b are the same /// /// returns false if the either shape is invalid -inline bool ArrayShape_is_same(ArrayShape shape_a, ArrayShape shape_b) { +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; } @@ -187,6 +187,19 @@ inline void drop_GridTableProps(GridTableProps* ptr) { /// 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 */ From 34bdcf64181fc315debc66ad212184d10f67d50c Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Wed, 3 Dec 2025 20:39:18 -0500 Subject: [PATCH 37/38] adjust a comment --- src/clib/initialize_cloudy_data.cpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index 867845642..e761a8d2d 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -193,13 +193,16 @@ int grackle::impl::initialize_cloudy_data( // Read Heating data. if (my_chemistry->UVbackground == 1) { - // Ideally we would parse the attributes describing how the table varies - // and confirm it's consistent with the Cooling table - // -> at this point, this is easy to do, but it introduces overhead - 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); @@ -218,12 +221,15 @@ int grackle::impl::initialize_cloudy_data( // Read MMW data. if (my_chemistry->primordial_chemistry == 0 && std::strcmp(group_name, "Primordial") == 0) { - // Ideally we would parse the attributes describing how the table varies - // and confirm it's consistent with the Cooling table - // -> at this point, this is easy to do, but it introduces overhead - 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); From b6966a2a035c9099c620ff0d9f38c9bb79a78f12 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 12 Dec 2025 12:13:49 -0500 Subject: [PATCH 38/38] move utils/h5io.[ch]pp -> support/h5io.[ch]pp --- src/clib/CMakeLists.txt | 2 +- src/clib/Make.config.objects | 2 +- src/clib/initialize_UVbackground_data.cpp | 2 +- src/clib/initialize_cloudy_data.cpp | 2 +- src/clib/{utils => support}/h5io.cpp | 0 src/clib/{utils => support}/h5io.hpp | 0 6 files changed, 4 insertions(+), 4 deletions(-) rename src/clib/{utils => support}/h5io.cpp (100%) rename src/clib/{utils => support}/h5io.hpp (100%) diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 7bcb5ecaa..c6f52dc03 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -124,8 +124,8 @@ 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/h5io.cpp utils/h5io.hpp utils-cpp.cpp utils-cpp.hpp utils-field.hpp fortran_func_wrappers.hpp diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 858d7d1ee..a762a8537 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -50,7 +50,7 @@ OBJS_CONFIG_LIB = \ rate_functions.lo \ rate_utils.lo \ gaussj_g.lo \ - utils/h5io.lo \ + support/h5io.lo \ utils.lo \ utils-cpp.lo \ internal_types.lo diff --git a/src/clib/initialize_UVbackground_data.cpp b/src/clib/initialize_UVbackground_data.cpp index 0ce07c74c..fafcfb813 100644 --- a/src/clib/initialize_UVbackground_data.cpp +++ b/src/clib/initialize_UVbackground_data.cpp @@ -16,7 +16,7 @@ #include "hdf5.h" #include "grackle.h" #include "grackle_macros.h" // FLOAT_UNDEFINED -#include "utils/h5io.hpp" +#include "support/h5io.hpp" #include "initialize_UVbackground_data.hpp" diff --git a/src/clib/initialize_cloudy_data.cpp b/src/clib/initialize_cloudy_data.cpp index e761a8d2d..fa0aab74d 100644 --- a/src/clib/initialize_cloudy_data.cpp +++ b/src/clib/initialize_cloudy_data.cpp @@ -18,7 +18,7 @@ #include "grackle.h" #include "grackle_macros.h" // TRUE #include "initialize_cloudy_data.hpp" -#include "utils/h5io.hpp" +#include "support/h5io.hpp" #define SMALL_LOG_VALUE (-99.0) #define MAX_PARAMETER_NAME_LENGTH (512) diff --git a/src/clib/utils/h5io.cpp b/src/clib/support/h5io.cpp similarity index 100% rename from src/clib/utils/h5io.cpp rename to src/clib/support/h5io.cpp diff --git a/src/clib/utils/h5io.hpp b/src/clib/support/h5io.hpp similarity index 100% rename from src/clib/utils/h5io.hpp rename to src/clib/support/h5io.hpp