Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:19:44 

0001 // Hydroinfo_MUSIC.cpp is a part of the MARTINI event generator.
0002 // Copyright (C) 2009-2010 Bjoern Schenke.
0003 // MARTINI is licenced under the GNU GPL version 2, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 
0006 // This file contains routines to read in hydro data from files and functions
0007 // that return interpolated data at a given space-time point
0008 
0009 #include <iostream>
0010 #include <fstream>
0011 #include <iomanip>
0012 #include <cstdlib>
0013 #include <cmath>
0014 #include <sstream>
0015 #include <vector>
0016 #include <string>
0017 
0018 #include "./Hydroinfo_MUSIC.h"
0019 
0020 using namespace std;
0021 
0022 Hydroinfo_MUSIC::Hydroinfo_MUSIC() {
0023     verbose_ = 9;
0024     hbarC = 0.19733;
0025     boost_invariant = false;
0026 }
0027 
0028 Hydroinfo_MUSIC::~Hydroinfo_MUSIC() {
0029     clean_hydro_event();
0030 }
0031 
0032 void Hydroinfo_MUSIC::clean_hydro_event() {
0033     if (boost_invariant) {
0034         lattice_2D.clear();
0035     } else {
0036         lattice_3D.clear();
0037         lattice_3D_ideal.clear();
0038     }
0039 }
0040 
0041 void Hydroinfo_MUSIC::readHydroData(int whichHydro, int nskip_tau_in,
0042                                     string input_filename_in,
0043                                     string hydro_ideal_filename_in,
0044                                     string hydro_shear_filename_in,
0045                                     string hydro_bulk_filename_in) {
0046     // all hydro data is stored in tau steps (not t)
0047     // evolution is converted to tau when accessing the hydro data
0048     lattice_2D.clear();
0049     lattice_3D.clear();
0050     lattice_3D_ideal.clear();
0051 
0052     input_filename = input_filename_in;
0053     hydro_ideal_filename = hydro_ideal_filename_in;
0054     hydro_shear_filename = hydro_shear_filename_in;
0055     hydro_bulk_filename = hydro_bulk_filename_in;
0056 
0057     // read in setups of the hydro simulation
0058     hydroWhichHydro = whichHydro;
0059     if (hydroWhichHydro < 10) {
0060         ostringstream config_file;
0061         config_file << input_filename;
0062         ifstream configuration;
0063         configuration.open(config_file.str().c_str(), ios::in);
0064         if (!configuration) {
0065             cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0066                  << "Unable to open file: " << config_file.str() << endl;
0067             exit(1);
0068         }
0069         string temp1;
0070         string temp_name;
0071         while (!configuration.eof()) {
0072             getline(configuration, temp1);
0073             stringstream ss(temp1);
0074             ss >> temp_name;
0075 
0076             // read in grid information
0077             if (temp_name == "Initial_time_tau_0") {
0078                 ss >> hydroTau0;
0079             } else if (temp_name == "Delta_Tau") {
0080                 ss >> hydroDtau;
0081             } else if (temp_name == "X_grid_size_in_fm") {
0082                 double temp;
0083                 ss >> temp;
0084                 hydroXmax = temp/2.;
0085             } else if (temp_name == "Grid_size_in_x") {
0086                 ss >> ixmax;
0087             } else if (temp_name == "Eta_grid_size") {
0088                 double temp;
0089                 ss >> temp;
0090                 hydro_eta_max = temp/2.;
0091             } else if (temp_name == "Grid_size_in_eta") {
0092                 ss >> ietamax;
0093             } else if (temp_name == "output_evolution_every_N_timesteps") {
0094                 ss >> nskip_tau;
0095             } else if (temp_name == "output_evolution_every_N_x") {
0096                 ss >> nskip_x;
0097             } else if (temp_name == "output_evolution_every_N_eta") {
0098                 ss >> nskip_eta;
0099             }
0100             // read in additioinal information
0101             if (temp_name == "Include_Rhob_Yes_1_No_0") {
0102                 ss >> turn_on_rhob;
0103             } else if (temp_name == "Include_Shear_Visc_Yes_1_No_0") {
0104                 ss >> turn_on_shear;
0105             } else if (temp_name == "Include_Bulk_Visc_Yes_1_No_0") {
0106                 ss >> turn_on_bulk;
0107             } else if (temp_name == "turn_on_baryon_diffusion") {
0108                 ss >> turn_on_diff;
0109             }
0110         }
0111         configuration.close();
0112         hydroDx = 2.*hydroXmax/(ixmax - 1.);
0113         hydroDeta = 2.*hydro_eta_max/(static_cast<double>(ietamax));
0114     }
0115 
0116     use_tau_eta_coordinate = 1;
0117     if (use_tau_eta_coordinate == 0) {
0118         cout << "Hydroinfo_MUSIC:: Warning hydro grid is set to "
0119              << "cartesian coordinates, please make sure this is correct!"
0120              << endl;
0121     }
0122 
0123     if (whichHydro == 6) {
0124         // 3+1D MUSIC hydro (Schenke, Jeon, Gale)
0125         cout << "Using 3+1D Jeon Schenke hydro reading data ..." << endl;
0126         boost_invariant = false;
0127 
0128         ixmax = static_cast<int>(2.*hydroXmax/hydroDx + 0.001);
0129         ietamax = static_cast<int>(2.*hydro_eta_max/hydroDeta + 0.001);
0130 
0131         // read in temperature, QGP fraction , flow velocity
0132         // The name of the evolution file: evolution_name
0133         string evolution_name = hydro_ideal_filename;
0134         cout << "Evolution file name = " << evolution_name << endl;
0135         ifstream fin;
0136         fin.open(evolution_name.c_str(), ios::in);
0137         if (!fin) {
0138             cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0139                  << "Unable to open file: " << evolution_name << endl;
0140             exit(1);
0141         }
0142 
0143         double T, vx, vy, vz, QGPfrac;
0144         fluidCell_3D newCell;
0145         int ik = 0;
0146         while (!fin.eof()) {
0147             ik++;
0148             fin >> T;
0149             fin >> QGPfrac;
0150             fin >> vx;
0151             fin >> vy;
0152             fin >> vz;
0153             newCell.temperature = T;
0154             newCell.vx = vx;
0155             newCell.vy = vy;
0156             newCell.vz = vz;
0157 
0158             lattice_3D.push_back(newCell);
0159             if (ik%50000 == 0)
0160                 cout << "o" << flush;
0161         }
0162         cout << ik << endl;
0163         fin.close();
0164     } else if (whichHydro == 8) {
0165         // event-by-event (2+1)-d MUSIC hydro from JF
0166         // there are two slices in medium in eta_s
0167         // one at eta_s = -15. and the other at eta_s = 0.0
0168         // only the medium at middle rapidity will be kept in the memory
0169         boost_invariant = true;
0170         cout << "Reading event-by-event hydro evolution data "
0171              << "from standard MUSIC ..." << endl;
0172 
0173         ixmax = static_cast<int>(2.*hydroXmax/hydroDx + 0.001);
0174         ietamax = 1;
0175         nskip_tau = nskip_tau_in;
0176 
0177         hydroDx *= nskip_x;
0178         hydroDtau *= nskip_tau;
0179         hydroDeta *= nskip_eta;
0180 
0181         int n_eta = 2;  // there are two slices in eta_s
0182         // number of fluid cell in the transverse plane
0183         int num_fluid_cell_trans = ixmax*ixmax;
0184 
0185         // read in hydro evolution
0186         string evolution_name = hydro_ideal_filename;
0187         string evolution_name_Wmunu = hydro_shear_filename;
0188         string evolution_name_Pi = hydro_bulk_filename;
0189 
0190         std::FILE *fin;
0191         string evolution_file_name = evolution_name;
0192         cout << "Evolution file name = " << evolution_file_name << endl;
0193         fin = std::fopen(evolution_file_name.c_str(), "rb");
0194 
0195         if (fin == NULL) {
0196             cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0197                  << "Unable to open file: " << evolution_file_name << endl;
0198             exit(1);
0199         }
0200 
0201         std::FILE *fin1 = NULL;
0202         if (turn_on_shear == 1) {
0203             fin1 = std::fopen(evolution_name_Wmunu.c_str(), "rb");
0204             if (fin1 == NULL) {
0205                 cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0206                      << "Unable to open file: " << evolution_name_Wmunu << endl;
0207                 exit(1);
0208             }
0209         }
0210 
0211         std::FILE *fin2 = NULL;
0212         if (turn_on_bulk == 1) {
0213             fin2 = std::fopen(evolution_name_Pi.c_str(), "rb");
0214             if (fin2 == NULL) {
0215                 cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0216                      << "Unable to open file: " << evolution_name_Pi << endl;
0217                 exit(1);
0218             }
0219         }
0220 
0221         int ik = 0;
0222         fluidCell_2D newCell;
0223         float T, QGPfrac, vx, vy, vz;
0224         double ux, uy, ueta;
0225         float pi00 = 0.0;
0226         float pi01 = 0.0;
0227         float pi02 = 0.0;
0228         float pi03 = 0.0;
0229         float pi11 = 0.0;
0230         float pi12 = 0.0;
0231         float pi13 = 0.0;
0232         float pi22 = 0.0;
0233         float pi23 = 0.0;
0234         float pi33 = 0.0;;
0235         float bulkPi = 0.0;
0236         float e_plus_P = 1e-15;
0237         float cs2 = 0.0;
0238         int size = sizeof(float);
0239         while (true) {
0240             int status = 0;
0241             status = std::fread(&T, size, 1, fin);
0242             status += std::fread(&QGPfrac, size, 1, fin);
0243             status += std::fread(&vx, size, 1, fin);
0244             status += std::fread(&vy, size, 1, fin);
0245             status += std::fread(&vz, size, 1, fin);
0246 
0247             if (status != 5) {  // this is the end of file
0248                 break;
0249             }
0250 
0251             int status_pi = 0;
0252             if (turn_on_shear == 1) {
0253                 status_pi = std::fread(&pi00, size, 1, fin1);
0254                 status_pi += std::fread(&pi01, size, 1, fin1);
0255                 status_pi += std::fread(&pi02, size, 1, fin1);
0256                 status_pi += std::fread(&pi03, size, 1, fin1);
0257                 status_pi += std::fread(&pi11, size, 1, fin1);
0258                 status_pi += std::fread(&pi12, size, 1, fin1);
0259                 status_pi += std::fread(&pi13, size, 1, fin1);
0260                 status_pi += std::fread(&pi22, size, 1, fin1);
0261                 status_pi += std::fread(&pi23, size, 1, fin1);
0262                 status_pi += std::fread(&pi33, size, 1, fin1);
0263 
0264                 if (status_pi != 10) {
0265                     cout << "Error:Hydroinfo_MUSIC::readHydroData: "
0266                          << "Wmunu file does not have the same number of "
0267                          << "fluid cells as the ideal file!" << endl;
0268                     exit(1);
0269                 }
0270             }
0271 
0272             int status_bulkPi = 0;
0273             if (turn_on_bulk == 1) {
0274                 status_bulkPi = std::fread(&bulkPi, size, 1, fin2);
0275                 status_bulkPi += std::fread(&e_plus_P, size, 1, fin2);
0276                 status_bulkPi += std::fread(&cs2, size, 1, fin2);
0277 
0278                 if (status_bulkPi != 3) {
0279                     cout << "Error:Hydroinfo_MUSIC::readHydroData: "
0280                          << "bulkPi file does not have the same number of "
0281                          << "fluid cells as the ideal file!" << endl;
0282                     exit(1);
0283                 }
0284             }
0285 
0286             int ieta_idx = static_cast<int>(ik/num_fluid_cell_trans) % n_eta;
0287             int itau_idx = static_cast<int>(ik/(num_fluid_cell_trans*n_eta));
0288             ik++;
0289             if (itau_idx%nskip_tau != 0)  // skip in tau
0290                 continue;
0291 
0292             // print out tau information
0293             double tau_local = hydroTau0 + itau_idx*hydroDtau/nskip_tau;
0294             if ((ik-1)%(num_fluid_cell_trans*n_eta) == 0) {
0295                 cout << "read in tau frame: " << itau_idx
0296                      << " tau_local = " << setprecision(3) << tau_local
0297                      << " fm ..."<< endl;
0298             }
0299 
0300             if (ieta_idx == (n_eta-1)) {
0301                 // store the hydro medium at eta_s = 0.0
0302                 double v2 = vx*vx + vy*vy + vz*vz;
0303                 if (v2 > 1.0) {
0304                     cerr << "[Hydroinfo_MUSIC::readHydroData:] Error: "
0305                          << "v > 1! vx = " << vx << ", vy = " << vy
0306                          << ", vz = " << vz << ", T = " << T << endl;
0307                     if (T > 0.01) {
0308                         exit(1);
0309                     } else {
0310                         v2 = 0.0;
0311                     }
0312                 }
0313                 double gamma = 1./sqrt(1. - v2);
0314                 ux = gamma*vx;
0315                 uy = gamma*vy;
0316                 ueta = gamma*vz;  // assuming eta = 0
0317 
0318                 newCell.temperature = T;
0319                 // convert vx and vy to longitudinal co-moving frame
0320                 newCell.ux = ux;
0321                 newCell.uy = uy;
0322                 newCell.ueta = ueta;
0323 
0324                 // pi^\mu\nu tensor
0325                 newCell.pi00 = pi00;
0326                 newCell.pi01 = pi01;
0327                 newCell.pi02 = pi02;
0328                 newCell.pi11 = pi11;
0329                 newCell.pi12 = pi12;
0330                 newCell.pi22 = pi22;
0331                 newCell.pi33 = pi33;
0332 
0333                 // bulk pressure
0334                 if (T > 0.18) {
0335                     // QGP phase prefactor is divided out here
0336                     newCell.bulkPi = bulkPi/(15.*(1./3. - cs2)*e_plus_P);
0337                 } else {
0338                     newCell.bulkPi = bulkPi;   // [1/fm^4]
0339                 }
0340                 lattice_2D.push_back(newCell);
0341             }
0342         }
0343         std::fclose(fin);
0344         if (turn_on_shear == 1) {
0345             std::fclose(fin1);
0346         }
0347         if (turn_on_bulk == 1) {
0348             std::fclose(fin2);
0349         }
0350         cout << endl;
0351         cout << "number of fluid cells: " << lattice_2D.size() << endl;
0352     } else if (whichHydro == 9) {
0353         // event-by-event (2+1)-d MUSIC hydro
0354         // the output medium is at middle rapidity
0355         boost_invariant = true;
0356         cout << "Reading event-by-event hydro evolution data "
0357              << "from (2+1)D MUSIC ..." << endl;
0358 
0359         ietamax = 1;
0360 
0361         hydroDx *= nskip_x;
0362         hydroDtau *= nskip_tau;
0363         hydroDeta *= nskip_eta;
0364 
0365         nskip_tau = nskip_tau_in;
0366         ixmax = static_cast<int>(2.*hydroXmax/hydroDx + 0.001);
0367 
0368         int n_eta = 1;
0369         // number of fluid cell in the transverse plane
0370         int num_fluid_cell_trans = ixmax*ixmax;
0371 
0372         // read in hydro evolution
0373         string evolution_name = hydro_ideal_filename;
0374         string evolution_name_Wmunu = hydro_shear_filename;
0375         string evolution_name_Pi = hydro_bulk_filename;
0376 
0377         std::FILE *fin;
0378         string evolution_file_name = evolution_name;
0379         cout << "Evolution file name = " << evolution_file_name << endl;
0380         fin = std::fopen(evolution_file_name.c_str(), "rb");
0381 
0382         if (fin == NULL) {
0383             cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0384                  << "Unable to open file: " << evolution_file_name << endl;
0385             exit(1);
0386         }
0387 
0388         std::FILE *fin1 = NULL;
0389         if (turn_on_shear == 1) {
0390             fin1 = std::fopen(evolution_name_Wmunu.c_str(), "rb");
0391             if (fin1 == NULL) {
0392                 cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0393                      << "Unable to open file: " << evolution_name_Wmunu << endl;
0394                 exit(1);
0395             }
0396         }
0397 
0398         std::FILE *fin2 = NULL;
0399         if (turn_on_bulk == 1) {
0400             fin2 = std::fopen(evolution_name_Pi.c_str(), "rb");
0401             if (fin2 == NULL) {
0402                 cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0403                      << "Unable to open file: " << evolution_name_Pi << endl;
0404                 exit(1);
0405             }
0406         }
0407 
0408         int ik = 0;
0409         fluidCell_2D newCell;
0410         double T, QGPfrac, ux, uy, ueta;
0411         double vx, vy, vz;
0412         float pi00 = 0.0;
0413         float pi01 = 0.0;
0414         float pi02 = 0.0;
0415         float pi03 = 0.0;
0416         float pi11 = 0.0;
0417         float pi12 = 0.0;
0418         float pi13 = 0.0;
0419         float pi22 = 0.0;
0420         float pi23 = 0.0;
0421         float pi33 = 0.0;;
0422         float bulkPi = 0.0;
0423         float e_plus_P = 1e-15;
0424         float cs2 = 0.0;
0425         int size = sizeof(double);
0426         while (true) {
0427             int status = 0;
0428             status = std::fread(&T, size, 1, fin);
0429             status += std::fread(&QGPfrac, size, 1, fin);
0430             status += std::fread(&vx, size, 1, fin);
0431             status += std::fread(&vy, size, 1, fin);
0432             status += std::fread(&vz, size, 1, fin);
0433             if (status != 5) {  // this is the end of file
0434                 break;
0435             }
0436 
0437             double v2 = vx*vx + vy*vy + vz*vz;
0438             if (v2 > 1.) {
0439                 cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0440                      << "v > 1! vx = " << vx << ", vy = " << vy
0441                      << ", vz = " << vz << endl;
0442                 exit(1);
0443             }
0444             double gamma = 1./sqrt(1. - v2);
0445             ux = vx*gamma;
0446             uy = vy*gamma;
0447             ueta = vz*gamma;  // assuming at the eta = 0
0448 
0449             int status_pi = 0;
0450             if (turn_on_shear == 1) {
0451                 status_pi = std::fread(&pi00, size, 1, fin1);
0452                 status_pi += std::fread(&pi01, size, 1, fin1);
0453                 status_pi += std::fread(&pi02, size, 1, fin1);
0454                 status_pi += std::fread(&pi03, size, 1, fin1);
0455                 status_pi += std::fread(&pi11, size, 1, fin1);
0456                 status_pi += std::fread(&pi12, size, 1, fin1);
0457                 status_pi += std::fread(&pi13, size, 1, fin1);
0458                 status_pi += std::fread(&pi22, size, 1, fin1);
0459                 status_pi += std::fread(&pi23, size, 1, fin1);
0460                 status_pi += std::fread(&pi33, size, 1, fin1);
0461 
0462                 if (status_pi != 10) {
0463                     cout << "Error:Hydroinfo_MUSIC::readHydroData: "
0464                          << "Wmunu file does not have the same number of "
0465                          << "fluid cells as the ideal file!" << endl;
0466                     exit(1);
0467                 }
0468             }
0469 
0470             int status_bulkPi = 0;
0471             if (turn_on_bulk == 1) {
0472                 status_bulkPi = std::fread(&bulkPi, size, 1, fin2);
0473                 status_bulkPi += std::fread(&e_plus_P, size, 1, fin2);
0474                 status_bulkPi += std::fread(&cs2, size, 1, fin2);
0475 
0476                 if (status_bulkPi != 3) {
0477                     cout << "Error:Hydroinfo_MUSIC::readHydroData: "
0478                          << "bulkPi file does not have the same number of "
0479                          << "fluid cells as the ideal file!" << endl;
0480                     exit(1);
0481                 }
0482             }
0483 
0484             int ieta_idx = static_cast<int>(ik/num_fluid_cell_trans) % n_eta;
0485             int itau_idx = static_cast<int>(ik/(num_fluid_cell_trans*n_eta));
0486             ik++;
0487             if (itau_idx%nskip_tau != 0)  // skip in tau
0488                 continue;
0489 
0490             // print out tau information
0491             double tau_local = hydroTau0 + itau_idx*hydroDtau/nskip_tau;
0492             if ((ik-1)%(num_fluid_cell_trans*n_eta) == 0) {
0493                 cout << "read in tau frame: " << itau_idx
0494                      << " tau_local = " << setprecision(3) << tau_local
0495                      << " fm ..."<< endl;
0496             }
0497 
0498             if (ieta_idx == (n_eta-1)) {
0499                 newCell.temperature = T;
0500                 newCell.ux = ux;
0501                 newCell.uy = uy;
0502                 newCell.ueta = ueta;
0503 
0504                 // pi^\mu\nu tensor
0505                 newCell.pi00 = pi00;
0506                 newCell.pi01 = pi01;
0507                 newCell.pi02 = pi02;
0508                 newCell.pi11 = pi11;
0509                 newCell.pi12 = pi12;
0510                 newCell.pi22 = pi22;
0511                 newCell.pi33 = pi33;
0512 
0513                 // bulk pressure
0514                 if (T > 0.18) {
0515                     // QGP phase prefactor is divided out here
0516                     newCell.bulkPi = bulkPi/(15.*(1./3. - cs2)*e_plus_P);
0517                 } else {
0518                     newCell.bulkPi = bulkPi;   // [1/fm^4]
0519                 }
0520                 lattice_2D.push_back(newCell);
0521             }
0522         }
0523         std::fclose(fin);
0524         if (turn_on_shear == 1) {
0525             std::fclose(fin1);
0526         }
0527         if (turn_on_bulk == 1) {
0528             std::fclose(fin2);
0529         }
0530         cout << endl;
0531         cout << "number of fluid cells: " << lattice_2D.size() << endl;
0532     } else if (whichHydro == 10 || whichHydro == 11) {
0533         // new MUSIC hydro (no regular grid)
0534         cout << "Using 3+1D new MUSIC hydro reading data ..." << endl;
0535         if (whichHydro == 10) {
0536             boost_invariant = false;
0537         } else {
0538             boost_invariant = true;
0539         }
0540 
0541         // read in temperature and flow velocity
0542         // The name of the evolution file: evolution_name
0543         string evolution_name = hydro_ideal_filename;
0544         cout << "Evolution file name = " << evolution_name << endl;
0545         std::FILE *fin;
0546         fin = std::fopen(evolution_name.c_str(), "rb");
0547         if (fin == NULL) {
0548             cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0549                  << "Unable to open file: " << evolution_name << endl;
0550             exit(1);
0551         }
0552 
0553         float header[16];
0554         int status = std::fread(&header, sizeof(float), 16, fin);
0555         if (status == 0) {
0556             cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0557                  << "Can not read the evolution file header" << endl;
0558             exit(1);
0559         }
0560 
0561         hydroTau0 = header[0];
0562         hydroDtau = header[1];
0563         ixmax = static_cast<int>(header[2]);
0564         hydroDx = header[3];
0565         hydroXmax = std::abs(header[4]);
0566         ietamax = static_cast<int>(header[8]);
0567         hydroDeta = header[9];
0568         hydro_eta_max = std::abs(header[10]);
0569         turn_on_rhob = static_cast<int>(header[11]);
0570         turn_on_shear = static_cast<int>(header[12]);
0571         turn_on_bulk = static_cast<int>(header[13]);
0572         turn_on_diff = static_cast<int>(header[14]);
0573         const int nVar_per_cell = static_cast<int>(header[15]);
0574 
0575         float cell_info[nVar_per_cell];
0576 
0577         int itau_max = 0;
0578         fluidCell_3D_ideal zeroCell;
0579         zeroCell.itau = 0;
0580         zeroCell.ix = 0;
0581         zeroCell.iy = 0;
0582         zeroCell.ieta = 0;
0583         zeroCell.temperature = 0.;
0584         zeroCell.ed = 0.;
0585         zeroCell.pressure = 0.;
0586         zeroCell.ux = 0.;
0587         zeroCell.uy = 0.;
0588         zeroCell.uz = 0.;
0589         lattice_3D_ideal.push_back(zeroCell);
0590         int ik = 0;
0591         while (true) {
0592             status = 0;
0593             status = std::fread(&cell_info, sizeof(float), nVar_per_cell, fin);
0594             if (status == 0) break;
0595             if (status != nVar_per_cell) {
0596                 cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
0597                      << "the evolution file format is not correct" << endl;
0598                 exit(1);
0599             }
0600 
0601             if (itau_max < static_cast<int>(cell_info[0]))
0602                 itau_max = static_cast<int>(cell_info[0]);
0603             fluidCell_3D_ideal newCell;
0604             newCell.itau = static_cast<int>(cell_info[0]);
0605             newCell.ix   = static_cast<int>(cell_info[1]);
0606             newCell.iy   = static_cast<int>(cell_info[2]);
0607             newCell.ieta = static_cast<int>(cell_info[3]);
0608             newCell.temperature = cell_info[6];
0609             newCell.ed = cell_info[4];
0610             newCell.pressure = cell_info[5];
0611             newCell.ux = cell_info[8];
0612             newCell.uy = cell_info[9];
0613             newCell.uz = cell_info[10];
0614             lattice_3D_ideal.push_back(newCell);
0615             ik++;
0616             if (verbose_ > 7) {
0617                 if (ik%50000 == 0)
0618                     cout << "o" << flush;
0619             }
0620         }
0621         cout << endl;
0622         std::fclose(fin);
0623         itaumax = itau_max;
0624         // create the index map
0625         long long ncells = (itaumax + 1)*ixmax*ixmax*ietamax;
0626         idx_map_.resize(ncells, 0);
0627         for (int i = 0; i < lattice_3D_ideal.size(); i++) {
0628             const auto cell_i = lattice_3D_ideal[i];
0629             long long cell_idx = (
0630                 (  (cell_i.itau*ietamax + cell_i.ieta)*ixmax
0631                  + cell_i.iy)*ixmax + cell_i.ix);
0632             idx_map_[cell_idx] = i;
0633         }
0634         hydroTauMax = hydroTau0 + hydroDtau*itaumax;
0635     } else {
0636         cout << "Hydroinfo_MUSIC:: This option is obsolete! whichHydro = "
0637              << whichHydro << endl;
0638         exit(1);
0639     }
0640 
0641     // One final step for easy automation of MARTINI:
0642     // hydroTauMax is reset for the case where writing to evolution.dat
0643     // ended early (due to all cells freezing out):
0644     if (whichHydro == 6) {
0645         hydroTauMax = (
0646             hydroTau0 + hydroDtau*static_cast<int>(
0647                         static_cast<double>(lattice_3D.size())
0648                         /((2.*hydroXmax/hydroDx+1.)*(2.*hydroXmax/hydroDx+1.)
0649                         *2.*(hydro_eta_max/hydroDeta))));
0650         itaumax = static_cast<int>((hydroTauMax-hydroTau0)/hydroDtau+0.001);
0651     }
0652     if (whichHydro == 8 || whichHydro == 9) {
0653         hydroTauMax = (
0654             hydroTau0 + hydroDtau*static_cast<int>(
0655                         static_cast<double>(lattice_2D.size())
0656                         /((2.*hydroXmax/hydroDx)*(2.*hydroXmax/hydroDx)) - 1));
0657         itaumax = static_cast<int>((hydroTauMax - hydroTau0)/hydroDtau);
0658     }
0659 
0660     cout << "hydro_tau0 = " << hydroTau0 << " fm"<< endl;
0661     cout << "hydro_tau_max = " << hydroTauMax << " fm" << endl;
0662     cout << "hydry_dtau = " << hydroDtau << " fm" << endl;
0663     cout << "hydro_Xmax = " << hydroXmax << " fm" << endl;
0664     cout << "hydro_dx = " << hydroDx << " fm" << endl;
0665     cout << "hydro_eta_max = " << hydro_eta_max << endl;
0666     cout << "hydro_deta = " << hydroDeta << endl;
0667 }
0668 
0669 
0670 void Hydroinfo_MUSIC::getHydroValues(double x, double y,
0671                                      double z, double t, hydrofluidCell* info) {
0672 // For interpolation of evolution files in tau-eta coordinates. Only the
0673 // reading of MUSIC's evolution_xyeta.dat file is implemented here.
0674 // For simplicity, hydro_eta_max refers to MUSIC's eta_size, and similarly for
0675 // hydroDeta; however, x, y, z, and t are as usual to stay compatible with
0676 // MARTINI.
0677     double tau, eta;
0678     if (use_tau_eta_coordinate == 1) {
0679         if (t*t > z*z) {
0680             tau = sqrt(t*t-z*z);
0681             eta = 0.5*log((t+z)/(t-z));
0682         } else {
0683             tau = 0.;
0684             eta = 0.;
0685         }
0686     } else {
0687         // if the medium is given in cartesian coordinates
0688         // set tau and eta to t and z
0689         tau = t;
0690         eta = z;
0691     }
0692 
0693     int ieta = floor((hydro_eta_max+eta)/hydroDeta + 0.0001);
0694     int itau = floor((tau-hydroTau0)/hydroDtau + 0.0001);
0695     int ix = floor((hydroXmax+x)/hydroDx + 0.0001);
0696     int iy = floor((hydroXmax+y)/hydroDx + 0.0001);
0697 
0698     double xfrac = (x - (static_cast<double>(ix)*hydroDx - hydroXmax))/hydroDx;
0699     double yfrac = (y - (static_cast<double>(iy)*hydroDx - hydroXmax))/hydroDx;
0700     double etafrac = (eta/hydroDeta - static_cast<double>(ieta)
0701                       + 0.5*static_cast<double>(ietamax));
0702     double taufrac = (tau - hydroTau0)/hydroDtau - static_cast<double>(itau);
0703 
0704     if (boost_invariant) {
0705         ieta = 0;
0706         etafrac = 0.;
0707     }
0708 
0709     if (ix < 0 || ix >= ixmax) {
0710         if (verbose_ > 7) {
0711             cout << "[Hydroinfo_MUSIC::getHydroValues]: "
0712                  << "WARNING - x out of range x=" << x
0713                  << ", ix=" << ix << ", ixmax=" << ixmax << endl;
0714             cout << "x=" << x << " y=" << y << " eta=" << eta
0715                  << " ix=" << ix << " iy=" << iy << " ieta=" << ieta << endl;
0716             cout << "t=" << t << " tau=" << tau
0717                  << " itau=" << itau << " itaumax=" << itaumax << endl;
0718         }
0719 
0720         info->temperature = 0.0;
0721         info->vx = 0.0;
0722         info->vy = 0.0;
0723         info->vz = 0.0;
0724         return;
0725     }
0726     if (iy < 0 || iy >= ixmax) {
0727         if (verbose_ > 7) {
0728             cout << "[Hydroinfo_MUSIC::getHydroValues]: "
0729                  << "WARNING - y out of range, y=" << y << ", iy="  << iy
0730                  << ", iymax=" << ixmax << endl;
0731             cout << "x=" << x << " y=" << y << " eta=" << eta
0732                  << " ix=" << ix << " iy=" << iy << " ieta=" << ieta << endl;
0733             cout << "t=" << t << " tau=" << tau
0734                  << " itau=" << itau << " itaumax=" << itaumax << endl;
0735         }
0736 
0737         info->temperature = 0.0;
0738         info->vx = 0.0;
0739         info->vy = 0.0;
0740         info->vz = 0.0;
0741         return;
0742     }
0743     if (itau < 0 || itau > itaumax) {
0744         if (verbose_ > 7) {
0745             cout << "[Hydroinfo_MUSIC::getHydroValues]: WARNING - "
0746                  << "tau out of range, itau=" << itau
0747                  << ", itaumax=" << itaumax << endl;
0748             cout << "[MARTINI:Hydroinfo_MUSIC::getHydroValues]: tau= " << tau
0749                  << ", hydroTauMax = " << hydroTauMax
0750                  << ", hydroDtau = " << hydroDtau << endl;
0751         }
0752 
0753         info->temperature = 0.0;
0754         info->vx = 0.0;
0755         info->vy = 0.0;
0756         info->vz = 0.0;
0757         return;
0758     }
0759     if (ieta < 0 || ieta >= ietamax) {
0760         if (verbose_ > 7) {
0761             cout << "[Hydroinfo_MUSIC::getHydroValues]: WARNING - "
0762                  << "eta out of range, ieta=" << ieta
0763                  << ", ietamax=" << ietamax
0764                  << endl;
0765         }
0766         info->temperature = 0.0;
0767         info->vx = 0.0;
0768         info->vy = 0.0;
0769         info->vz = 0.0;
0770         return;
0771     }
0772 
0773     // The array of positions on the 4-dimensional rectangle:
0774     int position[2][2][2][2];
0775     for (int ipx = 0; ipx < 2; ipx++) {
0776         int px;
0777         if (ipx == 0 || ix == ixmax-1) {
0778             px = ix;
0779         } else {
0780             px = ix + 1;
0781         }
0782         for (int ipy = 0; ipy < 2; ipy++) {
0783             int py;
0784             if (ipy == 0 || iy == ixmax-1) {
0785                 py = iy;
0786             } else {
0787                 py = iy + 1;
0788             }
0789             for (int ipeta = 0; ipeta < 2; ipeta++) {
0790                 int peta;
0791                 if (ipeta == 0 || ieta == ietamax-1) {
0792                     peta = ieta;
0793                 } else {
0794                     peta = ieta + 1;
0795                 }
0796                 for (int iptau = 0; iptau < 2; iptau++) {
0797                     int ptau;
0798                     if (iptau == 0 || itau == itaumax) {
0799                         ptau = itau;
0800                     } else {
0801                         ptau = itau + 1;
0802                     }
0803                     position[ipx][ipy][ipeta][iptau] = (
0804                                 px + ixmax*(py + ixmax*(peta + ietamax*ptau)));
0805                 }
0806             }
0807         }
0808     }
0809 
0810     // And now, the interpolation:
0811     double T = 0.0;
0812     double ed = 0.;
0813     double p = 0.;
0814     double vx = 0.0;
0815     double vy = 0.0;
0816     double vz = 0.0;
0817     double ux = 0.0;
0818     double uy = 0.0;
0819     double uz = 0.0;
0820     double ueta = 0.0;
0821     double pi00 = 0.0;
0822     double pi01 = 0.0;
0823     double pi02 = 0.0;
0824     double pi03 = 0.0;
0825     double pi11 = 0.0;
0826     double pi12 = 0.0;
0827     double pi13 = 0.0;
0828     double pi22 = 0.0;
0829     double pi23 = 0.0;
0830     double pi33 = 0.0;
0831     double bulkPi = 0.0;
0832 
0833     fluidCell_2D *HydroCell_2D_ptr1, *HydroCell_2D_ptr2;
0834     fluidCell_3D *HydroCell_3D_ptr1, *HydroCell_3D_ptr2;
0835     fluidCell_3D_ideal *HydroCell_3D_ideal_ptr1, *HydroCell_3D_ideal_ptr2;
0836     for (int iptau = 0; iptau < 2; iptau++) {
0837         double taufactor;
0838         if (iptau == 0)
0839             taufactor = 1. - taufrac;
0840         else
0841             taufactor = taufrac;
0842         for (int ipeta = 0; ipeta < 2; ipeta++) {
0843             double etafactor;
0844             if (ipeta == 0)
0845                 etafactor = 1. - etafrac;
0846             else
0847                 etafactor = etafrac;
0848             for (int ipy = 0; ipy < 2; ipy++) {
0849                 double yfactor;
0850                 if (ipy == 0)
0851                     yfactor = 1. - yfrac;
0852                 else
0853                     yfactor = yfrac;
0854 
0855                 double prefrac = yfactor*etafactor*taufactor;
0856 
0857                 if (hydroWhichHydro == 8 || hydroWhichHydro == 9) {
0858                     HydroCell_2D_ptr1 = (
0859                             &lattice_2D[position[0][ipy][ipeta][iptau]]);
0860                     HydroCell_2D_ptr2 = (
0861                             &lattice_2D[position[1][ipy][ipeta][iptau]]);
0862                     T += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->temperature
0863                                   + xfrac*HydroCell_2D_ptr2->temperature);
0864                     ux += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->ux
0865                                     + xfrac*HydroCell_2D_ptr2->ux);
0866                     uy += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->uy
0867                                     + xfrac*HydroCell_2D_ptr2->uy);
0868                     ueta += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->ueta
0869                                     + xfrac*HydroCell_2D_ptr2->ueta);
0870                     pi00 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi00
0871                                     + xfrac*HydroCell_2D_ptr2->pi00);
0872                     pi01 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi01
0873                                     + xfrac*HydroCell_2D_ptr2->pi01);
0874                     pi02 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi02
0875                                     + xfrac*HydroCell_2D_ptr2->pi02);
0876                     pi11 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi11
0877                                     + xfrac*HydroCell_2D_ptr2->pi11);
0878                     pi12 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi12
0879                                     + xfrac*HydroCell_2D_ptr2->pi12);
0880                     pi22 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi22
0881                                     + xfrac*HydroCell_2D_ptr2->pi22);
0882                     pi33 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi33
0883                                     + xfrac*HydroCell_2D_ptr2->pi33);
0884                     bulkPi += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->bulkPi
0885                                     + xfrac*HydroCell_2D_ptr2->bulkPi);
0886                 } else if (hydroWhichHydro == 6) {
0887                     HydroCell_3D_ptr1 = (
0888                             &lattice_3D[position[0][ipy][ipeta][iptau]]);
0889                     HydroCell_3D_ptr2 = (
0890                             &lattice_3D[position[1][ipy][ipeta][iptau]]);
0891                     T += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->temperature
0892                                   + xfrac*HydroCell_3D_ptr2->temperature);
0893                     vx += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->vx
0894                                     + xfrac*HydroCell_3D_ptr2->vx);
0895                     vy += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->vy
0896                                     + xfrac*HydroCell_3D_ptr2->vy);
0897                     vz += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->vz
0898                                     + xfrac*HydroCell_3D_ptr2->vz);
0899                     pi00 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi00
0900                                     + xfrac*HydroCell_3D_ptr2->pi00);
0901                     pi01 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi01
0902                                     + xfrac*HydroCell_3D_ptr2->pi01);
0903                     pi02 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi02
0904                                     + xfrac*HydroCell_3D_ptr2->pi02);
0905                     pi03 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi03
0906                                     + xfrac*HydroCell_3D_ptr2->pi03);
0907                     pi11 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi11
0908                                     + xfrac*HydroCell_3D_ptr2->pi11);
0909                     pi12 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi12
0910                                     + xfrac*HydroCell_3D_ptr2->pi12);
0911                     pi13 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi13
0912                                     + xfrac*HydroCell_3D_ptr2->pi13);
0913                     pi22 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi22
0914                                     + xfrac*HydroCell_3D_ptr2->pi22);
0915                     pi23 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi23
0916                                     + xfrac*HydroCell_3D_ptr2->pi23);
0917                     pi33 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi33
0918                                     + xfrac*HydroCell_3D_ptr2->pi33);
0919                     bulkPi += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->bulkPi
0920                                     + xfrac*HydroCell_3D_ptr2->bulkPi);
0921                 } else if (hydroWhichHydro == 10 || hydroWhichHydro == 11) {
0922                     HydroCell_3D_ideal_ptr1 = (
0923                         &lattice_3D_ideal[idx_map_[position[0][ipy][ipeta][iptau]]]);
0924                     HydroCell_3D_ideal_ptr2 = (
0925                         &lattice_3D_ideal[idx_map_[position[1][ipy][ipeta][iptau]]]);
0926                     T += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->temperature
0927                                   + xfrac*HydroCell_3D_ideal_ptr2->temperature);
0928                     ed += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->ed
0929                                   + xfrac*HydroCell_3D_ideal_ptr2->ed);
0930                     p += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->pressure
0931                                   + xfrac*HydroCell_3D_ideal_ptr2->pressure);
0932                     ux += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->ux
0933                                     + xfrac*HydroCell_3D_ideal_ptr2->ux);
0934                     uy += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->uy
0935                                     + xfrac*HydroCell_3D_ideal_ptr2->uy);
0936                     uz += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->uz
0937                                    + xfrac*HydroCell_3D_ideal_ptr2->uz);
0938                 }
0939             }
0940         }
0941     }
0942 
0943     if (hydroWhichHydro == 10) {
0944         double ut = sqrt(1. + ux*ux + uy*uy + uz*uz);
0945         vx = ux/ut;
0946         vy = uy/ut;
0947         vz = uz/ut;
0948     } else if (hydroWhichHydro == 11) {
0949         double eta_local = 0.5*log((t + z)/(t - z));
0950         double sinh_eta = sinh(eta_local);
0951         double cosh_eta = cosh(eta_local);
0952         ueta = uz;
0953         double utau = sqrt(1. + ux*ux + uy*uy + ueta*ueta);
0954         uz = utau*sinh_eta + ueta*cosh_eta;
0955         double ut = utau*cosh_eta + ueta*sinh_eta;
0956         vx = ux/ut;
0957         vy = uy/ut;
0958         vz = uz/ut;
0959     } else if (hydroWhichHydro == 8 || hydroWhichHydro == 9) {
0960         double eta_local = 0.5*log((t + z)/(t - z));
0961         double sinh_eta = sinh(eta_local);
0962         double cosh_eta = cosh(eta_local);
0963         double utau = sqrt(1. + ux*ux + uy*uy + ueta*ueta);
0964         uz = utau*sinh_eta + ueta*cosh_eta;
0965         double ut = utau*cosh_eta + ueta*sinh_eta;
0966         vx = ux/ut;
0967         vy = uy/ut;
0968         vz = uz/ut;
0969     }
0970 
0971     info->temperature = T;
0972     info->vx = vx;
0973     info->vy = vy;
0974     info->vz = vz;
0975 
0976     if (hydroWhichHydro < 10) {
0977         info->ed = 1.0;                 // pi's are already divided by e+P
0978         info->sd = 0.0;
0979         info->pressure = 0.0;
0980     } else {
0981         info->ed = ed;
0982         info->sd = (ed + p)/(T + 1e-16);
0983         info->pressure = p;
0984     }
0985 
0986     info->pi[0][0] = pi00;
0987     info->pi[0][1] = pi01;
0988     info->pi[0][2] = pi02;
0989     info->pi[0][3] = pi03;
0990     info->pi[1][0] = pi01;
0991     info->pi[1][1] = pi11;
0992     info->pi[1][2] = pi12;
0993     info->pi[1][3] = pi13;
0994     info->pi[2][0] = pi02;
0995     info->pi[2][1] = pi12;
0996     info->pi[2][2] = pi22;
0997     info->pi[2][3] = pi23;
0998     info->pi[3][0] = pi03;
0999     info->pi[3][1] = pi13;
1000     info->pi[3][2] = pi23;
1001     info->pi[3][3] = pi33;
1002 
1003     info->bulkPi = bulkPi;
1004     return;
1005 }
1006 
1007 
1008 void Hydroinfo_MUSIC::output_temperature_evolution(string filename_base) {
1009     hydrofluidCell *hydroInfo = new hydrofluidCell;
1010     for (int i = 0; i < itaumax; i++) {
1011         double tau = hydroTau0 + i*hydroDtau;
1012         ostringstream filename;
1013         filename << filename_base << "_tau_" << tau << ".dat";
1014         ofstream temp_evo(filename.str().c_str());
1015         for (int ix = 0; ix < ixmax; ix++) {
1016             double x_local = -hydroXmax + ix*hydroDx;
1017             for (int iy = 0; iy < ixmax; iy++) {
1018                 double y_local = -hydroXmax + iy*hydroDx;
1019                 getHydroValues(x_local, y_local, 0.0, tau, hydroInfo);
1020                 double temp_local = hydroInfo->temperature;
1021                 temp_evo << scientific << setw(16) << setprecision(8)
1022                          << temp_local << "   ";
1023             }
1024             temp_evo << endl;
1025         }
1026         temp_evo.close();
1027     }
1028     delete hydroInfo;
1029 }
1030 
1031 
1032 void Hydroinfo_MUSIC::update_grid_info(
1033     double tau0, double tau_max, double dtau,
1034     double x_max, double dx, double eta_max, double deta) {
1035     hydroTau0 = tau0;
1036     hydroTauMax = tau_max;
1037     hydroDtau = dtau;
1038     hydroXmax = x_max;
1039     hydroDx = dx;
1040     hydro_eta_max = eta_max;
1041     hydroDeta = deta;
1042     if (hydroWhichHydro == 8) {
1043         itaumax = static_cast<int>((tau_max-tau0)/dtau+0.001);
1044         ixmax = static_cast<int>(2*x_max/dx+0.001);
1045         ietamax = static_cast<int>(2*eta_max/deta+0.001);
1046     }
1047     if (hydroWhichHydro == 6) {
1048         itaumax = static_cast<int>((tau_max-tau0)/dtau+0.001);
1049         ixmax = static_cast<int>(2*x_max/dx+0.001);
1050         ietamax = static_cast<int>(2*eta_max/deta+0.001);
1051     }
1052 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!