src/framework/FluidEvolutionHistory.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 // This is a general basic class for hydrodynamics
0016
0017 #include <string>
0018 #include <MakeUniqueHelper.h>
0019 #include "FluidEvolutionHistory.h"
0020 #include "FluidCellInfo.h"
0021 #include "LinearInterpolation.h"
0022 #include "JetScapeLogger.h"
0023
0024 namespace Jetscape {
0025
0026 // convert the string type entry name to enum type EntryNames
0027 EntryName ResolveEntryName(std::string input) {
0028   static const std::map<std::string, EntryName> optionStrings = {
0029       {"energy_density", ENTRY_ENERGY_DENSITY},
0030       {"entropy_density", ENTRY_ENTROPY_DENSITY},
0031       {"temperature", ENTRY_TEMPERATURE},
0032       {"pressure", ENTRY_PRESSURE},
0033       {"qgp_fraction", ENTRY_QGP_FRACTION},
0034       {"mu_b", ENTRY_MU_B},
0035       {"mu_c", ENTRY_MU_C},
0036       {"mu_s", ENTRY_MU_S},
0037       {"vx", ENTRY_VX},
0038       {"vy", ENTRY_VY},
0039       {"vz", ENTRY_VZ},
0040       {"pi00", ENTRY_PI00},
0041       {"pi01", ENTRY_PI01},
0042       {"pi02", ENTRY_PI02},
0043       {"pi03", ENTRY_PI03},
0044       {"pi11", ENTRY_PI11},
0045       {"pi12", ENTRY_PI12},
0046       {"pi13", ENTRY_PI13},
0047       {"pi22", ENTRY_PI22},
0048       {"pi23", ENTRY_PI23},
0049       {"pi33", ENTRY_PI33},
0050       {"bulk_pi", ENTRY_BULK_PI},
0051   };
0052   auto itr = optionStrings.find(input);
0053   if (itr != optionStrings.end()) {
0054     return itr->second;
0055   } else {
0056     return ENTRY_INVALID;
0057   }
0058 }
0059
0060 // It checks whether a space-time point (tau, x, y, eta) is inside evolution
0061 // history or outside.
0062 int EvolutionHistory::CheckInRange(Jetscape::real tau, Jetscape::real x,
0063                                    Jetscape::real y, Jetscape::real eta) const {
0064   int status = 1;
0065   if (tau < tau_min || tau > TauMax()) {
0066     std::string warn_message =
0067         ("tau=" + std::to_string(tau) + " is not in range [" +
0068          std::to_string(tau_min) + "," + std::to_string(TauMax()) + "]");
0069     //throw InvalidSpaceTimeRange(warn_message);
0070     //JSWARN << warn_message;
0071     status = 0;
0072   }
0073   if (x < x_min || x > XMax()) {
0074     std::string warn_message =
0075         ("x=" + std::to_string(x) + " is not in range [" +
0076          std::to_string(x_min) + "," + std::to_string(XMax()) + "]");
0077     //throw InvalidSpaceTimeRange(warn_message);
0078     //JSWARN << warn_message;
0079     status = 0;
0080   }
0081   if (y < y_min || y > YMax()) {
0082     std::string warn_message =
0083         ("y=" + std::to_string(y) + " is not in range [" +
0084          std::to_string(y_min) + "," + std::to_string(YMax()) + "]");
0085     //throw InvalidSpaceTimeRange(warn_message);
0086     //JSWARN << warn_message;
0087     status = 0;
0088   }
0089   if (!boost_invariant) {
0090     if (eta < eta_min || eta > EtaMax()) {
0091       std::string warn_message =
0092           ("eta=" + std::to_string(eta) + " is not in range [" +
0093            std::to_string(eta_min) + "," + std::to_string(EtaMax()) + "]");
0094       //throw InvalidSpaceTimeRange(warn_message);
0095       //JSWARN << warn_message;
0096       status = 0;
0097     }
0098   }
0099   return (status);
0100 }
0101
0102 /** Construct evolution history given the bulk_data and the data_info */
0103 void EvolutionHistory::FromVector(const std::vector<float> &data_,
0104                                   const std::vector<std::string> &data_info_,
0105                                   float tau_min_, float dtau_, float x_min_,
0106                                   float dx_, int nx_, float y_min_, float dy_,
0107                                   int ny_, float eta_min_, float deta_,
0108                                   int neta_, bool tau_eta_is_tz_) {
0109   data_vector = data_;
0110   data_info = data_info_;
0111   tau_min = tau_min_;
0112   x_min = x_min_;
0113   y_min = y_min_;
0114   eta_min = eta_min_;
0115   dtau = dtau_;
0116   dx = dx_;
0117   dy = dy_;
0118   deta = deta_;
0119   nx = nx_;
0120   ny = ny_;
0121   neta = neta_;
0122   tau_eta_is_tz = tau_eta_is_tz_;
0123   ntau = data_.size() / (data_info_.size() * nx * ny * neta);
0124 }
0125
0126 /* This function will read the sparse data stored in data_ with associated
0127  * information data_info_ into to FluidCellInfo object */
0128 FluidCellInfo EvolutionHistory::GetFluidCell(int id_tau, int id_x, int id_y,
0129                                              int id_eta) const {
0130   int entries_per_record = data_info.size();
0131   int id_eta_corrected = id_eta;
0132   // set id_eta=0 if hydro is in 2+1D mode
0133   if (neta == 0 || neta == 1) {
0134     id_eta_corrected = 0;
0135   }
0136
0137   int record_starting_id = CellIndex(id_tau, id_x, id_y, id_eta_corrected);
0138
0139   // if data_vector and data_info are not used to construct evolution history
0140   // then the data should have the format of vector<FluidCellInfo>.
0141   if (entries_per_record == 0) {
0142     return data.at(record_starting_id);
0143   }
0144   // otherwise construct the fluid cell info from data_vector and data_info
0145   auto fluid_cell_ptr = make_unique<FluidCellInfo>();
0146
0147   record_starting_id *= entries_per_record;
0148   for (int i = 0; i < entries_per_record; i++) {
0149     auto entry_name = ResolveEntryName(data_info.at(i));
0150     auto entry_data = data_vector.at(record_starting_id + i);
0151     switch (entry_name) {
0152     case ENTRY_ENERGY_DENSITY:
0153       fluid_cell_ptr->energy_density = entry_data;
0154       break;
0155     case ENTRY_ENTROPY_DENSITY:
0156       fluid_cell_ptr->entropy_density = entry_data;
0157       break;
0158     case ENTRY_TEMPERATURE:
0159       fluid_cell_ptr->temperature = entry_data;
0160       break;
0161     case ENTRY_PRESSURE:
0162       fluid_cell_ptr->pressure = entry_data;
0163       break;
0164     case ENTRY_QGP_FRACTION:
0165       fluid_cell_ptr->qgp_fraction = entry_data;
0166       break;
0167     case ENTRY_MU_B:
0168       fluid_cell_ptr->mu_B = entry_data;
0169       break;
0170     case ENTRY_MU_C:
0171       fluid_cell_ptr->mu_C = entry_data;
0172       break;
0173     case ENTRY_MU_S:
0174       fluid_cell_ptr->mu_S = entry_data;
0175       break;
0176     case ENTRY_VX:
0177       fluid_cell_ptr->vx = entry_data;
0178       break;
0179     case ENTRY_VY:
0180       fluid_cell_ptr->vy = entry_data;
0181       break;
0182     case ENTRY_VZ:
0183       fluid_cell_ptr->vz = entry_data;
0184       break;
0185     case ENTRY_PI00:
0186       fluid_cell_ptr->pi[0][0] = entry_data;
0187       break;
0188     case ENTRY_PI01:
0189       fluid_cell_ptr->pi[0][1] = entry_data;
0190       fluid_cell_ptr->pi[1][0] = entry_data;
0191       break;
0192     case ENTRY_PI02:
0193       fluid_cell_ptr->pi[0][2] = entry_data;
0194       fluid_cell_ptr->pi[2][0] = entry_data;
0195       break;
0196     case ENTRY_PI03:
0197       fluid_cell_ptr->pi[0][3] = entry_data;
0198       fluid_cell_ptr->pi[3][0] = entry_data;
0199       break;
0200     case ENTRY_PI11:
0201       fluid_cell_ptr->pi[1][1] = entry_data;
0202       break;
0203     case ENTRY_PI12:
0204       fluid_cell_ptr->pi[1][2] = entry_data;
0205       fluid_cell_ptr->pi[2][1] = entry_data;
0206       break;
0207     case ENTRY_PI13:
0208       fluid_cell_ptr->pi[1][3] = entry_data;
0209       fluid_cell_ptr->pi[3][1] = entry_data;
0210       break;
0211     case ENTRY_PI22:
0212       fluid_cell_ptr->pi[2][2] = entry_data;
0213       break;
0214     case ENTRY_PI23:
0215       fluid_cell_ptr->pi[2][3] = entry_data;
0216       fluid_cell_ptr->pi[3][2] = entry_data;
0217       break;
0218     case ENTRY_PI33:
0219       fluid_cell_ptr->pi[3][3] = entry_data;
0220       break;
0221     case ENTRY_BULK_PI:
0222       fluid_cell_ptr->bulk_Pi = entry_data;
0223       break;
0224     default:
0225       JSWARN << "The entry name in data_info_ must be one of the \
0226                         energy_density, entropy_density, temperature, pressure, qgp_fraction, \
0227                         mu_b, mu_c, mu_s, vx, vy, vz, pi00, pi01, pi02, pi03, pi11, pi12, \
0228                         pi13, pi22, pi23, pi33, bulk_pi";
0229       break;
0230     }
0231   }
0232
0233   return *fluid_cell_ptr;
0234 }
0235
0236 /** For one given time step id_tau,
0237  * get FluidCellInfo at spatial point (x, y, eta)*/
0238 FluidCellInfo EvolutionHistory::GetAtTimeStep(int id_tau, Jetscape::real x,
0239                                               Jetscape::real y,
0240                                               Jetscape::real eta) const {
0241   int id_x = GetIdX(x);
0242   int id_y = GetIdY(y);
0243   int id_eta = 0;
0244   if (!boost_invariant)
0245     id_eta = GetIdEta(eta);
0246
0247   auto c000 = GetFluidCell(id_tau, id_x, id_y, id_eta);
0248   auto c001 = GetFluidCell(id_tau, id_x, id_y, id_eta + 1);
0249   auto c010 = GetFluidCell(id_tau, id_x, id_y + 1, id_eta);
0250   auto c011 = GetFluidCell(id_tau, id_x, id_y + 1, id_eta + 1);
0251   auto c100 = GetFluidCell(id_tau, id_x + 1, id_y, id_eta);
0252   auto c101 = GetFluidCell(id_tau, id_x + 1, id_y, id_eta + 1);
0253   auto c110 = GetFluidCell(id_tau, id_x + 1, id_y + 1, id_eta);
0254   auto c111 = GetFluidCell(id_tau, id_x + 1, id_y + 1, id_eta + 1);
0255   real x0 = XCoord(id_x);
0256   real x1 = XCoord(id_x + 1);
0257   real y0 = YCoord(id_y);
0258   real y1 = YCoord(id_y + 1);
0259   real eta0 = EtaCoord(id_eta);
0260   auto eta1 = 0.0;
0261   if (!boost_invariant)
0262     eta1 = EtaCoord(id_eta + 1);
0263
0264   return TrilinearInt(x0, x1, y0, y1, eta0, eta1, c000, c001, c010, c011, c100,
0265                       c101, c110, c111, x, y, eta);
0266 }
0267
0268 // do interpolation along time direction; we may also need high order
0269 // interpolation functions
0270 FluidCellInfo EvolutionHistory::get(Jetscape::real tau, Jetscape::real x,
0271                                     Jetscape::real y,
0272                                     Jetscape::real eta) const {
0273   int status = CheckInRange(tau, x, y, eta);
0274   if (status == 0) {
0275     FluidCellInfo zero_cell;
0276     return (zero_cell);
0277   }
0278   int id_tau = GetIdTau(tau);
0279   auto tau0 = TauCoord(id_tau);
0280   auto tau1 = TauCoord(id_tau + 1);
0281   auto bulk0 = GetAtTimeStep(id_tau, x, y, eta);
0282   auto bulk1 = GetAtTimeStep(id_tau + 1, x, y, eta);
0283   return (LinearInt(tau0, tau1, bulk0, bulk1, tau));
0284 }
0285
0286 FluidCellInfo EvolutionHistory::get_tz(Jetscape::real t, Jetscape::real x,
0287                                        Jetscape::real y,
0288                                        Jetscape::real z) const {
0289   Jetscape::real tau = 0.0;
0290   Jetscape::real eta = 0.0;
0291   if (t * t > z * z) {
0292     tau = sqrt(t * t - z * z);
0293     eta = 0.5 * log((t + z) / (t - z));
0294   } else {
0295     VERBOSE(4) << "the quest point is outside the light cone! "
0296                << "t = " << t << ", z = " << z;
0297   }
0298   return (get(tau, x, y, eta));
0299 }
0300
0301 } // end namespace Jetscape