src/hydro/HydroFromFile.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #include <stdio.h>
0017 #include <sys/stat.h>
0018 #include <MakeUniqueHelper.h>
0019
0020 #include <cstring>
0021 #include <sstream>
0022 #include <cmath>
0023 #include <iostream>
0024
0025 #include "JetScapeLogger.h"
0026
0027 #include "HydroFromFile.h"
0028
0029 using namespace Jetscape;
0030
0031 // Register the module with the base class
0032 RegisterJetScapeModule<HydroFromFile> HydroFromFile::reg("HydroFromFile");
0033
0034 HydroFromFile::HydroFromFile() {
0035   hydro_status = NOT_START;
0036   SetId("hydroFromFile");
0037   PreEq_tau0_ = 0.;
0038   PreEq_tauf_ = 0.;
0039 }
0040
0041 HydroFromFile::~HydroFromFile() { clean_hydro_event(); }
0042
0043 //! this function loads the hydro files
0044 void HydroFromFile::InitializeHydro(Parameter parameter_list) {
0045   JSDEBUG << "Initialize hydro from file (Test) ...";
0046   VERBOSE(8);
0047
0048   string s = GetXMLElementText({"Hydro", "hydro_from_file", "name"});
0049   JSDEBUG << s << " to be initialized ...";
0050
0051   hydro_type_ = GetXMLElementInt({"Hydro", "hydro_from_file", "hydro_type"});
0052   int boost_inv = GetXMLElementInt(
0053           {"Hydro", "hydro_from_file", "boost_invariant_"});
0054   if (boost_inv == 0) {
0055       boost_invariant_ = false;
0056   } else {
0057       boost_invariant_ = true;
0058   }
0059   load_viscous_ =
0060       GetXMLElementInt({"Hydro", "hydro_from_file", "load_viscous_info"});
0061   nskip_tau_ =
0062       GetXMLElementInt({"Hydro", "hydro_from_file", "read_hydro_every_ntau"});
0063   T_c_ = GetXMLElementDouble({"Hydro", "hydro_from_file", "T_c"});
0064   flag_read_in_multiple_hydro_ =
0065       GetXMLElementInt({"Hydro", "hydro_from_file", "read_in_multiple_hydro"});
0066
0067   hydro_event_idx_ = 0;
0068
0069   if (hydro_type_ == 1) {
0070 #ifdef USE_HDF5
0071     hydroinfo_h5_ptr = new HydroinfoH5();
0072 #else
0073     JSWARN << " : hydro_type == 1 requires the hdf5 library~";
0074     JSWARN << " : please check your inputs~";
0075     exit(-1);
0076 #endif
0077   } else if (hydro_type_ == 2 || hydro_type_ == 3
0078              || hydro_type_ == 4 || hydro_type_ == 5) {
0079     hydroinfo_MUSIC_ptr = new Hydroinfo_MUSIC();
0080     int verbose = GetXMLElementInt({"vlevel"});
0081     hydroinfo_MUSIC_ptr->set_verbose(verbose);
0082   } else if (hydro_type_ == 7) {
0083     hydroinfo_MUSIC_ptr = new Hydroinfo_MUSIC();
0084     hydroinfo_PreEq_ptr = new Hydroinfo_MUSIC();
0085     int verbose = GetXMLElementInt({"vlevel"});
0086     hydroinfo_MUSIC_ptr->set_verbose(verbose);
0087     hydroinfo_PreEq_ptr->set_verbose(verbose);
0088   }
0089
0090   hydro_status = INITIALIZED;
0091 }
0092
0093 //! This function load a VISHNew hydro event
0094 void HydroFromFile::read_in_hydro_event(string VISH_filename, int buffer_size,
0095                                         int load_viscous) {
0096   JSINFO << "read in a VISHNew hydro event from file " << VISH_filename;
0097   if (hydro_type_ == 1) {
0098 #ifdef USE_HDF5
0099     hydroinfo_h5_ptr->readHydroinfoH5(VISH_filename, buffer_size, load_viscous);
0100 #endif
0101   }
0102   hydro_status = FINISHED;
0103 }
0104
0105 //! This function load a MUSIC hydro event
0106 void HydroFromFile::read_in_hydro_event(string MUSIC_input_file,
0107                                         string MUSIC_hydro_ideal_file,
0108                                         int nskip_tau) {
0109   JSINFO << "read in a MUSIC hydro event from file " << MUSIC_hydro_ideal_file;
0110   string hydro_shear_file = "";
0111   string hydro_bulk_file = "";
0112   int hydro_mode = 0;
0113   if (hydro_type_ == 2) {
0114     hydro_mode = 8;
0115   } else if (hydro_type_ == 3) {
0116     hydro_mode = 9;
0117   } else if (hydro_type_ == 4) {
0118     hydro_mode = 10;
0119   } else if (hydro_type_ == 5) {
0120     hydro_mode = 11;
0121   }
0122   hydroinfo_MUSIC_ptr->readHydroData(hydro_mode, nskip_tau, MUSIC_input_file,
0123                                      MUSIC_hydro_ideal_file, hydro_shear_file,
0124                                      hydro_bulk_file);
0125   hydro_status = FINISHED;
0126 }
0127
0128
0129   //! This function load a PreEq event and a MUSIC hydro event
0130 void HydroFromFile::read_in_hydro_event(string MUSIC_input_file,
0131                                         string PreEq_file,
0132                                         string MUSIC_hydro_ideal_file,
0133                                         int nskip_tau) {
0134   int hydro_mode = 10;
0135   if (hydro_type_ == 8) {
0136     hydro_mode = 11;  // 3D
0137   }
0138   string hydro_shear_file = "";
0139   string hydro_bulk_file = "";
0140   JSINFO << "read in a PreEq event from file " << PreEq_file;
0141   hydroinfo_PreEq_ptr->readHydroData(hydro_mode, nskip_tau, MUSIC_input_file,
0142                                      PreEq_file, hydro_shear_file,
0143                                      hydro_bulk_file);
0144   JSINFO << "read in a MUSIC hydro event from file " << MUSIC_hydro_ideal_file;
0145   hydroinfo_MUSIC_ptr->readHydroData(hydro_mode, nskip_tau, MUSIC_input_file,
0146                                      MUSIC_hydro_ideal_file, hydro_shear_file,
0147                                      hydro_bulk_file);
0148   hydro_status = FINISHED;
0149 }
0150
0151 void HydroFromFile::EvolveHydro() {
0152   if (hydro_status == FINISHED) {
0153     clean_hydro_event();
0154     hydro_event_idx_ = ini->GetEventId();
0155   }
0156
0157   if (hydro_type_ == 1) {
0158     string filename;
0159     if (flag_read_in_multiple_hydro_ == 0) {
0160       filename = GetXMLElementText({"Hydro", "hydro_from_file", "VISH_file"});
0161     } else {
0162       string folder =
0163           GetXMLElementText({"Hydro", "hydro_from_file", "hydro_files_folder"});
0164       std::ostringstream hydro_filename;
0165       hydro_filename << folder << "/event-" << hydro_event_idx_
0166                      << "/JetData.h5";
0167       filename = hydro_filename.str();
0168     }
0169 #ifdef USE_HDF5
0170     read_in_hydro_event(filename, 500, load_viscous_);
0171     hydro_tau_0 = hydroinfo_h5_ptr->getHydrogridTau0();
0172     hydro_tau_max = hydroinfo_h5_ptr->getHydrogridTaumax();
0173 #endif
0174     hydro_status = FINISHED;
0175   } else if (hydro_type_ < 6) {
0176     string input_file;
0177     string hydro_ideal_file;
0178     if (flag_read_in_multiple_hydro_ == 0) {
0179       input_file = GetXMLElementText(
0180               {"Hydro", "hydro_from_file", "MUSIC_input_file"});
0181       hydro_ideal_file = GetXMLElementText(
0182               {"Hydro", "hydro_from_file", "MUSIC_file"});
0183     } else {
0184       string folder = GetXMLElementText(
0185               {"Hydro", "hydro_from_file", "hydro_files_folder"});
0186       std::ostringstream input_filename;
0187       std::ostringstream hydro_filename;
0188       input_filename << folder << "/event-" << hydro_event_idx_
0189                      << "/MUSIC_input";
0190       hydro_filename << folder << "/event-" << hydro_event_idx_
0191                      << "/MUSIC_evo.dat";
0192       input_file = input_filename.str();
0193       hydro_ideal_file = hydro_filename.str();
0194     }
0195     read_in_hydro_event(input_file, hydro_ideal_file, nskip_tau_);
0196     hydro_tau_0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
0197     hydro_tau_max = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
0198   } else if (hydro_type_ < 9) {
0199     string input_file = "music";
0200     string PreEq_file;
0201     string hydro_ideal_file;
0202     if (flag_read_in_multiple_hydro_ == 0) {
0203       PreEq_file = GetXMLElementText(
0204               {"Hydro", "hydro_from_file", "PreEq_file"});
0205       hydro_ideal_file = GetXMLElementText(
0206               {"Hydro", "hydro_from_file", "MUSIC_file"});
0207     } else {
0208       string folder = GetXMLElementText(
0209               {"Hydro", "hydro_from_file", "hydro_files_folder"});
0210       std::ostringstream preEq_filename;
0211       std::ostringstream hydro_filename;
0212       preEq_filename << folder << "/event-" << hydro_event_idx_
0213                      << "/PreEq_evo.dat";
0214       hydro_filename << folder << "/event-" << hydro_event_idx_
0215                      << "/MUSIC_evo.dat";
0216       PreEq_file = preEq_filename.str();
0217       hydro_ideal_file = hydro_filename.str();
0218     }
0219     read_in_hydro_event(input_file, PreEq_file, hydro_ideal_file, 1);
0220     PreEq_tau0_ = hydroinfo_PreEq_ptr->get_hydro_tau0();
0221     PreEq_tauf_ = hydroinfo_PreEq_ptr->get_hydro_tau_max();
0222     hydro_tau_0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
0223     hydro_tau_max = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
0224     if (std::abs(PreEq_tauf_ - hydro_tau_0) > 1e-6) {
0225         JSWARN << __PRETTY_FUNCTION__
0226                << "Preequilibrium medium end time is not the same as "
0227                << "the hydro starting time! "
0228                << "PreEq: tau_f = " << PreEq_tauf_ << " fm/c, "
0229                << "hydro: tau_0 = " << hydro_tau_0 << " fm/c.";
0230     }
0231   } else {
0232     JSWARN << __PRETTY_FUNCTION__
0233            << "unrecognized hydro_type = " << hydro_type_;
0234     exit(1);
0235   }
0236 }
0237
0238 //! clean up hydro event
0239 void HydroFromFile::clean_hydro_event() {
0240   JSINFO << " clean up the loaded hydro event ...";
0241   if (hydro_type_ == 1) {
0242 #ifdef USE_HDF5
0243     hydroinfo_h5_ptr->clean_hydro_event();
0244 #endif
0245   } else {
0246     hydroinfo_MUSIC_ptr->clean_hydro_event();
0247   }
0248
0249   if (hydro_type_ == 7) {
0250     hydroinfo_PreEq_ptr->clean_hydro_event();
0251   }
0252   hydro_status = NOT_START;
0253 }
0254
0255 //! this function returns the thermodynamic and dynamical information at
0256 //! the given space-time point
0257 void HydroFromFile::GetHydroInfo(
0258     Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
0259     std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
0260   if (hydro_status != FINISHED) {
0261     JSWARN << "Hydro not run yet ...";
0262     exit(-1);
0263   }
0264
0265   double t_local = static_cast<double>(t);
0266   double x_local = static_cast<double>(x);
0267   double y_local = static_cast<double>(y);
0268   double z_local = static_cast<double>(z);
0269
0270   if (t_local < z_local) {
0271       JSWARN << "Request a space-like point, t = " << t_local
0272              << ", z = " << z_local;
0273       exit(1);
0274   }
0275   double tau_local = sqrt(t * t - z * z);
0276   double eta_local = 0.5 * log((t + z) / (t - z + 1e-15));
0277   if (boost_invariant_) {
0278       eta_local = 0.;
0279   }
0280
0281   // initialize the fluid cell pointer
0282   hydrofluidCell *temp_fluid_cell_ptr = new hydrofluidCell;
0283   if (hydro_type_ == 1) { // for OSU 2+1d hydro
0284 #ifdef USE_HDF5
0285     eta_local = 0.5 * log((t + z) / (t - z + 1e-15));
0286     hydroinfo_h5_ptr->getHydroinfo(tau_local, x_local, y_local,
0287                                    temp_fluid_cell_ptr);
0288     // compute the flow velocity in the lab frame
0289     double u0_perp =
0290         (1. / sqrt(1. - temp_fluid_cell_ptr->vx * temp_fluid_cell_ptr->vx -
0291                    temp_fluid_cell_ptr->vy * temp_fluid_cell_ptr->vy));
0292     double u0 = u0_perp * cosh(eta_local);
0293     temp_fluid_cell_ptr->vx *= u0_perp / u0;
0294     temp_fluid_cell_ptr->vy *= u0_perp / u0;
0295     temp_fluid_cell_ptr->vz = z / (t + 1e-15);
0296 #endif
0297   } else if (hydro_type_ < 6) {
0298     t_local = tau_local*cosh(eta_local);
0299     z_local = tau_local*sinh(eta_local);
0300     hydroinfo_MUSIC_ptr->getHydroValues(x_local, y_local, z_local, t_local,
0301                                         temp_fluid_cell_ptr);
0302   } else if (hydro_type_ < 9) {
0303     t_local = tau_local*cosh(eta_local);
0304     z_local = tau_local*sinh(eta_local);
0305     if (tau_local < PreEq_tauf_) {
0306         hydroinfo_PreEq_ptr->getHydroValues(x_local, y_local, z_local, t_local,
0307                                             temp_fluid_cell_ptr);
0308     } else {
0309         hydroinfo_MUSIC_ptr->getHydroValues(x_local, y_local, z_local, t_local,
0310                                             temp_fluid_cell_ptr);
0311     }
0312   }
0313
0314   // assign all the quantites to JETSCAPE output
0315   // thermodyanmic quantities
0316   fluid_cell_info_ptr = make_unique<FluidCellInfo>();
0317   fluid_cell_info_ptr->energy_density =
0318       (static_cast<Jetscape::real>(temp_fluid_cell_ptr->ed));
0319   fluid_cell_info_ptr->entropy_density =
0320       (static_cast<Jetscape::real>(temp_fluid_cell_ptr->sd));
0321   fluid_cell_info_ptr->temperature =
0322       (static_cast<Jetscape::real>(temp_fluid_cell_ptr->temperature));
0323   fluid_cell_info_ptr->pressure =
0324       (static_cast<Jetscape::real>(temp_fluid_cell_ptr->pressure));
0325   // QGP fraction
0326   double qgp_fraction_local = 1.0;
0327   if (temp_fluid_cell_ptr->temperature < T_c_) {
0328     qgp_fraction_local = 0.0;
0329   }
0330   fluid_cell_info_ptr->qgp_fraction =
0331       static_cast<Jetscape::real>(qgp_fraction_local);
0332   // chemical potentials
0333   fluid_cell_info_ptr->mu_B = 0.0;
0334   fluid_cell_info_ptr->mu_C = 0.0;
0335   fluid_cell_info_ptr->mu_S = 0.0;
0336   // dynamical quantites
0337   fluid_cell_info_ptr->vx =
0338       static_cast<Jetscape::real>(temp_fluid_cell_ptr->vx);
0339   fluid_cell_info_ptr->vy =
0340       static_cast<Jetscape::real>(temp_fluid_cell_ptr->vy);
0341   fluid_cell_info_ptr->vz =
0342       static_cast<Jetscape::real>(temp_fluid_cell_ptr->vz);
0343   for (int i = 0; i < 4; i++) {
0344     for (int j = 0; j < 4; j++) {
0345       fluid_cell_info_ptr->pi[i][j] =
0346           (static_cast<Jetscape::real>(temp_fluid_cell_ptr->pi[i][j]));
0347     }
0348   }
0349   fluid_cell_info_ptr->bulk_Pi =
0350       (static_cast<Jetscape::real>(temp_fluid_cell_ptr->bulkPi));
0351   delete temp_fluid_cell_ptr;
0352 }
0353
0354 double HydroFromFile::GetEventPlaneAngle() {
0355   double v2 = 0.0;
0356   double psi_2 = 0.0;
0357   if (hydro_type_ == 1) {
0358     std::ostringstream angle_filename;
0359     string folder =
0360         GetXMLElementText({"Hydro", "hydro_from_file", "hydro_files_folder"});
0361     angle_filename << folder << "/event-" << hydro_event_idx_
0362                    << "/EventPlanesFrzout.dat";
0363     std::ifstream inputfile(angle_filename.str().c_str());
0364     string dummy;
0365     std::getline(inputfile, dummy);
0366     std::getline(inputfile, dummy);
0367     std::getline(inputfile, dummy);
0368     inputfile >> dummy >> v2 >> dummy >> psi_2;
0369     inputfile.close();
0370   }
0371   return (psi_2);
0372 }