#include #define H5_USE_16_API 1 #include "hdf5.h" #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" // mini NoMPI #ifdef HAVE_CAPABILITY_MPI #include #define BCAST(buffer, size) MPI_Bcast(buffer, size, MPI_BYTE, my_reader_process, MPI_COMM_WORLD) #else #define BCAST(buffer, size) do { /* do nothing */ } while(0) #endif // Interface of function for fortran module eostable void CCTK_FNAME(setup_eosmodule) (CCTK_INT* , CCTK_INT* , CCTK_INT*, CCTK_REAL*, CCTK_REAL*, CCTK_REAL*, CCTK_REAL*, CCTK_REAL*); // call HDF5 if we are an IO processor, handle errors from HDF5 #define HDF5_ERROR(fn_call) \ if (doIO) { \ 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); \ } // 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) { DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS CCTK_Info(CCTK_THORNSTRING,"*******************************"); CCTK_Info(CCTK_THORNSTRING,"Reading nuc_eos table file:"); CCTK_Info(CCTK_THORNSTRING,nuceos_table_name); CCTK_Info(CCTK_THORNSTRING,"*******************************"); CCTK_INT my_reader_process = reader_process; if (my_reader_process < 0 || my_reader_process >= CCTK_nProcs(cctkGH)) { CCTK_VWarn(CCTK_WARN_COMPLAIN, __LINE__, __FILE__, CCTK_THORNSTRING, "Requested IO process %d out of range. Reverting to process 0.", my_reader_process); my_reader_process = 0; } const int doIO = !read_table_on_single_process || CCTK_MyProc(cctkGH) == my_reader_process; hid_t file = -1; if (doIO && ! H5Fis_hdf5(nuceos_table_name) > 0) 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_BCAST_EOS_HDF5(NAME,VAR,TYPE,NELEMS) \ do { \ hid_t dataset; \ HDF5_ERROR(dataset = H5Dopen(file, NAME)); \ HDF5_ERROR(H5Dread(dataset, TYPE, H5S_ALL, H5S_ALL, H5P_DEFAULT, VAR)); \ if (read_table_on_single_process) \ BCAST (VAR, sizeof(*(VAR))*(NELEMS)); \ HDF5_ERROR(H5Dclose(dataset)); \ } while (0) // The second reads a given variable into a hyperslab of the alltables array #define READ_BCAST_EOSTABLE_HDF5(NAME,OFF,DIMS) \ do { \ READ_BCAST_EOS_HDF5(NAME,&alltables[(OFF)*(DIMS)[1]],H5T_NATIVE_DOUBLE,(DIMS)[1]); \ } while (0) // Read size of tables CCTK_INT nrho, ntemp, nye; READ_BCAST_EOS_HDF5("pointsrho", &nrho, H5T_NATIVE_INT, 1); READ_BCAST_EOS_HDF5("pointstemp", &ntemp, H5T_NATIVE_INT, 1); READ_BCAST_EOS_HDF5("pointsye", &nye, H5T_NATIVE_INT, 1); // Allocate memory for tables CCTK_REAL *alltables, *logrho, *logtemp, *ye, energy_shift; // protect our use of H5T_NATIVE_DOUBLE and H5T_NATIVE_INT assert (sizeof(CCTK_REAL) == sizeof(double)); assert (sizeof(CCTK_INT) == sizeof(int)); 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}; // Read alltables READ_BCAST_EOSTABLE_HDF5("logpress", 0, table_dims); READ_BCAST_EOSTABLE_HDF5("logenergy", 1, table_dims); READ_BCAST_EOSTABLE_HDF5("entropy", 2, table_dims); READ_BCAST_EOSTABLE_HDF5("munu", 3, table_dims); READ_BCAST_EOSTABLE_HDF5("cs2", 4, table_dims); READ_BCAST_EOSTABLE_HDF5("dedt", 5, table_dims); READ_BCAST_EOSTABLE_HDF5("dpdrhoe", 6, table_dims); READ_BCAST_EOSTABLE_HDF5("dpderho", 7, table_dims); // chemical potentials READ_BCAST_EOSTABLE_HDF5("muhat", 8, table_dims); READ_BCAST_EOSTABLE_HDF5("mu_e", 9, table_dims); READ_BCAST_EOSTABLE_HDF5("mu_p", 10, table_dims); READ_BCAST_EOSTABLE_HDF5("mu_n", 11, table_dims); // compositions READ_BCAST_EOSTABLE_HDF5("Xa", 12, table_dims); READ_BCAST_EOSTABLE_HDF5("Xh", 13, table_dims); READ_BCAST_EOSTABLE_HDF5("Xn", 14, table_dims); READ_BCAST_EOSTABLE_HDF5("Xp", 15, table_dims); // average nucleus READ_BCAST_EOSTABLE_HDF5("Abar", 16, table_dims); READ_BCAST_EOSTABLE_HDF5("Zbar", 17, table_dims); // Gamma READ_BCAST_EOSTABLE_HDF5("gamma", 18, table_dims); // Read additional tables and variables READ_BCAST_EOS_HDF5("logrho", logrho, H5T_NATIVE_DOUBLE, nrho); READ_BCAST_EOS_HDF5("logtemp", logtemp, H5T_NATIVE_DOUBLE, ntemp); READ_BCAST_EOS_HDF5("ye", ye, H5T_NATIVE_DOUBLE, nye); READ_BCAST_EOS_HDF5("energy_shift", &energy_shift, H5T_NATIVE_DOUBLE, 1); HDF5_ERROR(H5Fclose(file)); // Give all values to fortran - which will store pointers to them, so don't // free these arrays. CCTK_FNAME(setup_eosmodule) (&nrho, &ntemp, &nye, alltables, logrho, logtemp, ye, &energy_shift); }