Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Hydroinfo_h5.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Hydroinfo_h5.cpp
1 #include<iostream>
2 #include<sstream>
3 #include<fstream>
4 #include<cmath>
5 #include<iomanip>
6 #include<string>
7 #include<stdlib.h>
8 
9 #include "hdf5.h"
10 #include "Hydroinfo_h5.h"
11 
12 using namespace std;
13 
14 HydroinfoH5::HydroinfoH5() {
15  readinFlag = 0;
16  outputFlag = 0;
17 }
18 
19 HydroinfoH5::HydroinfoH5(string filename_in,
20  int bufferSize_in, int Visflag_in) {
21  readHydroinfoH5(filename_in, bufferSize_in, Visflag_in);
22 }
23 
24 HydroinfoH5::HydroinfoH5(int XL_in, int XH_in, double DX_in, int LSX_in,
25  int YL_in, int YH_in, double DY_in, int LSY_in,
26  double Tau0_in, double dTau_in, double LST_in,
27  int Visflag_in, string filename_in) {
28  setHydroFiles(XL_in, XH_in, DX_in, LSX_in, YL_in, YH_in, DY_in, LSY_in,
29  Tau0_in, dTau_in, LST_in, Visflag_in, filename_in);
30 }
31 
32 HydroinfoH5::~HydroinfoH5() {
33  if (readinFlag == 1) {
34  clean_hydro_event();
35  }
36 }
37 
38 void HydroinfoH5::clean_hydro_event() {
39  for (int i=0; i<Buffersize; i++) {
40  for (int j=0; j<dimensionX; j++) {
41  delete[] ed[i][j];
42  delete[] sd[i][j];
43  delete[] vx[i][j];
44  delete[] vy[i][j];
45  delete[] Temperature[i][j];
46  delete[] Pressure[i][j];
47  if (Visflag == 1) {
48  delete[] pi00[i][j];
49  delete[] pi01[i][j];
50  delete[] pi02[i][j];
51  delete[] pi03[i][j];
52  delete[] pi11[i][j];
53  delete[] pi12[i][j];
54  delete[] pi13[i][j];
55  delete[] pi22[i][j];
56  delete[] pi23[i][j];
57  delete[] pi33[i][j];
58  delete[] BulkPi[i][j];
59  }
60  }
61  delete[] ed[i];
62  delete[] sd[i];
63  delete[] vx[i];
64  delete[] vy[i];
65  delete[] Temperature[i];
66  delete[] Pressure[i];
67  if (Visflag == 1) {
68  delete[] pi00[i];
69  delete[] pi01[i];
70  delete[] pi02[i];
71  delete[] pi03[i];
72  delete[] pi11[i];
73  delete[] pi12[i];
74  delete[] pi13[i];
75  delete[] pi22[i];
76  delete[] pi23[i];
77  delete[] pi33[i];
78  delete[] BulkPi[i];
79  }
80  }
81  delete[] ed;
82  delete[] sd;
83  delete[] vx;
84  delete[] vy;
85  delete[] Temperature;
86  delete[] Pressure;
87  if (Visflag == 1) {
88  delete[] pi00;
89  delete[] pi01;
90  delete[] pi02;
91  delete[] pi03;
92  delete[] pi11;
93  delete[] pi12;
94  delete[] pi13;
95  delete[] pi22;
96  delete[] pi23;
97  delete[] pi33;
98  delete[] BulkPi;
99  }
100  readinFlag = 0;
101 }
102 
103 void HydroinfoH5::setHydroFiles(int XL_in, int XH_in, double DX_in, int LSX_in, int YL_in, int YH_in, double DY_in, int LSY_in, double Tau0_in, double dTau_in, double LST_in, int Visflag_in, string filename_in)
104 {
105  outputFlag = 1;
106  grid_XL = XL_in;
107  grid_XH = XH_in;
108  grid_dx = DX_in;
109  grid_YL = YL_in;
110  grid_YH = YH_in;
111  grid_dy = DY_in;
112  grid_Tau0 = Tau0_in;
113  grid_dTau = dTau_in;
114  grid_LSX = LSX_in;
115  grid_LSY = LSY_in;
116  grid_LST = LST_in;
117  LST_cur = grid_LST - 1;
118  Visflag = Visflag_in;
119 
120  int XShift = abs(grid_XL%grid_LSX);
121  int YShift = abs(grid_YL%grid_LSY);
122  dimensionX = (int) (grid_XH - grid_XL - 2*XShift)/grid_LSX + 1;
123  dimensionY = (int) (grid_YH - grid_YL - 2*YShift)/grid_LSY + 1;
124 
125  filename = filename_in;
126  herr_t status;
127  /* Create a new file using default properties. */
128  H5file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
129 
130  /* Create a group named "/Event" in the file. */
131  H5groupEventid = H5Gcreate(H5file_id, "/Event", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
132 
133  writeGroupattribute(H5groupEventid);
134 
135  /* Close the group. */
136  status = H5Gclose(H5groupEventid);
137  /* Terminate access to the file. */
138  status = H5Fclose(H5file_id);
139 }
140 
141 void HydroinfoH5::writeGroupattribute(hid_t H5groupEventid)
142 {
143  int XShift = abs(grid_XL%grid_LSX);
144  int YShift = abs(grid_YL%grid_LSY);
145 
146  addGroupattributeInt(H5groupEventid, "XL", (grid_XL + XShift)/grid_LSX);
147  addGroupattributeInt(H5groupEventid, "XH", (grid_XH - XShift)/grid_LSX);
148  addGroupattributeInt(H5groupEventid, "YL", (grid_YL + YShift)/grid_LSY);
149  addGroupattributeInt(H5groupEventid, "YH", (grid_YH - YShift)/grid_LSY);
150  addGroupattributeDouble(H5groupEventid, "DX", grid_dx*grid_LSX);
151  addGroupattributeDouble(H5groupEventid, "DY", grid_dy*grid_LSY);
152  addGroupattributeDouble(H5groupEventid, "Tau0", grid_Tau0);
153  addGroupattributeDouble(H5groupEventid, "dTau", grid_dTau*grid_LST);
154  addGroupattributeInt(H5groupEventid, "OutputViscousFlag", Visflag);
155 }
156 
157 void HydroinfoH5::addGroupattributeInt(hid_t H5groupEventid, string attName, int attValue)
158 {
159  herr_t status;
160  hsize_t dims;
161  hid_t attribute_id, dataspace_id;
162 
163  /* Create the data space for the attribute. */
164  dims = 1;
165  dataspace_id = H5Screate_simple(1, &dims, NULL);
166  /* Create a dataset attribute. */
167  attribute_id = H5Acreate (H5groupEventid, attName.c_str(), H5T_STD_I32BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
168 
169  /* Write the attribute data. */
170  status = H5Awrite(attribute_id, H5T_NATIVE_INT, &attValue);
171 
172  /* Close the attribute. */
173  status = H5Aclose(attribute_id);
174  /* Close the dataspace. */
175  status = H5Sclose(dataspace_id);
176 }
177 
178 void HydroinfoH5::addGroupattributeDouble(hid_t H5groupEventid, string attName, double attValue)
179 {
180  herr_t status;
181  hsize_t dims;
182  hid_t attribute_id, dataspace_id;
183 
184  /* Create the data space for the attribute. */
185  dims = 1;
186  dataspace_id = H5Screate_simple(1, &dims, NULL);
187  /* Create a dataset attribute. */
188  attribute_id = H5Acreate (H5groupEventid, attName.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
189 
190  /* Write the attribute data. */
191  status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &attValue);
192 
193  /* Close the attribute. */
194  status = H5Aclose(attribute_id);
195  /* Close the dataspace. */
196  status = H5Sclose(dataspace_id);
197 }
198 
199 void HydroinfoH5::writeHydroBlock(int Time_id, double **ed_in, double **sd_in, double **p_in, double **Temp_in, double **Vx_in, double **Vy_in, double **Pi00_in, double **Pi01_in, double **Pi02_in, double **Pi03_in, double **Pi11_in, double **Pi12_in, double **Pi13_in, double **Pi22_in, double **Pi23_in, double **Pi33_in, double ** BulkPi_in)
200 {
201  if(LST_cur != 0)
202  {
203  LST_cur--;
204  return;
205  }
206  else
207  LST_cur = grid_LST - 1;
208 
209  herr_t status;
210  hid_t groupFrameid;
211  hsize_t dim_x = dimensionX;
212  hsize_t dim_y = dimensionY;
213  const hsize_t dims[2] = {dim_x, dim_y};
214 
215  int XShift = abs(grid_XL%grid_LSX);
216  int YShift = abs(grid_YL%grid_LSY);
217 
218  int Frame_id = (int) Time_id/grid_LST;
219  ostringstream frameName;
220  frameName << "Frame_" << setfill('0') << setw(4) << Frame_id;
221 
222  /* Open an existing file. */
223  H5file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
224  H5groupEventid = H5Gopen(H5file_id, "/Event", H5P_DEFAULT);
225 
226  /* Create a group named "/Event/FramName" in the file. */
227  groupFrameid = H5Gcreate(H5groupEventid, frameName.str().c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
228  addGroupattributeDouble(groupFrameid, "Time", grid_Tau0 + grid_dTau*grid_LST*Frame_id);
229 
230  //Dump data into h5 file
231  CSH5dumpBlockdata(groupFrameid, dims, "e", ed_in);
232  CSH5dumpBlockdata(groupFrameid, dims, "s", sd_in);
233  CSH5dumpBlockdata(groupFrameid, dims, "p", p_in);
234  CSH5dumpBlockdata(groupFrameid, dims, "Temp", Temp_in);
235  CSH5dumpBlockdata(groupFrameid, dims, "Vx", Vx_in);
236  CSH5dumpBlockdata(groupFrameid, dims, "Vy", Vy_in);
237  if(Visflag == 1)
238  {
239  CSH5dumpBlockdata(groupFrameid, dims, "Pi00", Pi00_in);
240  CSH5dumpBlockdata(groupFrameid, dims, "Pi01", Pi01_in);
241  CSH5dumpBlockdata(groupFrameid, dims, "Pi02", Pi02_in);
242  CSH5dumpBlockdata(groupFrameid, dims, "Pi03", Pi03_in);
243  CSH5dumpBlockdata(groupFrameid, dims, "Pi11", Pi11_in);
244  CSH5dumpBlockdata(groupFrameid, dims, "Pi12", Pi12_in);
245  CSH5dumpBlockdata(groupFrameid, dims, "Pi13", Pi13_in);
246  CSH5dumpBlockdata(groupFrameid, dims, "Pi22", Pi22_in);
247  CSH5dumpBlockdata(groupFrameid, dims, "Pi23", Pi23_in);
248  CSH5dumpBlockdata(groupFrameid, dims, "Pi33", Pi33_in);
249  CSH5dumpBlockdata(groupFrameid, dims, "BulkPi", BulkPi_in);
250  }
251 
252  status = H5Gclose(H5groupEventid);
253  status = H5Fclose(H5file_id);
254 }
255 
256 void HydroinfoH5::CSH5dumpBlockdata(hid_t group_id, const hsize_t * dims, string DatasetName, double **Dataset)
257 {
258  hid_t dataset_id, dataspace_id;
259  herr_t status;
260 
261  int XShift = abs(grid_XL%grid_LSX);
262  int YShift = abs(grid_YL%grid_LSY);
263 
264  double subDataset[dims[0]][dims[1]];
265  for(int i = 0; i < dims[0]; i++)
266  {
267  int idx_x = XShift + i*grid_LSX;
268  for(int j = 0; j < dims[1]; j++)
269  {
270  int idx_y = YShift + j*grid_LSY;
271  subDataset[i][j] = Dataset[idx_x][idx_y];
272  }
273  }
274 
275  dataspace_id = H5Screate_simple(2, dims, NULL);
276  /* Create the dataset. */
277  dataset_id = H5Dcreate(group_id, DatasetName.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
278 
279  status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, subDataset);
280 
281  /* End access to the dataset and release resources used by it. */
282  status = H5Dclose(dataset_id);
283  /* Terminate access to the data space. */
284  status = H5Sclose(dataspace_id);
285 
286 }
287 
288 void HydroinfoH5::readHydroinfoH5(string filename_in, int bufferSize_in, int Visflag_in)
289 {
290  readinFlag = 1;
291  // flag to determine whether to read evolutions for viscous variables
292  Visflag = Visflag_in;
293 
294  herr_t status;
295  filename = filename_in;
296  const char *fileptr = (char*) filename.c_str();
297  H5file_id = H5Fopen(fileptr, H5F_ACC_RDONLY, H5P_DEFAULT);
298  H5groupEventid = H5Gopen(H5file_id, "/Event", H5P_DEFAULT);
299 
300  readHydrogridInfo();
301  printHydrogridInfo();
302 
303  Buffersize = bufferSize_in;
304  dimensionX = grid_XH - grid_XL + 1;
305  dimensionY = grid_YH - grid_YL + 1;
306 
307  if(Buffersize < grid_Framenum)
308  {
309  cout << "Buffersize is too small, increase it to at lease to " << grid_Framenum << endl;
310  exit(1);
311  }
312  //initialize all matrices
313  ed = new double** [Buffersize];
314  sd = new double** [Buffersize];
315  vx = new double** [Buffersize];
316  vy = new double** [Buffersize];
317  Temperature = new double** [Buffersize];
318  Pressure = new double** [Buffersize];
319  if (Visflag == 1) {
320  pi00 = new double** [Buffersize];
321  pi01 = new double** [Buffersize];
322  pi02 = new double** [Buffersize];
323  pi03 = new double** [Buffersize];
324  pi11 = new double** [Buffersize];
325  pi12 = new double** [Buffersize];
326  pi13 = new double** [Buffersize];
327  pi22 = new double** [Buffersize];
328  pi23 = new double** [Buffersize];
329  pi33 = new double** [Buffersize];
330  BulkPi = new double** [Buffersize];
331  }
332  for(int i=0; i<Buffersize; i++)
333  {
334  ed[i] = new double* [dimensionX];
335  sd[i] = new double* [dimensionX];
336  vx[i] = new double* [dimensionX];
337  vy[i] = new double* [dimensionX];
338  Temperature[i] = new double* [dimensionX];
339  Pressure[i] = new double* [dimensionX];
340  if (Visflag == 1) {
341  pi00[i] = new double* [dimensionX];
342  pi01[i] = new double* [dimensionX];
343  pi02[i] = new double* [dimensionX];
344  pi03[i] = new double* [dimensionX];
345  pi11[i] = new double* [dimensionX];
346  pi12[i] = new double* [dimensionX];
347  pi13[i] = new double* [dimensionX];
348  pi22[i] = new double* [dimensionX];
349  pi23[i] = new double* [dimensionX];
350  pi33[i] = new double* [dimensionX];
351  BulkPi[i] = new double* [dimensionX];
352  }
353  for(int j=0; j<dimensionX; j++)
354  {
355  ed[i][j] = new double [dimensionY];
356  sd[i][j] = new double [dimensionY];
357  vx[i][j] = new double [dimensionY];
358  vy[i][j] = new double [dimensionY];
359  Temperature[i][j] = new double [dimensionY];
360  Pressure[i][j] = new double [dimensionY];
361  if (Visflag == 1) {
362  pi00[i][j] = new double [dimensionY];
363  pi01[i][j] = new double [dimensionY];
364  pi02[i][j] = new double [dimensionY];
365  pi03[i][j] = new double [dimensionY];
366  pi11[i][j] = new double [dimensionY];
367  pi12[i][j] = new double [dimensionY];
368  pi13[i][j] = new double [dimensionY];
369  pi22[i][j] = new double [dimensionY];
370  pi23[i][j] = new double [dimensionY];
371  pi33[i][j] = new double [dimensionY];
372  BulkPi[i][j] = new double [dimensionY];
373  }
374  }
375  }
376 
377  readHydroinfoBuffered_total();
378 
379  status = H5Gclose(H5groupEventid);
380  status = H5Fclose(H5file_id);
381 }
382 
383 
384 void HydroinfoH5::readHydrogridInfo()
385 {
386  herr_t status;
387 
388  grid_XL = readH5Attribute_int(H5groupEventid, "XL");
389  grid_XH = readH5Attribute_int(H5groupEventid, "XH");
390  grid_YL = readH5Attribute_int(H5groupEventid, "YL");
391  grid_YH = readH5Attribute_int(H5groupEventid, "YH");
392  grid_Tau0 = readH5Attribute_double(H5groupEventid, "Tau0");
393  grid_dTau = readH5Attribute_double(H5groupEventid, "dTau");
394  grid_dx = readH5Attribute_double(H5groupEventid, "DX");
395  grid_dy = readH5Attribute_double(H5groupEventid, "DY");
396 
397  grid_X0 = grid_XL * grid_dx;
398  grid_Y0 = grid_YL * grid_dy;
399  grid_Xmax = grid_XH * grid_dx;
400  grid_Ymax = grid_YH * grid_dy;
401 
402  hsize_t tempFramenum;
403  status = H5Gget_num_objs(H5groupEventid, &tempFramenum);
404  grid_Framenum = (int) tempFramenum;
405  grid_Taumax = grid_Tau0 + (grid_Framenum - 1)*grid_dTau;
406 
407  int tempflag = readH5Attribute_int(H5groupEventid, "OutputViscousFlag");
408  Visflag = Visflag*tempflag;
409 }
410 
411 void HydroinfoH5::printHydrogridInfo()
412 {
413  cout << "-----------------------------------------" << endl;
414  cout << "-----------hydro grid info---------------" << endl;
415  cout << "-----------------------------------------" << endl;
416  cout << "XL = " << grid_XL << endl;
417  cout << "XH = " << grid_XH << endl;
418  cout << "DX = " << grid_dx << " fm" << endl;
419  cout << "YL = " << grid_YL << endl;
420  cout << "YH = " << grid_YH << endl;
421  cout << "DY = " << grid_dy << " fm" << endl;
422  cout << "Tau0 = " << grid_Tau0 << " fm" << endl;
423  cout << "dTau = " << grid_dTau << " fm" << endl;
424  cout << "Number of Frames: " << grid_Framenum << endl;
425  cout << "Taumax = " << grid_Taumax << " fm" << endl;
426  cout << "Read in viscous information? ";
427  if(Visflag == 1)
428  cout << " Yes!" << endl;
429  else
430  cout << " No!" << endl;
431  cout << "-----------------------------------------" << endl;
432 }
433 
434 int HydroinfoH5::readH5Attribute_int(hid_t id, string attributeName)
435 {
436  int attributeValue;
437  hid_t attr;
438  herr_t ret;
439  attr = H5Aopen_name(id, attributeName.c_str());
440  ret = H5Aread(attr, H5T_NATIVE_INT, &attributeValue);
441  ret = H5Aclose(attr);
442  return(attributeValue);
443 }
444 
445 double HydroinfoH5::readH5Attribute_double(hid_t id, string attributeName)
446 {
447  double attributeValue;
448  hid_t attr;
449  herr_t ret;
450  attr = H5Aopen_name(id, attributeName.c_str());
451  ret = H5Aread(attr, H5T_NATIVE_DOUBLE, &attributeValue);
452  ret = H5Aclose(attr);
453  return(attributeValue);
454 }
455 
456 void HydroinfoH5::readHydroinfoBuffered_total()
457 {
458  hid_t group_id;
459  herr_t status;
460 
461  int frameIdx;
462  for(int i=0; i<Buffersize; i++)
463  {
464  frameIdx = i;
465  if(frameIdx < (int) grid_Framenum)
466  {
467  stringstream frameName;
468  frameName << "Frame_" << setw(4) << setfill('0') << frameIdx;
469  group_id = H5Gopen(H5groupEventid, frameName.str().c_str(), H5P_DEFAULT);
470 
471  readH5Dataset_double(group_id, "e", ed[i]);
472  readH5Dataset_double(group_id, "s", sd[i]);
473  readH5Dataset_double(group_id, "Vx", vx[i]);
474  readH5Dataset_double(group_id, "Vy", vy[i]);
475  readH5Dataset_double(group_id, "Temp", Temperature[i]);
476  readH5Dataset_double(group_id, "P", Pressure[i]);
477  if(Visflag == 1)
478  {
479  readH5Dataset_double(group_id, "Pi00", pi00[i]);
480  readH5Dataset_double(group_id, "Pi01", pi01[i]);
481  readH5Dataset_double(group_id, "Pi02", pi02[i]);
482  readH5Dataset_double(group_id, "Pi03", pi03[i]);
483  readH5Dataset_double(group_id, "Pi11", pi11[i]);
484  readH5Dataset_double(group_id, "Pi12", pi12[i]);
485  readH5Dataset_double(group_id, "Pi13", pi13[i]);
486  readH5Dataset_double(group_id, "Pi22", pi22[i]);
487  readH5Dataset_double(group_id, "Pi23", pi23[i]);
488  readH5Dataset_double(group_id, "Pi33", pi33[i]);
489  readH5Dataset_double(group_id, "BulkPi", BulkPi[i]);
490  }
491  status = H5Gclose(group_id);
492  }
493  else
494  break;
495  }
496 }
497 
498 void HydroinfoH5::readHydroinfoSingleframe(int frameIdx)
499 {
500  hid_t group_id;
501  herr_t status;
502 
503  if(frameIdx < (int) grid_Framenum)
504  {
505  stringstream frameName;
506  frameName << "Frame_" << setw(4) << setfill('0') << frameIdx;
507  group_id = H5Gopen(H5groupEventid, frameName.str().c_str(), H5P_DEFAULT);
508 
509  int Idx = frameIdx;
510 
511  readH5Dataset_double(group_id, "e", ed[Idx]);
512  readH5Dataset_double(group_id, "s", sd[Idx]);
513  readH5Dataset_double(group_id, "Vx", vx[Idx]);
514  readH5Dataset_double(group_id, "Vy", vy[Idx]);
515  readH5Dataset_double(group_id, "Temp", Temperature[Idx]);
516  readH5Dataset_double(group_id, "P", Pressure[Idx]);
517  if(Visflag == 1)
518  {
519  readH5Dataset_double(group_id, "Pi00", pi00[Idx]);
520  readH5Dataset_double(group_id, "Pi01", pi01[Idx]);
521  readH5Dataset_double(group_id, "Pi02", pi02[Idx]);
522  readH5Dataset_double(group_id, "Pi03", pi03[Idx]);
523  readH5Dataset_double(group_id, "Pi11", pi11[Idx]);
524  readH5Dataset_double(group_id, "Pi12", pi12[Idx]);
525  readH5Dataset_double(group_id, "Pi13", pi13[Idx]);
526  readH5Dataset_double(group_id, "Pi22", pi22[Idx]);
527  readH5Dataset_double(group_id, "Pi23", pi23[Idx]);
528  readH5Dataset_double(group_id, "Pi33", pi33[Idx]);
529  readH5Dataset_double(group_id, "BulkPi", BulkPi[Idx]);
530  }
531  status = H5Gclose(group_id);
532  }
533  else
534  {
535  cout << "Error: readHydroinfoSingleframe :: frameIdx exceed maximum frame number from hydro" << endl;
536  cout << "frameIdx = " << frameIdx << endl;
537  exit(1);
538  }
539 }
540 
541 void HydroinfoH5::readH5Dataset_double(hid_t id, string datasetName, double** dset_data)
542 {
543  herr_t status;
544  hid_t dataset_id;
545 
546  double temp_data[dimensionX][dimensionY];
547  dataset_id = H5Dopen(id, datasetName.c_str(), H5P_DEFAULT);
548  status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data);
549  for(int i=0; i<dimensionX; i++)
550  for(int j=0; j<dimensionY; j++)
551  dset_data[i][j] = temp_data[i][j];
552  status = H5Dclose(dataset_id);
553 }
554 
555 void HydroinfoH5::getHydroinfoOnlattice(int frameIdx, int xIdx, int yIdx, hydrofluidCell* fluidCellptr)
556 {
557  if(frameIdx < 0 || frameIdx > grid_Framenum || xIdx < 0 || xIdx > (grid_XH - grid_XL) || yIdx < 0 || yIdx > (grid_YH - grid_YL))
558  {
559  cout << "Error: getHydroinfoOnlattice:: Index is wrong" << endl;
560  cout << "frameIdx = " << frameIdx << " xIdx = " << xIdx
561  << "yIdx = " << yIdx << endl;
562  exit(1);
563  }
564  fluidCellptr->ed = ed[frameIdx][xIdx][yIdx];
565  fluidCellptr->sd = sd[frameIdx][xIdx][yIdx];
566  fluidCellptr->vx = vx[frameIdx][xIdx][yIdx];
567  fluidCellptr->vy = vy[frameIdx][xIdx][yIdx];
568  fluidCellptr->temperature = Temperature[frameIdx][xIdx][yIdx];
569  fluidCellptr->pressure = Pressure[frameIdx][xIdx][yIdx];
570  fluidCellptr->pi[0][0] = pi00[frameIdx][xIdx][yIdx];
571  fluidCellptr->pi[0][1] = pi01[frameIdx][xIdx][yIdx];
572  fluidCellptr->pi[0][2] = pi02[frameIdx][xIdx][yIdx];
573  fluidCellptr->pi[0][3] = pi03[frameIdx][xIdx][yIdx];
574  fluidCellptr->pi[1][0] = fluidCellptr->pi[0][1];
575  fluidCellptr->pi[1][1] = pi11[frameIdx][xIdx][yIdx];
576  fluidCellptr->pi[1][2] = pi12[frameIdx][xIdx][yIdx];
577  fluidCellptr->pi[1][3] = pi13[frameIdx][xIdx][yIdx];
578  fluidCellptr->pi[2][0] = fluidCellptr->pi[0][2];
579  fluidCellptr->pi[2][1] = fluidCellptr->pi[1][2];
580  fluidCellptr->pi[2][2] = pi22[frameIdx][xIdx][yIdx];
581  fluidCellptr->pi[2][3] = pi23[frameIdx][xIdx][yIdx];
582  fluidCellptr->pi[3][0] = fluidCellptr->pi[0][3];
583  fluidCellptr->pi[3][1] = fluidCellptr->pi[1][3];
584  fluidCellptr->pi[3][2] = fluidCellptr->pi[2][3];
585  fluidCellptr->pi[3][3] = pi33[frameIdx][xIdx][yIdx];
586  fluidCellptr->bulkPi = BulkPi[frameIdx][xIdx][yIdx];
587 }
588 
589 
590 void HydroinfoH5::getHydroinfo(double tau, double x, double y, hydrofluidCell* fluidCellptr)
591 {
592  double eps = 1e-10;
593  if(tau < grid_Tau0 || tau > grid_Taumax-eps || x < grid_X0 || x > grid_Xmax-eps || y < grid_Y0 || y > grid_Ymax-eps)
594  {
595  setZero_fluidCell(fluidCellptr);
596  return;
597  }
598  int frameIdx, xIdx, yIdx;
599  double tauInc, xInc, yInc;
600  double temp;
601 
602  temp = (tau - grid_Tau0)/grid_dTau;
603  frameIdx = (int) floor(temp);
604  tauInc = temp - frameIdx;
605 
606  temp = (x - grid_X0)/grid_dx;
607  xIdx = (int) floor(temp);
608  xInc = temp - xIdx;
609 
610  temp = (y - grid_Y0)/grid_dy;
611  yIdx = (int) floor(temp);
612  yInc = temp - yIdx;
613 
614  fluidCellptr->ed = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, ed);
615  fluidCellptr->sd = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, sd);
616  fluidCellptr->vx = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, vx);
617  fluidCellptr->vy = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, vy);
618  fluidCellptr->temperature = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, Temperature);
619  fluidCellptr->pressure = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, Pressure);
620  if(Visflag == 1)
621  {
622  fluidCellptr->pi[0][0] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi00);
623  fluidCellptr->pi[0][1] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi01);
624  fluidCellptr->pi[0][2] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi02);
625  fluidCellptr->pi[0][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi03);
626  fluidCellptr->pi[1][0] = fluidCellptr->pi[0][1];
627  fluidCellptr->pi[1][1] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi11);
628  fluidCellptr->pi[1][2] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi12);
629  fluidCellptr->pi[1][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi13);
630  fluidCellptr->pi[2][0] = fluidCellptr->pi[0][2];
631  fluidCellptr->pi[2][1] = fluidCellptr->pi[1][2];
632  fluidCellptr->pi[2][2] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi22);
633  fluidCellptr->pi[2][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi23);
634  fluidCellptr->pi[3][0] = fluidCellptr->pi[0][3];
635  fluidCellptr->pi[3][1] = fluidCellptr->pi[1][3];
636  fluidCellptr->pi[3][2] = fluidCellptr->pi[2][3];
637  fluidCellptr->pi[3][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi33);
638  fluidCellptr->bulkPi = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, BulkPi);
639  }
640  else
641  {
642  for(int i=0; i<4; i++)
643  for(int j=0; j<4; j++)
644  fluidCellptr->pi[i][j] = 0.0e0;
645  fluidCellptr->bulkPi = 0.0e0;
646  }
647 }
648 
649 void HydroinfoH5::setZero_fluidCell(hydrofluidCell* fluidCellptr)
650 {
651  fluidCellptr->ed = 0.0e0;
652  fluidCellptr->sd = 0.0e0;
653  fluidCellptr->vx = 0.0e0;
654  fluidCellptr->vy = 0.0e0;
655  fluidCellptr->temperature = 0.0e0;
656  fluidCellptr->pressure = 0.0e0;
657  for(int i=0; i<4; i++)
658  for(int j=0; j<4; j++)
659  fluidCellptr->pi[i][j] = 0.0e0;
660  fluidCellptr->bulkPi = 0.0e0;
661 }
662 
663 double HydroinfoH5::cubeInterpShell(int idx_x, int idx_y, int idx_z, double x, double y, double z, double ***dataset)
664 {
665  double result = cubeInterp(x, y, z,
666  dataset[idx_z][idx_x][idx_y], dataset[idx_z][idx_x+1][idx_y], dataset[idx_z][idx_x][idx_y+1], dataset[idx_z][idx_x+1][idx_y+1],
667  dataset[idx_z+1][idx_x][idx_y], dataset[idx_z+1][idx_x+1][idx_y], dataset[idx_z+1][idx_x][idx_y+1], dataset[idx_z+1][idx_x+1][idx_y+1]);
668  return(result);
669 }
670 
671 double HydroinfoH5::cubeInterp(double x, double y, double z, double A000, double A100, double A010, double A110, double A001, double A101, double A011, double A111)
672 // Perform a 3d interpolation. The known data are A### located at the 8 corners,
673 // labels using the xyz order. Therefore A000 is value at the origin and A010
674 // is the value at (x=0,y=1,z=0). Note that the coordinate (x,y,z) must be
675 // constrained to the unit cube. Axyz is the return value.
676 {
677  /* for debug
678  cout << A000 << " " << A100 << " " << A010 << " " << A110 << endl;
679  cout << A001 << " " << A101 << " " << A011 << " " << A111 << endl; */
680 
681  double Axyz = A000*(1-x)*(1-y)*(1-z) + A100*x*(1-y)*(1-z)
682  + A010*(1-x)*y*(1-z) + A001*(1-x)*(1-y)*z
683  + A101*x*(1-y)*z + A011*(1-x)*y*z + A110*x*y*(1-z)
684  + A111*x*y*z;
685  return(Axyz);
686 }