Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FluidEvolutionHistory.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FluidEvolutionHistory.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 // This is a general basic class for hydrodynamics
16 
17 #include <string>
18 #include <MakeUniqueHelper.h>
19 #include "FluidEvolutionHistory.h"
20 #include "FluidCellInfo.h"
21 #include "LinearInterpolation.h"
22 #include "JetScapeLogger.h"
23 
24 namespace Jetscape {
25 
26 // convert the string type entry name to enum type EntryNames
27 EntryName ResolveEntryName(std::string input) {
28  static const std::map<std::string, EntryName> optionStrings = {
29  {"energy_density", ENTRY_ENERGY_DENSITY},
30  {"entropy_density", ENTRY_ENTROPY_DENSITY},
31  {"temperature", ENTRY_TEMPERATURE},
32  {"pressure", ENTRY_PRESSURE},
33  {"qgp_fraction", ENTRY_QGP_FRACTION},
34  {"mu_b", ENTRY_MU_B},
35  {"mu_c", ENTRY_MU_C},
36  {"mu_s", ENTRY_MU_S},
37  {"vx", ENTRY_VX},
38  {"vy", ENTRY_VY},
39  {"vz", ENTRY_VZ},
40  {"pi00", ENTRY_PI00},
41  {"pi01", ENTRY_PI01},
42  {"pi02", ENTRY_PI02},
43  {"pi03", ENTRY_PI03},
44  {"pi11", ENTRY_PI11},
45  {"pi12", ENTRY_PI12},
46  {"pi13", ENTRY_PI13},
47  {"pi22", ENTRY_PI22},
48  {"pi23", ENTRY_PI23},
49  {"pi33", ENTRY_PI33},
50  {"bulk_pi", ENTRY_BULK_PI},
51  };
52  auto itr = optionStrings.find(input);
53  if (itr != optionStrings.end()) {
54  return itr->second;
55  } else {
56  return ENTRY_INVALID;
57  }
58 }
59 
60 // It checks whether a space-time point (tau, x, y, eta) is inside evolution
61 // history or outside.
62 int EvolutionHistory::CheckInRange(Jetscape::real tau, Jetscape::real x,
63  Jetscape::real y, Jetscape::real eta) const {
64  int status = 1;
65  if (tau < tau_min || tau > TauMax()) {
66  std::string warn_message =
67  ("tau=" + std::to_string(tau) + " is not in range [" +
68  std::to_string(tau_min) + "," + std::to_string(TauMax()) + "]");
69  //throw InvalidSpaceTimeRange(warn_message);
70  //JSWARN << warn_message;
71  status = 0;
72  }
73  if (x < x_min || x > XMax()) {
74  std::string warn_message =
75  ("x=" + std::to_string(x) + " is not in range [" +
76  std::to_string(x_min) + "," + std::to_string(XMax()) + "]");
77  //throw InvalidSpaceTimeRange(warn_message);
78  //JSWARN << warn_message;
79  status = 0;
80  }
81  if (y < y_min || y > YMax()) {
82  std::string warn_message =
83  ("y=" + std::to_string(y) + " is not in range [" +
84  std::to_string(y_min) + "," + std::to_string(YMax()) + "]");
85  //throw InvalidSpaceTimeRange(warn_message);
86  //JSWARN << warn_message;
87  status = 0;
88  }
89  if (!boost_invariant) {
90  if (eta < eta_min || eta > EtaMax()) {
91  std::string warn_message =
92  ("eta=" + std::to_string(eta) + " is not in range [" +
93  std::to_string(eta_min) + "," + std::to_string(EtaMax()) + "]");
94  //throw InvalidSpaceTimeRange(warn_message);
95  //JSWARN << warn_message;
96  status = 0;
97  }
98  }
99  return (status);
100 }
101 
103 void EvolutionHistory::FromVector(const std::vector<float> &data_,
104  const std::vector<std::string> &data_info_,
105  float tau_min_, float dtau_, float x_min_,
106  float dx_, int nx_, float y_min_, float dy_,
107  int ny_, float eta_min_, float deta_,
108  int neta_, bool tau_eta_is_tz_) {
109  data_vector = data_;
110  data_info = data_info_;
111  tau_min = tau_min_;
112  x_min = x_min_;
113  y_min = y_min_;
114  eta_min = eta_min_;
115  dtau = dtau_;
116  dx = dx_;
117  dy = dy_;
118  deta = deta_;
119  nx = nx_;
120  ny = ny_;
121  neta = neta_;
122  tau_eta_is_tz = tau_eta_is_tz_;
123  ntau = data_.size() / (data_info_.size() * nx * ny * neta);
124 }
125 
126 /* This function will read the sparse data stored in data_ with associated
127  * information data_info_ into to FluidCellInfo object */
128 FluidCellInfo EvolutionHistory::GetFluidCell(int id_tau, int id_x, int id_y,
129  int id_eta) const {
130  int entries_per_record = data_info.size();
131  int id_eta_corrected = id_eta;
132  // set id_eta=0 if hydro is in 2+1D mode
133  if (neta == 0 || neta == 1) {
134  id_eta_corrected = 0;
135  }
136 
137  int record_starting_id = CellIndex(id_tau, id_x, id_y, id_eta_corrected);
138 
139  // if data_vector and data_info are not used to construct evolution history
140  // then the data should have the format of vector<FluidCellInfo>.
141  if (entries_per_record == 0) {
142  return data.at(record_starting_id);
143  }
144  // otherwise construct the fluid cell info from data_vector and data_info
145  auto fluid_cell_ptr = make_unique<FluidCellInfo>();
146 
147  record_starting_id *= entries_per_record;
148  for (int i = 0; i < entries_per_record; i++) {
149  auto entry_name = ResolveEntryName(data_info.at(i));
150  auto entry_data = data_vector.at(record_starting_id + i);
151  switch (entry_name) {
152  case ENTRY_ENERGY_DENSITY:
153  fluid_cell_ptr->energy_density = entry_data;
154  break;
155  case ENTRY_ENTROPY_DENSITY:
156  fluid_cell_ptr->entropy_density = entry_data;
157  break;
158  case ENTRY_TEMPERATURE:
159  fluid_cell_ptr->temperature = entry_data;
160  break;
161  case ENTRY_PRESSURE:
162  fluid_cell_ptr->pressure = entry_data;
163  break;
164  case ENTRY_QGP_FRACTION:
165  fluid_cell_ptr->qgp_fraction = entry_data;
166  break;
167  case ENTRY_MU_B:
168  fluid_cell_ptr->mu_B = entry_data;
169  break;
170  case ENTRY_MU_C:
171  fluid_cell_ptr->mu_C = entry_data;
172  break;
173  case ENTRY_MU_S:
174  fluid_cell_ptr->mu_S = entry_data;
175  break;
176  case ENTRY_VX:
177  fluid_cell_ptr->vx = entry_data;
178  break;
179  case ENTRY_VY:
180  fluid_cell_ptr->vy = entry_data;
181  break;
182  case ENTRY_VZ:
183  fluid_cell_ptr->vz = entry_data;
184  break;
185  case ENTRY_PI00:
186  fluid_cell_ptr->pi[0][0] = entry_data;
187  break;
188  case ENTRY_PI01:
189  fluid_cell_ptr->pi[0][1] = entry_data;
190  fluid_cell_ptr->pi[1][0] = entry_data;
191  break;
192  case ENTRY_PI02:
193  fluid_cell_ptr->pi[0][2] = entry_data;
194  fluid_cell_ptr->pi[2][0] = entry_data;
195  break;
196  case ENTRY_PI03:
197  fluid_cell_ptr->pi[0][3] = entry_data;
198  fluid_cell_ptr->pi[3][0] = entry_data;
199  break;
200  case ENTRY_PI11:
201  fluid_cell_ptr->pi[1][1] = entry_data;
202  break;
203  case ENTRY_PI12:
204  fluid_cell_ptr->pi[1][2] = entry_data;
205  fluid_cell_ptr->pi[2][1] = entry_data;
206  break;
207  case ENTRY_PI13:
208  fluid_cell_ptr->pi[1][3] = entry_data;
209  fluid_cell_ptr->pi[3][1] = entry_data;
210  break;
211  case ENTRY_PI22:
212  fluid_cell_ptr->pi[2][2] = entry_data;
213  break;
214  case ENTRY_PI23:
215  fluid_cell_ptr->pi[2][3] = entry_data;
216  fluid_cell_ptr->pi[3][2] = entry_data;
217  break;
218  case ENTRY_PI33:
219  fluid_cell_ptr->pi[3][3] = entry_data;
220  break;
221  case ENTRY_BULK_PI:
222  fluid_cell_ptr->bulk_Pi = entry_data;
223  break;
224  default:
225  JSWARN << "The entry name in data_info_ must be one of the \
226  energy_density, entropy_density, temperature, pressure, qgp_fraction, \
227  mu_b, mu_c, mu_s, vx, vy, vz, pi00, pi01, pi02, pi03, pi11, pi12, \
228  pi13, pi22, pi23, pi33, bulk_pi";
229  break;
230  }
231  }
232 
233  return *fluid_cell_ptr;
234 }
235 
238 FluidCellInfo EvolutionHistory::GetAtTimeStep(int id_tau, Jetscape::real x,
239  Jetscape::real y,
240  Jetscape::real eta) const {
241  int id_x = GetIdX(x);
242  int id_y = GetIdY(y);
243  int id_eta = 0;
244  if (!boost_invariant)
245  id_eta = GetIdEta(eta);
246 
247  auto c000 = GetFluidCell(id_tau, id_x, id_y, id_eta);
248  auto c001 = GetFluidCell(id_tau, id_x, id_y, id_eta + 1);
249  auto c010 = GetFluidCell(id_tau, id_x, id_y + 1, id_eta);
250  auto c011 = GetFluidCell(id_tau, id_x, id_y + 1, id_eta + 1);
251  auto c100 = GetFluidCell(id_tau, id_x + 1, id_y, id_eta);
252  auto c101 = GetFluidCell(id_tau, id_x + 1, id_y, id_eta + 1);
253  auto c110 = GetFluidCell(id_tau, id_x + 1, id_y + 1, id_eta);
254  auto c111 = GetFluidCell(id_tau, id_x + 1, id_y + 1, id_eta + 1);
255  real x0 = XCoord(id_x);
256  real x1 = XCoord(id_x + 1);
257  real y0 = YCoord(id_y);
258  real y1 = YCoord(id_y + 1);
259  real eta0 = EtaCoord(id_eta);
260  auto eta1 = 0.0;
261  if (!boost_invariant)
262  eta1 = EtaCoord(id_eta + 1);
263 
264  return TrilinearInt(x0, x1, y0, y1, eta0, eta1, c000, c001, c010, c011, c100,
265  c101, c110, c111, x, y, eta);
266 }
267 
268 // do interpolation along time direction; we may also need high order
269 // interpolation functions
270 FluidCellInfo EvolutionHistory::get(Jetscape::real tau, Jetscape::real x,
271  Jetscape::real y,
272  Jetscape::real eta) const {
273  int status = CheckInRange(tau, x, y, eta);
274  if (status == 0) {
275  FluidCellInfo zero_cell;
276  return (zero_cell);
277  }
278  int id_tau = GetIdTau(tau);
279  auto tau0 = TauCoord(id_tau);
280  auto tau1 = TauCoord(id_tau + 1);
281  auto bulk0 = GetAtTimeStep(id_tau, x, y, eta);
282  auto bulk1 = GetAtTimeStep(id_tau + 1, x, y, eta);
283  return (LinearInt(tau0, tau1, bulk0, bulk1, tau));
284 }
285 
286 FluidCellInfo EvolutionHistory::get_tz(Jetscape::real t, Jetscape::real x,
287  Jetscape::real y,
288  Jetscape::real z) const {
289  Jetscape::real tau = 0.0;
290  Jetscape::real eta = 0.0;
291  if (t * t > z * z) {
292  tau = sqrt(t * t - z * z);
293  eta = 0.5 * log((t + z) / (t - z));
294  } else {
295  VERBOSE(4) << "the quest point is outside the light cone! "
296  << "t = " << t << ", z = " << z;
297  }
298  return (get(tau, x, y, eta));
299 }
300 
301 } // end namespace Jetscape