external_packages/hydro_from_external_file/Hydroinfo_h5.cpp

0001 #include<iostream>
0002 #include<sstream>
0003 #include<fstream>
0004 #include<cmath>
0005 #include<iomanip>
0006 #include<string>
0007 #include<stdlib.h>
0008
0009 #include "hdf5.h"
0010 #include "Hydroinfo_h5.h"
0011
0012 using namespace std;
0013
0014 HydroinfoH5::HydroinfoH5() {
0015    readinFlag = 0;
0016    outputFlag = 0;
0017 }
0018
0019 HydroinfoH5::HydroinfoH5(string filename_in,
0020                          int bufferSize_in, int Visflag_in) {
0021    readHydroinfoH5(filename_in, bufferSize_in, Visflag_in);
0022 }
0023
0024 HydroinfoH5::HydroinfoH5(int XL_in, int XH_in, double DX_in, int LSX_in,
0025                          int YL_in, int YH_in, double DY_in, int LSY_in,
0026                          double Tau0_in, double dTau_in, double LST_in,
0027                          int Visflag_in, string filename_in) {
0028    setHydroFiles(XL_in, XH_in, DX_in, LSX_in, YL_in, YH_in, DY_in, LSY_in,
0029                  Tau0_in, dTau_in, LST_in, Visflag_in, filename_in);
0030 }
0031
0032 HydroinfoH5::~HydroinfoH5() {
0033    if (readinFlag == 1) {
0034        clean_hydro_event();
0035    }
0036 }
0037
0038 void HydroinfoH5::clean_hydro_event() {
0039     for (int i=0; i<Buffersize; i++) {
0040         for (int j=0; j<dimensionX; j++) {
0041            delete[] ed[i][j];
0042            delete[] sd[i][j];
0043            delete[] vx[i][j];
0044            delete[] vy[i][j];
0045            delete[] Temperature[i][j];
0046            delete[] Pressure[i][j];
0047            if (Visflag == 1) {
0048               delete[] pi00[i][j];
0049               delete[] pi01[i][j];
0050               delete[] pi02[i][j];
0051               delete[] pi03[i][j];
0052               delete[] pi11[i][j];
0053               delete[] pi12[i][j];
0054               delete[] pi13[i][j];
0055               delete[] pi22[i][j];
0056               delete[] pi23[i][j];
0057               delete[] pi33[i][j];
0058               delete[] BulkPi[i][j];
0059            }
0060         }
0061         delete[] ed[i];
0062         delete[] sd[i];
0063         delete[] vx[i];
0064         delete[] vy[i];
0065         delete[] Temperature[i];
0066         delete[] Pressure[i];
0067         if (Visflag == 1) {
0068            delete[] pi00[i];
0069            delete[] pi01[i];
0070            delete[] pi02[i];
0071            delete[] pi03[i];
0072            delete[] pi11[i];
0073            delete[] pi12[i];
0074            delete[] pi13[i];
0075            delete[] pi22[i];
0076            delete[] pi23[i];
0077            delete[] pi33[i];
0078            delete[] BulkPi[i];
0079         }
0080     }
0081     delete[] ed;
0082     delete[] sd;
0083     delete[] vx;
0084     delete[] vy;
0085     delete[] Temperature;
0086     delete[] Pressure;
0087     if (Visflag == 1) {
0088        delete[] pi00;
0089        delete[] pi01;
0090        delete[] pi02;
0091        delete[] pi03;
0092        delete[] pi11;
0093        delete[] pi12;
0094        delete[] pi13;
0095        delete[] pi22;
0096        delete[] pi23;
0097        delete[] pi33;
0098        delete[] BulkPi;
0099     }
0100     readinFlag = 0;
0101 }
0102
0103 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)
0104 {
0105     outputFlag = 1;
0106     grid_XL = XL_in;
0107     grid_XH = XH_in;
0108     grid_dx = DX_in;
0109     grid_YL = YL_in;
0110     grid_YH = YH_in;
0111     grid_dy = DY_in;
0112     grid_Tau0 = Tau0_in;
0113     grid_dTau = dTau_in;
0114     grid_LSX = LSX_in;
0115     grid_LSY = LSY_in;
0116     grid_LST = LST_in;
0117     LST_cur = grid_LST - 1;
0118     Visflag = Visflag_in;
0119
0120     int XShift = abs(grid_XL%grid_LSX);
0121     int YShift = abs(grid_YL%grid_LSY);
0122     dimensionX = (int) (grid_XH - grid_XL - 2*XShift)/grid_LSX + 1;
0123     dimensionY = (int) (grid_YH - grid_YL - 2*YShift)/grid_LSY + 1;
0124
0125     filename = filename_in;
0126     herr_t status;
0127     /* Create a new file using default properties. */
0128     H5file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
0129
0130     /* Create a group named "/Event" in the file. */
0131     H5groupEventid = H5Gcreate(H5file_id, "/Event", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
0132
0133     writeGroupattribute(H5groupEventid);
0134
0135     /* Close the group. */
0136     status = H5Gclose(H5groupEventid);
0137     /* Terminate access to the file. */
0138     status = H5Fclose(H5file_id);
0139 }
0140
0141 void HydroinfoH5::writeGroupattribute(hid_t H5groupEventid)
0142 {
0143     int XShift = abs(grid_XL%grid_LSX);
0144     int YShift = abs(grid_YL%grid_LSY);
0145
0146     addGroupattributeInt(H5groupEventid, "XL", (grid_XL + XShift)/grid_LSX);
0147     addGroupattributeInt(H5groupEventid, "XH", (grid_XH - XShift)/grid_LSX);
0148     addGroupattributeInt(H5groupEventid, "YL", (grid_YL + YShift)/grid_LSY);
0149     addGroupattributeInt(H5groupEventid, "YH", (grid_YH - YShift)/grid_LSY);
0150     addGroupattributeDouble(H5groupEventid, "DX", grid_dx*grid_LSX);
0151     addGroupattributeDouble(H5groupEventid, "DY", grid_dy*grid_LSY);
0152     addGroupattributeDouble(H5groupEventid, "Tau0", grid_Tau0);
0153     addGroupattributeDouble(H5groupEventid, "dTau", grid_dTau*grid_LST);
0154     addGroupattributeInt(H5groupEventid, "OutputViscousFlag", Visflag);
0155 }
0156
0157 void HydroinfoH5::addGroupattributeInt(hid_t H5groupEventid, string attName, int attValue)
0158 {
0159    herr_t status;
0160    hsize_t dims;
0161    hid_t attribute_id, dataspace_id;
0162
0163    /* Create the data space for the attribute. */
0164    dims = 1;
0165    dataspace_id = H5Screate_simple(1, &dims, NULL);
0166    /* Create a dataset attribute. */
0167    attribute_id = H5Acreate (H5groupEventid, attName.c_str(), H5T_STD_I32BE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
0168
0169    /* Write the attribute data. */
0170    status = H5Awrite(attribute_id, H5T_NATIVE_INT, &attValue);
0171
0172    /* Close the attribute. */
0173    status = H5Aclose(attribute_id);
0174    /* Close the dataspace. */
0175    status = H5Sclose(dataspace_id);
0176 }
0177
0178 void HydroinfoH5::addGroupattributeDouble(hid_t H5groupEventid, string attName, double attValue)
0179 {
0180    herr_t status;
0181    hsize_t dims;
0182    hid_t attribute_id, dataspace_id;
0183
0184    /* Create the data space for the attribute. */
0185    dims = 1;
0186    dataspace_id = H5Screate_simple(1, &dims, NULL);
0187    /* Create a dataset attribute. */
0188    attribute_id = H5Acreate (H5groupEventid, attName.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
0189
0190    /* Write the attribute data. */
0191    status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &attValue);
0192
0193    /* Close the attribute. */
0194    status = H5Aclose(attribute_id);
0195    /* Close the dataspace. */
0196    status = H5Sclose(dataspace_id);
0197 }
0198
0199 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)
0200 {
0201    if(LST_cur != 0)
0202    {
0203       LST_cur--;
0204       return;
0205    }
0206    else
0207       LST_cur = grid_LST - 1;
0208
0209    herr_t status;
0210    hid_t groupFrameid;
0211    hsize_t dim_x = dimensionX;
0212    hsize_t dim_y = dimensionY;
0213    const hsize_t dims[2] = {dim_x, dim_y};
0214
0215    int XShift = abs(grid_XL%grid_LSX);
0216    int YShift = abs(grid_YL%grid_LSY);
0217
0218    int Frame_id = (int) Time_id/grid_LST;
0219    ostringstream frameName;
0220    frameName << "Frame_" << setfill('0') << setw(4) << Frame_id;
0221
0222    /* Open an existing file. */
0223    H5file_id = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
0224    H5groupEventid = H5Gopen(H5file_id, "/Event", H5P_DEFAULT);
0225
0226     /* Create a group named "/Event/FramName" in the file. */
0227     groupFrameid = H5Gcreate(H5groupEventid, frameName.str().c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
0228     addGroupattributeDouble(groupFrameid, "Time", grid_Tau0 + grid_dTau*grid_LST*Frame_id);
0229
0230     //Dump data into h5 file
0231     CSH5dumpBlockdata(groupFrameid, dims, "e", ed_in);
0232     CSH5dumpBlockdata(groupFrameid, dims, "s", sd_in);
0233     CSH5dumpBlockdata(groupFrameid, dims, "p", p_in);
0234     CSH5dumpBlockdata(groupFrameid, dims, "Temp", Temp_in);
0235     CSH5dumpBlockdata(groupFrameid, dims, "Vx", Vx_in);
0236     CSH5dumpBlockdata(groupFrameid, dims, "Vy", Vy_in);
0237     if(Visflag == 1)
0238     {
0239        CSH5dumpBlockdata(groupFrameid, dims, "Pi00", Pi00_in);
0240        CSH5dumpBlockdata(groupFrameid, dims, "Pi01", Pi01_in);
0241        CSH5dumpBlockdata(groupFrameid, dims, "Pi02", Pi02_in);
0242        CSH5dumpBlockdata(groupFrameid, dims, "Pi03", Pi03_in);
0243        CSH5dumpBlockdata(groupFrameid, dims, "Pi11", Pi11_in);
0244        CSH5dumpBlockdata(groupFrameid, dims, "Pi12", Pi12_in);
0245        CSH5dumpBlockdata(groupFrameid, dims, "Pi13", Pi13_in);
0246        CSH5dumpBlockdata(groupFrameid, dims, "Pi22", Pi22_in);
0247        CSH5dumpBlockdata(groupFrameid, dims, "Pi23", Pi23_in);
0248        CSH5dumpBlockdata(groupFrameid, dims, "Pi33", Pi33_in);
0249        CSH5dumpBlockdata(groupFrameid, dims, "BulkPi", BulkPi_in);
0250     }
0251
0252     status = H5Gclose(H5groupEventid);
0253     status = H5Fclose(H5file_id);
0254 }
0255
0256 void HydroinfoH5::CSH5dumpBlockdata(hid_t group_id, const hsize_t * dims, string DatasetName, double **Dataset)
0257 {
0258    hid_t dataset_id, dataspace_id;
0259    herr_t status;
0260
0261    int XShift = abs(grid_XL%grid_LSX);
0262    int YShift = abs(grid_YL%grid_LSY);
0263
0264    double subDataset[dims[0]][dims[1]];
0265    for(int i = 0; i < dims[0]; i++)
0266    {
0267       int idx_x = XShift + i*grid_LSX;
0268       for(int j = 0; j < dims[1]; j++)
0269       {
0270          int idx_y = YShift + j*grid_LSY;
0271          subDataset[i][j] = Dataset[idx_x][idx_y];
0272       }
0273    }
0274
0275    dataspace_id = H5Screate_simple(2, dims, NULL);
0276    /* Create the dataset. */
0277    dataset_id = H5Dcreate(group_id, DatasetName.c_str(), H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
0278
0279    status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, subDataset);
0280
0281    /* End access to the dataset and release resources used by it. */
0282    status = H5Dclose(dataset_id);
0283    /* Terminate access to the data space. */
0284    status = H5Sclose(dataspace_id);
0285
0286 }
0287
0288 void HydroinfoH5::readHydroinfoH5(string filename_in, int bufferSize_in, int Visflag_in)
0289 {
0290    readinFlag = 1;
0291    // flag to determine whether to read evolutions for viscous variables
0292    Visflag = Visflag_in;
0293
0294    herr_t status;
0295    filename = filename_in;
0296    const char *fileptr = (char*) filename.c_str();
0297    H5file_id = H5Fopen(fileptr, H5F_ACC_RDONLY, H5P_DEFAULT);  ////H5F_ACC_RDWR, H5F_ACC_RDONLY
0298    H5groupEventid = H5Gopen(H5file_id, "/Event", H5P_DEFAULT);
0299
0300    readHydrogridInfo();
0301    printHydrogridInfo();
0302
0303    Buffersize = bufferSize_in;
0304    dimensionX = grid_XH - grid_XL + 1;
0305    dimensionY = grid_YH - grid_YL + 1;
0306
0307    if(Buffersize < grid_Framenum)
0308    {
0309       cout << "Buffersize is too small, increase it to at lease to " << grid_Framenum << endl;
0310       exit(1);
0311    }
0312    //initialize all matrices
0313    ed = new double** [Buffersize];
0314    sd = new double** [Buffersize];
0315    vx = new double** [Buffersize];
0316    vy = new double** [Buffersize];
0317    Temperature = new double** [Buffersize];
0318    Pressure = new double** [Buffersize];
0319    if (Visflag == 1) {
0320       pi00 = new double** [Buffersize];
0321       pi01 = new double** [Buffersize];
0322       pi02 = new double** [Buffersize];
0323       pi03 = new double** [Buffersize];
0324       pi11 = new double** [Buffersize];
0325       pi12 = new double** [Buffersize];
0326       pi13 = new double** [Buffersize];
0327       pi22 = new double** [Buffersize];
0328       pi23 = new double** [Buffersize];
0329       pi33 = new double** [Buffersize];
0330       BulkPi = new double** [Buffersize];
0331    }
0332    for(int i=0; i<Buffersize; i++)
0333    {
0334       ed[i] = new double* [dimensionX];
0335       sd[i] = new double* [dimensionX];
0336       vx[i] = new double* [dimensionX];
0337       vy[i] = new double* [dimensionX];
0338       Temperature[i] = new double* [dimensionX];
0339       Pressure[i] = new double* [dimensionX];
0340       if (Visflag == 1) {
0341          pi00[i] = new double* [dimensionX];
0342          pi01[i] = new double* [dimensionX];
0343          pi02[i] = new double* [dimensionX];
0344          pi03[i] = new double* [dimensionX];
0345          pi11[i] = new double* [dimensionX];
0346          pi12[i] = new double* [dimensionX];
0347          pi13[i] = new double* [dimensionX];
0348          pi22[i] = new double* [dimensionX];
0349          pi23[i] = new double* [dimensionX];
0350          pi33[i] = new double* [dimensionX];
0351          BulkPi[i] = new double* [dimensionX];
0352       }
0353       for(int j=0; j<dimensionX; j++)
0354       {
0355          ed[i][j] = new double [dimensionY];
0356          sd[i][j] = new double [dimensionY];
0357          vx[i][j] = new double [dimensionY];
0358          vy[i][j] = new double [dimensionY];
0359          Temperature[i][j] = new double [dimensionY];
0360          Pressure[i][j] = new double [dimensionY];
0361          if (Visflag == 1) {
0362             pi00[i][j] = new double [dimensionY];
0363             pi01[i][j] = new double [dimensionY];
0364             pi02[i][j] = new double [dimensionY];
0365             pi03[i][j] = new double [dimensionY];
0366             pi11[i][j] = new double [dimensionY];
0367             pi12[i][j] = new double [dimensionY];
0368             pi13[i][j] = new double [dimensionY];
0369             pi22[i][j] = new double [dimensionY];
0370             pi23[i][j] = new double [dimensionY];
0371             pi33[i][j] = new double [dimensionY];
0372             BulkPi[i][j] = new double [dimensionY];
0373          }
0374       }
0375    }
0376
0377    readHydroinfoBuffered_total();
0378
0379    status = H5Gclose(H5groupEventid);
0380    status = H5Fclose(H5file_id);
0381 }
0382
0383
0384 void HydroinfoH5::readHydrogridInfo()
0385 {
0386    herr_t status;
0387
0388    grid_XL = readH5Attribute_int(H5groupEventid, "XL");
0389    grid_XH = readH5Attribute_int(H5groupEventid, "XH");
0390    grid_YL = readH5Attribute_int(H5groupEventid, "YL");
0391    grid_YH = readH5Attribute_int(H5groupEventid, "YH");
0392    grid_Tau0 = readH5Attribute_double(H5groupEventid, "Tau0");
0393    grid_dTau = readH5Attribute_double(H5groupEventid, "dTau");
0394    grid_dx = readH5Attribute_double(H5groupEventid, "DX");
0395    grid_dy = readH5Attribute_double(H5groupEventid, "DY");
0396
0397    grid_X0 = grid_XL * grid_dx;
0398    grid_Y0 = grid_YL * grid_dy;
0399    grid_Xmax = grid_XH * grid_dx;
0400    grid_Ymax = grid_YH * grid_dy;
0401
0402    hsize_t tempFramenum;
0403    status = H5Gget_num_objs(H5groupEventid, &tempFramenum);
0404    grid_Framenum = (int) tempFramenum;
0405    grid_Taumax = grid_Tau0 + (grid_Framenum - 1)*grid_dTau;
0406
0407    int tempflag = readH5Attribute_int(H5groupEventid, "OutputViscousFlag");
0408    Visflag = Visflag*tempflag;
0409 }
0410
0411 void HydroinfoH5::printHydrogridInfo()
0412 {
0413    cout << "-----------------------------------------" << endl;
0414    cout << "-----------hydro grid info---------------" << endl;
0415    cout << "-----------------------------------------" << endl;
0416    cout << "XL = " << grid_XL << endl;
0417    cout << "XH = " << grid_XH << endl;
0418    cout << "DX = " << grid_dx << " fm" << endl;
0419    cout << "YL = " << grid_YL << endl;
0420    cout << "YH = " << grid_YH << endl;
0421    cout << "DY = " << grid_dy << " fm" << endl;
0422    cout << "Tau0 = " << grid_Tau0 << " fm" << endl;
0423    cout << "dTau = " << grid_dTau << " fm" << endl;
0424    cout << "Number of Frames: " << grid_Framenum << endl;
0425    cout << "Taumax = " << grid_Taumax << " fm" << endl;
0426    cout << "Read in viscous information? ";
0427    if(Visflag == 1)
0428       cout << " Yes!" << endl;
0429    else
0430       cout << " No!" << endl;
0431    cout << "-----------------------------------------" << endl;
0432 }
0433
0434 int HydroinfoH5::readH5Attribute_int(hid_t id, string attributeName)
0435 {
0436    int attributeValue;
0437    hid_t attr;
0438    herr_t ret;
0439    attr = H5Aopen_name(id, attributeName.c_str());
0440    ret  = H5Aread(attr, H5T_NATIVE_INT, &attributeValue);
0441    ret =  H5Aclose(attr);
0442    return(attributeValue);
0443 }
0444
0445 double HydroinfoH5::readH5Attribute_double(hid_t id, string attributeName)
0446 {
0447    double attributeValue;
0448    hid_t attr;
0449    herr_t ret;
0450    attr = H5Aopen_name(id, attributeName.c_str());
0451    ret  = H5Aread(attr, H5T_NATIVE_DOUBLE, &attributeValue);
0452    ret =  H5Aclose(attr);
0453    return(attributeValue);
0454 }
0455
0456 void HydroinfoH5::readHydroinfoBuffered_total()
0457 {
0458    hid_t group_id;
0459    herr_t status;
0460
0461    int frameIdx;
0462    for(int i=0; i<Buffersize; i++)
0463    {
0464       frameIdx = i;
0465       if(frameIdx < (int) grid_Framenum)
0466       {
0467          stringstream frameName;
0468          frameName << "Frame_" <<  setw(4) << setfill('0') << frameIdx;
0469          group_id = H5Gopen(H5groupEventid, frameName.str().c_str(), H5P_DEFAULT);
0470
0471          readH5Dataset_double(group_id, "e", ed[i]);
0472          readH5Dataset_double(group_id, "s", sd[i]);
0473          readH5Dataset_double(group_id, "Vx", vx[i]);
0474          readH5Dataset_double(group_id, "Vy", vy[i]);
0475          readH5Dataset_double(group_id, "Temp", Temperature[i]);
0476          readH5Dataset_double(group_id, "P", Pressure[i]);
0477          if(Visflag == 1)
0478          {
0479             readH5Dataset_double(group_id, "Pi00", pi00[i]);
0480             readH5Dataset_double(group_id, "Pi01", pi01[i]);
0481             readH5Dataset_double(group_id, "Pi02", pi02[i]);
0482             readH5Dataset_double(group_id, "Pi03", pi03[i]);
0483             readH5Dataset_double(group_id, "Pi11", pi11[i]);
0484             readH5Dataset_double(group_id, "Pi12", pi12[i]);
0485             readH5Dataset_double(group_id, "Pi13", pi13[i]);
0486             readH5Dataset_double(group_id, "Pi22", pi22[i]);
0487             readH5Dataset_double(group_id, "Pi23", pi23[i]);
0488             readH5Dataset_double(group_id, "Pi33", pi33[i]);
0489             readH5Dataset_double(group_id, "BulkPi", BulkPi[i]);
0490          }
0491          status = H5Gclose(group_id);
0492       }
0493       else
0494          break;
0495    }
0496 }
0497
0498 void HydroinfoH5::readHydroinfoSingleframe(int frameIdx)
0499 {
0500    hid_t group_id;
0501    herr_t status;
0502
0503    if(frameIdx < (int) grid_Framenum)
0504    {
0505       stringstream frameName;
0506       frameName << "Frame_" <<  setw(4) << setfill('0') << frameIdx;
0507       group_id = H5Gopen(H5groupEventid, frameName.str().c_str(), H5P_DEFAULT);
0508
0509       int Idx = frameIdx;
0510
0511       readH5Dataset_double(group_id, "e", ed[Idx]);
0512       readH5Dataset_double(group_id, "s", sd[Idx]);
0513       readH5Dataset_double(group_id, "Vx", vx[Idx]);
0514       readH5Dataset_double(group_id, "Vy", vy[Idx]);
0515       readH5Dataset_double(group_id, "Temp", Temperature[Idx]);
0516       readH5Dataset_double(group_id, "P", Pressure[Idx]);
0517       if(Visflag == 1)
0518       {
0519          readH5Dataset_double(group_id, "Pi00", pi00[Idx]);
0520          readH5Dataset_double(group_id, "Pi01", pi01[Idx]);
0521          readH5Dataset_double(group_id, "Pi02", pi02[Idx]);
0522          readH5Dataset_double(group_id, "Pi03", pi03[Idx]);
0523          readH5Dataset_double(group_id, "Pi11", pi11[Idx]);
0524          readH5Dataset_double(group_id, "Pi12", pi12[Idx]);
0525          readH5Dataset_double(group_id, "Pi13", pi13[Idx]);
0526          readH5Dataset_double(group_id, "Pi22", pi22[Idx]);
0527          readH5Dataset_double(group_id, "Pi23", pi23[Idx]);
0528          readH5Dataset_double(group_id, "Pi33", pi33[Idx]);
0529          readH5Dataset_double(group_id, "BulkPi", BulkPi[Idx]);
0530       }
0531       status = H5Gclose(group_id);
0532    }
0533    else
0534    {
0535       cout << "Error: readHydroinfoSingleframe :: frameIdx exceed maximum frame number from hydro" << endl;
0536       cout << "frameIdx = " << frameIdx << endl;
0537       exit(1);
0538    }
0539 }
0540
0541 void HydroinfoH5::readH5Dataset_double(hid_t id, string datasetName, double** dset_data)
0542 {
0543    herr_t status;
0544    hid_t dataset_id;
0545
0546    double temp_data[dimensionX][dimensionY];
0547    dataset_id = H5Dopen(id, datasetName.c_str(), H5P_DEFAULT);
0548    status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp_data);
0549    for(int i=0; i<dimensionX; i++)
0550       for(int j=0; j<dimensionY; j++)
0551          dset_data[i][j] = temp_data[i][j];
0552    status = H5Dclose(dataset_id);
0553 }
0554
0555 void HydroinfoH5::getHydroinfoOnlattice(int frameIdx, int xIdx, int yIdx, hydrofluidCell* fluidCellptr)
0556 {
0557    if(frameIdx < 0 || frameIdx > grid_Framenum || xIdx < 0 || xIdx > (grid_XH - grid_XL) || yIdx < 0 || yIdx > (grid_YH - grid_YL))
0558    {
0559       cout << "Error: getHydroinfoOnlattice:: Index is wrong" << endl;
0560       cout << "frameIdx = " << frameIdx << " xIdx = " << xIdx
0561            << "yIdx = " << yIdx << endl;
0562       exit(1);
0563    }
0564    fluidCellptr->ed = ed[frameIdx][xIdx][yIdx];
0565    fluidCellptr->sd = sd[frameIdx][xIdx][yIdx];
0566    fluidCellptr->vx = vx[frameIdx][xIdx][yIdx];
0567    fluidCellptr->vy = vy[frameIdx][xIdx][yIdx];
0568    fluidCellptr->temperature = Temperature[frameIdx][xIdx][yIdx];
0569    fluidCellptr->pressure = Pressure[frameIdx][xIdx][yIdx];
0570    fluidCellptr->pi[0][0] = pi00[frameIdx][xIdx][yIdx];
0571    fluidCellptr->pi[0][1] = pi01[frameIdx][xIdx][yIdx];
0572    fluidCellptr->pi[0][2] = pi02[frameIdx][xIdx][yIdx];
0573    fluidCellptr->pi[0][3] = pi03[frameIdx][xIdx][yIdx];
0574    fluidCellptr->pi[1][0] = fluidCellptr->pi[0][1];
0575    fluidCellptr->pi[1][1] = pi11[frameIdx][xIdx][yIdx];
0576    fluidCellptr->pi[1][2] = pi12[frameIdx][xIdx][yIdx];
0577    fluidCellptr->pi[1][3] = pi13[frameIdx][xIdx][yIdx];
0578    fluidCellptr->pi[2][0] = fluidCellptr->pi[0][2];
0579    fluidCellptr->pi[2][1] = fluidCellptr->pi[1][2];
0580    fluidCellptr->pi[2][2] = pi22[frameIdx][xIdx][yIdx];
0581    fluidCellptr->pi[2][3] = pi23[frameIdx][xIdx][yIdx];
0582    fluidCellptr->pi[3][0] = fluidCellptr->pi[0][3];
0583    fluidCellptr->pi[3][1] = fluidCellptr->pi[1][3];
0584    fluidCellptr->pi[3][2] = fluidCellptr->pi[2][3];
0585    fluidCellptr->pi[3][3] = pi33[frameIdx][xIdx][yIdx];
0586    fluidCellptr->bulkPi = BulkPi[frameIdx][xIdx][yIdx];
0587 }
0588
0589
0590 void HydroinfoH5::getHydroinfo(double tau, double x, double y, hydrofluidCell* fluidCellptr)
0591 {
0592    double eps = 1e-10;
0593    if(tau < grid_Tau0 || tau > grid_Taumax-eps || x < grid_X0 || x > grid_Xmax-eps || y < grid_Y0 || y > grid_Ymax-eps)
0594    {
0595       setZero_fluidCell(fluidCellptr);
0596       return;
0597    }
0598    int frameIdx, xIdx, yIdx;
0599    double tauInc, xInc, yInc;
0600    double temp;
0601
0602    temp = (tau - grid_Tau0)/grid_dTau;
0603    frameIdx = (int) floor(temp);
0604    tauInc = temp - frameIdx;
0605
0606    temp = (x - grid_X0)/grid_dx;
0607    xIdx = (int) floor(temp);
0608    xInc = temp - xIdx;
0609
0610    temp = (y - grid_Y0)/grid_dy;
0611    yIdx = (int) floor(temp);
0612    yInc = temp - yIdx;
0613
0614    fluidCellptr->ed = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, ed);
0615    fluidCellptr->sd = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, sd);
0616    fluidCellptr->vx = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, vx);
0617    fluidCellptr->vy = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, vy);
0618    fluidCellptr->temperature = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, Temperature);
0619    fluidCellptr->pressure = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, Pressure);
0620    if(Visflag == 1)
0621    {
0622       fluidCellptr->pi[0][0] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi00);
0623       fluidCellptr->pi[0][1] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi01);
0624       fluidCellptr->pi[0][2] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi02);
0625       fluidCellptr->pi[0][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi03);
0626       fluidCellptr->pi[1][0] = fluidCellptr->pi[0][1];
0627       fluidCellptr->pi[1][1] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi11);
0628       fluidCellptr->pi[1][2] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi12);
0629       fluidCellptr->pi[1][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi13);
0630       fluidCellptr->pi[2][0] = fluidCellptr->pi[0][2];
0631       fluidCellptr->pi[2][1] = fluidCellptr->pi[1][2];
0632       fluidCellptr->pi[2][2] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi22);
0633       fluidCellptr->pi[2][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi23);
0634       fluidCellptr->pi[3][0] = fluidCellptr->pi[0][3];
0635       fluidCellptr->pi[3][1] = fluidCellptr->pi[1][3];
0636       fluidCellptr->pi[3][2] = fluidCellptr->pi[2][3];
0637       fluidCellptr->pi[3][3] = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, pi33);
0638       fluidCellptr->bulkPi = cubeInterpShell(xIdx, yIdx, frameIdx, xInc, yInc, tauInc, BulkPi);
0639    }
0640    else
0641    {
0642       for(int i=0; i<4; i++)
0643          for(int j=0; j<4; j++)
0644             fluidCellptr->pi[i][j] = 0.0e0;
0645       fluidCellptr->bulkPi = 0.0e0;
0646    }
0647 }
0648
0649 void HydroinfoH5::setZero_fluidCell(hydrofluidCell* fluidCellptr)
0650 {
0651    fluidCellptr->ed = 0.0e0;
0652    fluidCellptr->sd = 0.0e0;
0653    fluidCellptr->vx = 0.0e0;
0654    fluidCellptr->vy = 0.0e0;
0655    fluidCellptr->temperature = 0.0e0;
0656    fluidCellptr->pressure = 0.0e0;
0657    for(int i=0; i<4; i++)
0658       for(int j=0; j<4; j++)
0659          fluidCellptr->pi[i][j] = 0.0e0;
0660    fluidCellptr->bulkPi = 0.0e0;
0661 }
0662
0663 double HydroinfoH5::cubeInterpShell(int idx_x, int idx_y, int idx_z, double x, double y, double z, double ***dataset)
0664 {
0665    double result = cubeInterp(x, y, z,
0666       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],
0667       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]);
0668    return(result);
0669 }
0670
0671 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)
0672 // Perform a 3d interpolation. The known data are A### located at the 8 corners,
0673 // labels using the xyz order. Therefore A000 is value at the origin and A010
0674 // is the value at (x=0,y=1,z=0). Note that the coordinate (x,y,z) must be
0675 // constrained to the unit cube. Axyz is the return value.
0676 {
0677    /* for debug
0678    cout << A000 << "  " << A100 << "  " << A010 << "  " << A110 << endl;
0679    cout << A001 << "  " << A101 << "  " << A011 << "  " << A111 << endl; */
0680
0681    double Axyz = A000*(1-x)*(1-y)*(1-z) + A100*x*(1-y)*(1-z)
0682                  + A010*(1-x)*y*(1-z) + A001*(1-x)*(1-y)*z
0683                  + A101*x*(1-y)*z + A011*(1-x)*y*z + A110*x*y*(1-z)
0684                  + A111*x*y*z;
0685    return(Axyz);
0686 }