src/framework/FluidDynamics.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 <iostream>
0018 #include <array>
0019 #include "FluidDynamics.h"
0020 #include "LinearInterpolation.h"
0021 #include "JetScapeSignalManager.h"
0022 #include "MakeUniqueHelper.h"
0023 #include "SurfaceFinder.h"
0024
0025 #define MAGENTA "\033[35m"
0026
0027 using namespace std;
0028
0029 namespace Jetscape {
0030
0031 FluidDynamics::FluidDynamics() {
0032   VERBOSE(8);
0033   eta = -99.99;
0034   boost_invariant_ = true;
0035   SetId("FluidDynamics");
0036 }
0037
0038 FluidDynamics::~FluidDynamics() {
0039   VERBOSE(8);
0040   disconnect_all();
0041 }
0042
0043 void FluidDynamics::Init() {
0044   JetScapeModuleBase::Init();
0045
0046   JSINFO << "Initialize FluidDynamics : " << GetId() << " ...";
0047
0048   VERBOSE(8);
0049   ini = JetScapeSignalManager::Instance()->GetInitialStatePointer().lock();
0050   if (!ini) {
0051     JSWARN << "No initialization module, "
0052            << "try: auto trento = make_shared<TrentoInitial>(); "
0053            << "jetscape->Add(trento);";
0054   }
0055
0056   pre_eq_ptr =
0057       JetScapeSignalManager::Instance()->GetPreEquilibriumPointer().lock();
0058   if (!pre_eq_ptr) {
0059     JSWARN << "No Pre-equilibrium module";
0060   }
0061
0062   InitializeHydro(parameter_list);
0063   InitTask();
0064
0065   JetScapeTask::InitTasks();
0066 }
0067
0068 void FluidDynamics::Exec() {
0069   VERBOSE(2) << "Run Hydro : " << GetId() << " ...";
0070   VERBOSE(8) << "Current Event #" << GetCurrentEvent();
0071
0072   if (ini) {
0073     VERBOSE(3) << "length of entropy density vector="
0074                << ini->GetEntropyDensityDistribution().size();
0075   }
0076
0077   EvolveHydro();
0078   JetScapeTask::ExecuteTasks();
0079 }
0080
0081 void FluidDynamics::Clear() {
0082   clear_up_evolution_data();
0083   if (!weak_ptr_is_uninitialized(liquefier_ptr)) {
0084     liquefier_ptr.lock()->Clear();
0085   }
0086 }
0087
0088 void FluidDynamics::CollectHeader(weak_ptr<JetScapeWriter> w) {
0089   auto f = w.lock();
0090   if (f) {
0091     auto &header = f->GetHeader();
0092     header.SetEventPlaneAngle(GetEventPlaneAngle());
0093   }
0094 }
0095
0096 void FluidDynamics::FindAConstantTemperatureSurface(
0097         Jetscape::real T_sw, std::vector<SurfaceCellInfo> &surface_cells) {
0098   std::unique_ptr<SurfaceFinder> surface_finder_ptr(
0099       new SurfaceFinder(T_sw, bulk_info));
0100   surface_finder_ptr->Find_full_hypersurface();
0101   surface_cells = surface_finder_ptr->get_surface_cells_vector();
0102   JSINFO << "number of surface cells: " << surface_cells.size();
0103 }
0104
0105 // this function returns the energy density [GeV] at a space time point
0106 // (time, x, y, z)
0107 Jetscape::real FluidDynamics::GetEnergyDensity(Jetscape::real time,
0108                                                Jetscape::real x,
0109                                                Jetscape::real y,
0110                                                Jetscape::real z) {
0111   std::unique_ptr<FluidCellInfo> fluid_cell_ptr;
0112   GetHydroInfo(time, x, y, z, fluid_cell_ptr);
0113   real energy_density = fluid_cell_ptr->energy_density;
0114   return (energy_density);
0115 }
0116
0117 // this function returns the entropy density [GeV] at a space time point
0118 // (time, x, y, z)
0119 Jetscape::real FluidDynamics::GetEntropyDensity(Jetscape::real time,
0120                                                 Jetscape::real x,
0121                                                 Jetscape::real y,
0122                                                 Jetscape::real z) {
0123   std::unique_ptr<FluidCellInfo> fluid_cell_ptr;
0124   GetHydroInfo(time, x, y, z, fluid_cell_ptr);
0125   real entropy_density = fluid_cell_ptr->entropy_density;
0126   return (entropy_density);
0127 }
0128
0129 // this function returns the temperature [GeV] at a space time point
0130 // (time, x, y, z)
0131 Jetscape::real FluidDynamics::GetTemperature(Jetscape::real time,
0132                                              Jetscape::real x, Jetscape::real y,
0133                                              Jetscape::real z) {
0134   std::unique_ptr<FluidCellInfo> fluid_cell_ptr;
0135   GetHydroInfo(time, x, y, z, fluid_cell_ptr);
0136   real temperature = fluid_cell_ptr->temperature;
0137   return (temperature);
0138 }
0139
0140 // this function returns the QGP fraction at a space time point
0141 // (time, x, y, z)
0142 Jetscape::real FluidDynamics::GetQgpFraction(Jetscape::real time,
0143                                              Jetscape::real x, Jetscape::real y,
0144                                              Jetscape::real z) {
0145   std::unique_ptr<FluidCellInfo> fluid_cell_ptr;
0146   GetHydroInfo(time, x, y, z, fluid_cell_ptr);
0147   real qgp_fraction = fluid_cell_ptr->qgp_fraction;
0148   return (qgp_fraction);
0149 }
0150
0151 void FluidDynamics::get_source_term(Jetscape::real tau, Jetscape::real x,
0152                                     Jetscape::real y, Jetscape::real eta,
0153                                     std::array<Jetscape::real, 4> jmu) const {
0154   liquefier_ptr.lock()->get_source(tau, x, y, eta, jmu);
0155 }
0156
0157 void FluidDynamics::PrintFluidCellInformation(
0158     FluidCellInfo *fluid_cell_info_ptr) {
0159   // this function print out the information of the fluid cell to the screen
0160   JSINFO << "=======================================================";
0161   JSINFO << "print out cell information:";
0162   JSINFO << "=======================================================";
0163   JSINFO << "energy density = " << fluid_cell_info_ptr->energy_density
0164          << " GeV/fm^3.";
0165   JSINFO << "entropy density = " << fluid_cell_info_ptr->entropy_density
0166          << " 1/fm^3.";
0167   JSINFO << "temperature = " << fluid_cell_info_ptr->temperature << " GeV.";
0168   JSINFO << "pressure = " << fluid_cell_info_ptr->pressure << " GeV/fm^3.";
0169   JSINFO << "QGP_fraction = " << fluid_cell_info_ptr->qgp_fraction;
0170   JSINFO << "mu_B = " << fluid_cell_info_ptr->mu_B << " GeV.";
0171   JSINFO << "mu_S = " << fluid_cell_info_ptr->mu_S << " GeV.";
0172   JSINFO << "mu_C = " << fluid_cell_info_ptr->mu_C << " GeV.";
0173   JSINFO << "vx = " << fluid_cell_info_ptr->vx;
0174   JSINFO << "vy = " << fluid_cell_info_ptr->vy;
0175   JSINFO << "vz = " << fluid_cell_info_ptr->vz;
0176   JSINFO << "shear viscous pi^{munu} (GeV/fm^3): ";
0177   for (int i = 0; i < 4; i++) {
0178     for (int j = 0; j < 4; j++) {
0179       JSINFO << fluid_cell_info_ptr->pi[i][j];
0180     }
0181   }
0182   JSINFO << "bulk_Pi = " << fluid_cell_info_ptr->bulk_Pi << " GeV/fm^3";
0183   JSINFO << "=======================================================";
0184 }
0185
0186 void FluidDynamics::UpdateEnergyDeposit(int t, double edop) {
0187   //sigslot::lock_block<multi_threaded_local> lock(this);
0188   JSDEBUG << MAGENTA << "Jet Signal received : " << t << " " << edop;
0189 }
0190
0191 } // end namespace Jetscape