diff options
Diffstat (limited to 'src/nuc_eos/readtable.c')
-rw-r--r-- | src/nuc_eos/readtable.c | 139 |
1 files changed, 139 insertions, 0 deletions
diff --git a/src/nuc_eos/readtable.c b/src/nuc_eos/readtable.c new file mode 100644 index 0000000..72077e1 --- /dev/null +++ b/src/nuc_eos/readtable.c @@ -0,0 +1,139 @@ +#include <stdlib.h> + +#include "hdf5.h" + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +// Interface of function for fortran module eostable +void CCTK_FNAME(allocate_eosmodule) + (CCTK_INT* , CCTK_INT* , CCTK_INT*, CCTK_REAL*, + CCTK_REAL*, CCTK_REAL*, CCTK_REAL*, CCTK_REAL*); + +// Catch HDF5 errors +#define HDF5_ERROR(fn_call) \ + do { \ + int _error_code = fn_call; \ + if (_error_code < 0) \ + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, \ + "HDF5 call '%s' returned error code %d", \ + #fn_call, _error_code); \ + } while (0) + +int file_is_readable(const char* filename); +int file_is_readable(const char* filename) +{ + FILE* fp = NULL; + fp = fopen(filename, "r"); + if(fp != NULL) + { + fclose(fp); + return 1; + } + return 0; +} + +// Cactus calls this function. It reads in the table and calls a fortran +// function to setup values for the fortran eos module +void EOS_OMNI_ReadTable(CCTK_ARGUMENTS); +void EOS_OMNI_ReadTable(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + + hid_t file; + if (!file_is_readable(nuceos_table_name)) + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Could not read nuceos_table_name'%s'", + nuceos_table_name); + HDF5_ERROR(file = H5Fopen(nuceos_table_name, H5F_ACC_RDONLY, H5P_DEFAULT)); + +// Use these two defines to easily read in a lot of variables in the same way +// The first reads in one variable of a given type completely +#define READ_EOS_HDF5(NAME,VAR,TYPE,MEM) \ + do { \ + hid_t dataset; \ + HDF5_ERROR(dataset = H5Dopen(file, NAME, H5P_DEFAULT)); \ + HDF5_ERROR(H5Dread(dataset, TYPE, MEM, H5S_ALL, H5P_DEFAULT, VAR)); \ + HDF5_ERROR(H5Dclose(dataset)); \ + } while (0) +// The second reads a given variable into a hyperslab of the alltables array +#define READ_EOSTABLE_HDF5(NAME,OFF) \ + do { \ + hsize_t offset[2] = {OFF,0}; \ + H5Sselect_hyperslab(mem3, H5S_SELECT_SET, offset, NULL, var3, NULL); \ + READ_EOS_HDF5(NAME,alltables,H5T_NATIVE_DOUBLE,mem3); \ + } while (0) + + // Read size of tables + CCTK_INT nrho, ntemp, nye; + + READ_EOS_HDF5("pointsrho", &nrho, H5T_NATIVE_INT, H5S_ALL); + READ_EOS_HDF5("pointstemp", &ntemp, H5T_NATIVE_INT, H5S_ALL); + READ_EOS_HDF5("pointsye", &nye, H5T_NATIVE_INT, H5S_ALL); + + // Allocate memory for tables + CCTK_REAL *alltables, *logrho, *logtemp, *ye, energy_shift; + + if (!(alltables = (CCTK_REAL*)malloc(nrho * ntemp * nye * 19 * sizeof(CCTK_REAL)))) + CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n"); + if (!(logrho = (CCTK_REAL*)malloc(nrho * sizeof(CCTK_REAL)))) + CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n"); + if (!(logtemp = (CCTK_REAL*)malloc(ntemp * sizeof(CCTK_REAL)))) + CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n"); + if (!(ye = (CCTK_REAL*)malloc(nye * sizeof(CCTK_REAL)))) + CCTK_WARN(CCTK_WARN_ABORT, "Cannot allocate memory for EOS table\n"); + + // Prepare HDF5 to read hyperslabs into alltables + hsize_t table_dims[2] = {19, nrho * ntemp * nye}; + hsize_t var3[2] = { 1, nrho * ntemp * nye}; + hid_t mem3 = H5Screate_simple(2, table_dims, NULL); + + // Read alltables + READ_EOSTABLE_HDF5("logpress", 0); + READ_EOSTABLE_HDF5("logenergy", 1); + READ_EOSTABLE_HDF5("entropy", 2); + READ_EOSTABLE_HDF5("munu", 3); + READ_EOSTABLE_HDF5("cs2", 4); + READ_EOSTABLE_HDF5("dedt", 5); + READ_EOSTABLE_HDF5("dpdrhoe", 6); + READ_EOSTABLE_HDF5("dpderho", 7); + // chemical potentials + READ_EOSTABLE_HDF5("muhat", 8); + READ_EOSTABLE_HDF5("mu_e", 9); + READ_EOSTABLE_HDF5("mu_p", 10); + READ_EOSTABLE_HDF5("mu_n", 11); + // compositions + READ_EOSTABLE_HDF5("Xa", 12); + READ_EOSTABLE_HDF5("Xh", 13); + READ_EOSTABLE_HDF5("Xn", 14); + READ_EOSTABLE_HDF5("Xp", 15); + // average nucleus + READ_EOSTABLE_HDF5("Abar", 16); + READ_EOSTABLE_HDF5("Zbar", 17); + // Gamma + READ_EOSTABLE_HDF5("gamma", 18); + + // Read additional tables and variables + READ_EOS_HDF5("logrho", logrho, H5T_NATIVE_DOUBLE, H5S_ALL); + READ_EOS_HDF5("logtemp", logtemp, H5T_NATIVE_DOUBLE, H5S_ALL); + READ_EOS_HDF5("ye", ye, H5T_NATIVE_DOUBLE, H5S_ALL); + READ_EOS_HDF5("energy_shift", &energy_shift, H5T_NATIVE_DOUBLE, H5S_ALL); + + H5Fclose(file); + + // Give all values to fortran - which will copy them until I can find out how + // I can tell Fortran to use those pointers inside the module + CCTK_FNAME(allocate_eosmodule) + (&nrho, &ntemp, &nye, alltables, logrho, logtemp, ye, &energy_shift); + + // Free the memory again because fortran copied the whole thing now + free(ye); + free(logtemp); + free(logrho); + free(alltables); + +} + |