Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HydroFromFile.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HydroFromFile.cc
1 /*******************************************************************************
2  * Copyright (c) The JETSCAPE Collaboration, 2018
3  *
4  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
5  *
6  * For the list of contributors see AUTHORS.
7  *
8  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
9  *
10  * or via email to bugs.jetscape@gmail.com
11  *
12  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
13  * See COPYING for details.
14  ******************************************************************************/
15 
16 #include <stdio.h>
17 #include <sys/stat.h>
18 #include <MakeUniqueHelper.h>
19 
20 #include <cstring>
21 #include <sstream>
22 #include <cmath>
23 #include <iostream>
24 
25 #include "JetScapeLogger.h"
26 
27 #include "HydroFromFile.h"
28 
29 using namespace Jetscape;
30 
31 // Register the module with the base class
32 RegisterJetScapeModule<HydroFromFile> HydroFromFile::reg("HydroFromFile");
33 
34 HydroFromFile::HydroFromFile() {
35  hydro_status = NOT_START;
36  SetId("hydroFromFile");
37  PreEq_tau0_ = 0.;
38  PreEq_tauf_ = 0.;
39 }
40 
41 HydroFromFile::~HydroFromFile() { clean_hydro_event(); }
42 
44 void HydroFromFile::InitializeHydro(Parameter parameter_list) {
45  JSDEBUG << "Initialize hydro from file (Test) ...";
46  VERBOSE(8);
47 
48  string s = GetXMLElementText({"Hydro", "hydro_from_file", "name"});
49  JSDEBUG << s << " to be initialized ...";
50 
51  hydro_type_ = GetXMLElementInt({"Hydro", "hydro_from_file", "hydro_type"});
52  int boost_inv = GetXMLElementInt(
53  {"Hydro", "hydro_from_file", "boost_invariant_"});
54  if (boost_inv == 0) {
55  boost_invariant_ = false;
56  } else {
57  boost_invariant_ = true;
58  }
59  load_viscous_ =
60  GetXMLElementInt({"Hydro", "hydro_from_file", "load_viscous_info"});
61  nskip_tau_ =
62  GetXMLElementInt({"Hydro", "hydro_from_file", "read_hydro_every_ntau"});
63  T_c_ = GetXMLElementDouble({"Hydro", "hydro_from_file", "T_c"});
64  flag_read_in_multiple_hydro_ =
65  GetXMLElementInt({"Hydro", "hydro_from_file", "read_in_multiple_hydro"});
66 
67  hydro_event_idx_ = 0;
68 
69  if (hydro_type_ == 1) {
70 #ifdef USE_HDF5
71  hydroinfo_h5_ptr = new HydroinfoH5();
72 #else
73  JSWARN << " : hydro_type == 1 requires the hdf5 library~";
74  JSWARN << " : please check your inputs~";
75  exit(-1);
76 #endif
77  } else if (hydro_type_ == 2 || hydro_type_ == 3
78  || hydro_type_ == 4 || hydro_type_ == 5) {
79  hydroinfo_MUSIC_ptr = new Hydroinfo_MUSIC();
80  int verbose = GetXMLElementInt({"vlevel"});
81  hydroinfo_MUSIC_ptr->set_verbose(verbose);
82  } else if (hydro_type_ == 7) {
83  hydroinfo_MUSIC_ptr = new Hydroinfo_MUSIC();
84  hydroinfo_PreEq_ptr = new Hydroinfo_MUSIC();
85  int verbose = GetXMLElementInt({"vlevel"});
86  hydroinfo_MUSIC_ptr->set_verbose(verbose);
87  hydroinfo_PreEq_ptr->set_verbose(verbose);
88  }
89 
90  hydro_status = INITIALIZED;
91 }
92 
94 void HydroFromFile::read_in_hydro_event(string VISH_filename, int buffer_size,
95  int load_viscous) {
96  JSINFO << "read in a VISHNew hydro event from file " << VISH_filename;
97  if (hydro_type_ == 1) {
98 #ifdef USE_HDF5
99  hydroinfo_h5_ptr->readHydroinfoH5(VISH_filename, buffer_size, load_viscous);
100 #endif
101  }
102  hydro_status = FINISHED;
103 }
104 
106 void HydroFromFile::read_in_hydro_event(string MUSIC_input_file,
107  string MUSIC_hydro_ideal_file,
108  int nskip_tau) {
109  JSINFO << "read in a MUSIC hydro event from file " << MUSIC_hydro_ideal_file;
110  string hydro_shear_file = "";
111  string hydro_bulk_file = "";
112  int hydro_mode = 0;
113  if (hydro_type_ == 2) {
114  hydro_mode = 8;
115  } else if (hydro_type_ == 3) {
116  hydro_mode = 9;
117  } else if (hydro_type_ == 4) {
118  hydro_mode = 10;
119  } else if (hydro_type_ == 5) {
120  hydro_mode = 11;
121  }
122  hydroinfo_MUSIC_ptr->readHydroData(hydro_mode, nskip_tau, MUSIC_input_file,
123  MUSIC_hydro_ideal_file, hydro_shear_file,
124  hydro_bulk_file);
125  hydro_status = FINISHED;
126 }
127 
128 
130 void HydroFromFile::read_in_hydro_event(string MUSIC_input_file,
131  string PreEq_file,
132  string MUSIC_hydro_ideal_file,
133  int nskip_tau) {
134  int hydro_mode = 10;
135  if (hydro_type_ == 8) {
136  hydro_mode = 11; // 3D
137  }
138  string hydro_shear_file = "";
139  string hydro_bulk_file = "";
140  JSINFO << "read in a PreEq event from file " << PreEq_file;
141  hydroinfo_PreEq_ptr->readHydroData(hydro_mode, nskip_tau, MUSIC_input_file,
142  PreEq_file, hydro_shear_file,
143  hydro_bulk_file);
144  JSINFO << "read in a MUSIC hydro event from file " << MUSIC_hydro_ideal_file;
145  hydroinfo_MUSIC_ptr->readHydroData(hydro_mode, nskip_tau, MUSIC_input_file,
146  MUSIC_hydro_ideal_file, hydro_shear_file,
147  hydro_bulk_file);
148  hydro_status = FINISHED;
149 }
150 
151 void HydroFromFile::EvolveHydro() {
152  if (hydro_status == FINISHED) {
153  clean_hydro_event();
154  hydro_event_idx_ = ini->GetEventId();
155  }
156 
157  if (hydro_type_ == 1) {
158  string filename;
159  if (flag_read_in_multiple_hydro_ == 0) {
160  filename = GetXMLElementText({"Hydro", "hydro_from_file", "VISH_file"});
161  } else {
162  string folder =
163  GetXMLElementText({"Hydro", "hydro_from_file", "hydro_files_folder"});
164  std::ostringstream hydro_filename;
165  hydro_filename << folder << "/event-" << hydro_event_idx_
166  << "/JetData.h5";
167  filename = hydro_filename.str();
168  }
169 #ifdef USE_HDF5
170  read_in_hydro_event(filename, 500, load_viscous_);
171  hydro_tau_0 = hydroinfo_h5_ptr->getHydrogridTau0();
172  hydro_tau_max = hydroinfo_h5_ptr->getHydrogridTaumax();
173 #endif
174  hydro_status = FINISHED;
175  } else if (hydro_type_ < 6) {
176  string input_file;
177  string hydro_ideal_file;
178  if (flag_read_in_multiple_hydro_ == 0) {
179  input_file = GetXMLElementText(
180  {"Hydro", "hydro_from_file", "MUSIC_input_file"});
181  hydro_ideal_file = GetXMLElementText(
182  {"Hydro", "hydro_from_file", "MUSIC_file"});
183  } else {
184  string folder = GetXMLElementText(
185  {"Hydro", "hydro_from_file", "hydro_files_folder"});
186  std::ostringstream input_filename;
187  std::ostringstream hydro_filename;
188  input_filename << folder << "/event-" << hydro_event_idx_
189  << "/MUSIC_input";
190  hydro_filename << folder << "/event-" << hydro_event_idx_
191  << "/MUSIC_evo.dat";
192  input_file = input_filename.str();
193  hydro_ideal_file = hydro_filename.str();
194  }
195  read_in_hydro_event(input_file, hydro_ideal_file, nskip_tau_);
196  hydro_tau_0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
197  hydro_tau_max = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
198  } else if (hydro_type_ < 9) {
199  string input_file = "music";
200  string PreEq_file;
201  string hydro_ideal_file;
202  if (flag_read_in_multiple_hydro_ == 0) {
203  PreEq_file = GetXMLElementText(
204  {"Hydro", "hydro_from_file", "PreEq_file"});
205  hydro_ideal_file = GetXMLElementText(
206  {"Hydro", "hydro_from_file", "MUSIC_file"});
207  } else {
208  string folder = GetXMLElementText(
209  {"Hydro", "hydro_from_file", "hydro_files_folder"});
210  std::ostringstream preEq_filename;
211  std::ostringstream hydro_filename;
212  preEq_filename << folder << "/event-" << hydro_event_idx_
213  << "/PreEq_evo.dat";
214  hydro_filename << folder << "/event-" << hydro_event_idx_
215  << "/MUSIC_evo.dat";
216  PreEq_file = preEq_filename.str();
217  hydro_ideal_file = hydro_filename.str();
218  }
219  read_in_hydro_event(input_file, PreEq_file, hydro_ideal_file, 1);
220  PreEq_tau0_ = hydroinfo_PreEq_ptr->get_hydro_tau0();
221  PreEq_tauf_ = hydroinfo_PreEq_ptr->get_hydro_tau_max();
222  hydro_tau_0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
223  hydro_tau_max = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
224  if (std::abs(PreEq_tauf_ - hydro_tau_0) > 1e-6) {
225  JSWARN << __PRETTY_FUNCTION__
226  << "Preequilibrium medium end time is not the same as "
227  << "the hydro starting time! "
228  << "PreEq: tau_f = " << PreEq_tauf_ << " fm/c, "
229  << "hydro: tau_0 = " << hydro_tau_0 << " fm/c.";
230  }
231  } else {
232  JSWARN << __PRETTY_FUNCTION__
233  << "unrecognized hydro_type = " << hydro_type_;
234  exit(1);
235  }
236 }
237 
239 void HydroFromFile::clean_hydro_event() {
240  JSINFO << " clean up the loaded hydro event ...";
241  if (hydro_type_ == 1) {
242 #ifdef USE_HDF5
243  hydroinfo_h5_ptr->clean_hydro_event();
244 #endif
245  } else {
246  hydroinfo_MUSIC_ptr->clean_hydro_event();
247  }
248 
249  if (hydro_type_ == 7) {
250  hydroinfo_PreEq_ptr->clean_hydro_event();
251  }
252  hydro_status = NOT_START;
253 }
254 
257 void HydroFromFile::GetHydroInfo(
258  Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
259  std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
260  if (hydro_status != FINISHED) {
261  JSWARN << "Hydro not run yet ...";
262  exit(-1);
263  }
264 
265  double t_local = static_cast<double>(t);
266  double x_local = static_cast<double>(x);
267  double y_local = static_cast<double>(y);
268  double z_local = static_cast<double>(z);
269 
270  if (t_local < z_local) {
271  JSWARN << "Request a space-like point, t = " << t_local
272  << ", z = " << z_local;
273  exit(1);
274  }
275  double tau_local = sqrt(t * t - z * z);
276  double eta_local = 0.5 * log((t + z) / (t - z + 1e-15));
277  if (boost_invariant_) {
278  eta_local = 0.;
279  }
280 
281  // initialize the fluid cell pointer
282  hydrofluidCell *temp_fluid_cell_ptr = new hydrofluidCell;
283  if (hydro_type_ == 1) { // for OSU 2+1d hydro
284 #ifdef USE_HDF5
285  eta_local = 0.5 * log((t + z) / (t - z + 1e-15));
286  hydroinfo_h5_ptr->getHydroinfo(tau_local, x_local, y_local,
287  temp_fluid_cell_ptr);
288  // compute the flow velocity in the lab frame
289  double u0_perp =
290  (1. / sqrt(1. - temp_fluid_cell_ptr->vx * temp_fluid_cell_ptr->vx -
291  temp_fluid_cell_ptr->vy * temp_fluid_cell_ptr->vy));
292  double u0 = u0_perp * cosh(eta_local);
293  temp_fluid_cell_ptr->vx *= u0_perp / u0;
294  temp_fluid_cell_ptr->vy *= u0_perp / u0;
295  temp_fluid_cell_ptr->vz = z / (t + 1e-15);
296 #endif
297  } else if (hydro_type_ < 6) {
298  t_local = tau_local*cosh(eta_local);
299  z_local = tau_local*sinh(eta_local);
300  hydroinfo_MUSIC_ptr->getHydroValues(x_local, y_local, z_local, t_local,
301  temp_fluid_cell_ptr);
302  } else if (hydro_type_ < 9) {
303  t_local = tau_local*cosh(eta_local);
304  z_local = tau_local*sinh(eta_local);
305  if (tau_local < PreEq_tauf_) {
306  hydroinfo_PreEq_ptr->getHydroValues(x_local, y_local, z_local, t_local,
307  temp_fluid_cell_ptr);
308  } else {
309  hydroinfo_MUSIC_ptr->getHydroValues(x_local, y_local, z_local, t_local,
310  temp_fluid_cell_ptr);
311  }
312  }
313 
314  // assign all the quantites to JETSCAPE output
315  // thermodyanmic quantities
316  fluid_cell_info_ptr = make_unique<FluidCellInfo>();
317  fluid_cell_info_ptr->energy_density =
318  (static_cast<Jetscape::real>(temp_fluid_cell_ptr->ed));
319  fluid_cell_info_ptr->entropy_density =
320  (static_cast<Jetscape::real>(temp_fluid_cell_ptr->sd));
321  fluid_cell_info_ptr->temperature =
322  (static_cast<Jetscape::real>(temp_fluid_cell_ptr->temperature));
323  fluid_cell_info_ptr->pressure =
324  (static_cast<Jetscape::real>(temp_fluid_cell_ptr->pressure));
325  // QGP fraction
326  double qgp_fraction_local = 1.0;
327  if (temp_fluid_cell_ptr->temperature < T_c_) {
328  qgp_fraction_local = 0.0;
329  }
330  fluid_cell_info_ptr->qgp_fraction =
331  static_cast<Jetscape::real>(qgp_fraction_local);
332  // chemical potentials
333  fluid_cell_info_ptr->mu_B = 0.0;
334  fluid_cell_info_ptr->mu_C = 0.0;
335  fluid_cell_info_ptr->mu_S = 0.0;
336  // dynamical quantites
337  fluid_cell_info_ptr->vx =
338  static_cast<Jetscape::real>(temp_fluid_cell_ptr->vx);
339  fluid_cell_info_ptr->vy =
340  static_cast<Jetscape::real>(temp_fluid_cell_ptr->vy);
341  fluid_cell_info_ptr->vz =
342  static_cast<Jetscape::real>(temp_fluid_cell_ptr->vz);
343  for (int i = 0; i < 4; i++) {
344  for (int j = 0; j < 4; j++) {
345  fluid_cell_info_ptr->pi[i][j] =
346  (static_cast<Jetscape::real>(temp_fluid_cell_ptr->pi[i][j]));
347  }
348  }
349  fluid_cell_info_ptr->bulk_Pi =
350  (static_cast<Jetscape::real>(temp_fluid_cell_ptr->bulkPi));
351  delete temp_fluid_cell_ptr;
352 }
353 
354 double HydroFromFile::GetEventPlaneAngle() {
355  double v2 = 0.0;
356  double psi_2 = 0.0;
357  if (hydro_type_ == 1) {
358  std::ostringstream angle_filename;
359  string folder =
360  GetXMLElementText({"Hydro", "hydro_from_file", "hydro_files_folder"});
361  angle_filename << folder << "/event-" << hydro_event_idx_
362  << "/EventPlanesFrzout.dat";
363  std::ifstream inputfile(angle_filename.str().c_str());
364  string dummy;
365  std::getline(inputfile, dummy);
366  std::getline(inputfile, dummy);
367  std::getline(inputfile, dummy);
368  inputfile >> dummy >> v2 >> dummy >> psi_2;
369  inputfile.close();
370  }
371  return (psi_2);
372 }