1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
#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
CCTK_Info(CCTK_THORNSTRING,"*******************************");
CCTK_Info(CCTK_THORNSTRING,"Reading nuc_eos table file:");
CCTK_Info(CCTK_THORNSTRING,nuceos_table_name);
CCTK_Info(CCTK_THORNSTRING,"*******************************");
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);
}
|
