src/hydro/MusicWrapper.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 <string>
0021 #include <sstream>
0022 #include <vector>
0023 #include <memory>
0024
0025 #include "JetScapeLogger.h"
0026 #include "MusicWrapper.h"
0027
0028 using namespace Jetscape;
0029
0030 // Register the module with the base class
0031 RegisterJetScapeModule<MpiMusic> MpiMusic::reg("MUSIC");
0032
0033 MpiMusic::MpiMusic() {
0034   hydro_status = NOT_START;
0035   freezeout_temperature = 0.0;
0036   doCooperFrye = 0;
0037   flag_output_evo_to_file = 0;
0038   has_source_terms = false;
0039   SetId("MUSIC");
0040   hydro_source_terms_ptr =
0041       std::shared_ptr<HydroSourceJETSCAPE>(new HydroSourceJETSCAPE());
0042 }
0043
0044 MpiMusic::~MpiMusic() {}
0045
0046 void MpiMusic::InitializeHydro(Parameter parameter_list) {
0047   JSINFO << "Initialize MUSIC ...";
0048   VERBOSE(8);
0049
0050   string input_file = GetXMLElementText({"Hydro", "MUSIC", "MUSIC_input_file"});
0051   doCooperFrye =
0052       GetXMLElementInt({"Hydro", "MUSIC", "Perform_CooperFrye_Feezeout"});
0053
0054   music_hydro_ptr = std::unique_ptr<MUSIC>(new MUSIC(input_file));
0055
0056   // overwrite input options
0057   int echoLevel = GetXMLElementInt({"vlevel"});
0058   music_hydro_ptr->set_parameter("JSechoLevel", echoLevel);
0059
0060   flag_output_evo_to_file = (
0061       GetXMLElementInt({"Hydro", "MUSIC", "output_evolution_to_file"}));
0062   if (flag_output_evo_to_file == 1) {
0063     music_hydro_ptr->set_parameter("output_evolution_to_file", 2);
0064   }
0065   double tau_hydro = (
0066           GetXMLElementDouble({"Hydro", "MUSIC", "Initial_time_tau_0"}));
0067   music_hydro_ptr->set_parameter("Initial_time_tau_0", tau_hydro);
0068   int beastMode = GetXMLElementInt({"Hydro", "MUSIC", "beastMode"});
0069   music_hydro_ptr->set_parameter("beastMode", beastMode);
0070
0071   double eta_over_s =
0072       GetXMLElementDouble({"Hydro", "MUSIC", "shear_viscosity_eta_over_s"});
0073   if (eta_over_s > 1e-6) {
0074     music_hydro_ptr->set_parameter("Viscosity_Flag_Yes_1_No_0", 1);
0075     music_hydro_ptr->set_parameter("Include_Shear_Visc_Yes_1_No_0", 1);
0076     music_hydro_ptr->set_parameter("Shear_to_S_ratio", eta_over_s);
0077   } else if (eta_over_s >= 0.) {
0078     music_hydro_ptr->set_parameter("Viscosity_Flag_Yes_1_No_0", 0);
0079     music_hydro_ptr->set_parameter("Include_Shear_Visc_Yes_1_No_0", 0);
0080   } else {
0081     JSWARN << "The input shear viscosity is negative! eta/s = " << eta_over_s;
0082     exit(1);
0083   }
0084
0085   int flag_shear_Tdep = (
0086       GetXMLElementInt({"Hydro", "MUSIC", "T_dependent_Shear_to_S_ratio"}));
0087   if (flag_shear_Tdep > 0) {
0088     music_hydro_ptr->set_parameter("Viscosity_Flag_Yes_1_No_0", 1);
0089     music_hydro_ptr->set_parameter("T_dependent_Shear_to_S_ratio",
0090                                    flag_shear_Tdep);
0091     if (flag_shear_Tdep == 3) {
0092       double shear_kinkT = (
0093         GetXMLElementDouble({"Hydro", "MUSIC", "eta_over_s_T_kink_in_GeV"}));
0094       music_hydro_ptr->set_parameter("shear_viscosity_3_T_kink_in_GeV",
0095                                      shear_kinkT);
0096       double shear_lowTslope = (
0097         GetXMLElementDouble({"Hydro", "MUSIC",
0098                              "eta_over_s_low_T_slope_in_GeV"}));
0099       music_hydro_ptr->set_parameter("shear_viscosity_3_low_T_slope_in_GeV",
0100                                      shear_lowTslope);
0101       double shear_highTslope = (
0102         GetXMLElementDouble({"Hydro", "MUSIC",
0103                              "eta_over_s_high_T_slope_in_GeV"}));
0104       music_hydro_ptr->set_parameter("shear_viscosity_3_high_T_slope_in_GeV",
0105                                      shear_highTslope);
0106       double shear_kink = (
0107         GetXMLElementDouble({"Hydro", "MUSIC", "eta_over_s_at_kink"}));
0108       music_hydro_ptr->set_parameter("shear_viscosity_3_at_kink", shear_kink);
0109     }
0110   }
0111
0112   int flag_bulkvis = GetXMLElementInt(
0113           {"Hydro", "MUSIC", "temperature_dependent_bulk_viscosity"});
0114   if (flag_bulkvis != 0) {
0115     music_hydro_ptr->set_parameter("Include_Bulk_Visc_Yes_1_No_0", 1);
0116     music_hydro_ptr->set_parameter("T_dependent_Bulk_to_S_ratio",
0117                                    flag_bulkvis);
0118     if (flag_bulkvis == 3) {
0119         double bulk_max = GetXMLElementDouble(
0120               {"Hydro", "MUSIC", "zeta_over_s_max"});
0121         music_hydro_ptr->set_parameter("bulk_viscosity_3_max", bulk_max);
0122         double bulk_peakT = GetXMLElementDouble(
0123               {"Hydro", "MUSIC", "zeta_over_s_T_peak_in_GeV"});
0124         music_hydro_ptr->set_parameter("bulk_viscosity_3_T_peak_in_GeV",
0125                                        bulk_peakT);
0126         double bulk_width = GetXMLElementDouble(
0127               {"Hydro", "MUSIC", "zeta_over_s_width_in_GeV"});
0128         music_hydro_ptr->set_parameter("bulk_viscosity_3_width_in_GeV",
0129                                        bulk_width);
0130         double bulk_asy = GetXMLElementDouble(
0131               {"Hydro", "MUSIC", "zeta_over_s_lambda_asymm"});
0132         music_hydro_ptr->set_parameter("bulk_viscosity_3_lambda_asymm",
0133                                        bulk_asy);
0134     }
0135   }
0136
0137   int flag_secondorderTerms = GetXMLElementInt(
0138           {"Hydro", "MUSIC", "Include_second_order_terms"});
0139   if (flag_secondorderTerms == 1) {
0140     music_hydro_ptr->set_parameter("Include_second_order_terms", 1);
0141   }
0142
0143   freezeout_temperature =
0144       GetXMLElementDouble({"Hydro", "MUSIC", "freezeout_temperature"});
0145   if (freezeout_temperature > 0.05) {
0146     music_hydro_ptr->set_parameter("T_freeze", freezeout_temperature);
0147   } else {
0148     JSWARN << "The input freeze-out temperature is too low! T_frez = "
0149            << freezeout_temperature << " GeV!";
0150     exit(1);
0151   }
0152
0153   music_hydro_ptr->add_hydro_source_terms(hydro_source_terms_ptr);
0154 }
0155
0156 void MpiMusic::EvolveHydro() {
0157   VERBOSE(8);
0158   JSINFO << "Initialize density profiles in MUSIC ...";
0159   std::vector<double> entropy_density = ini->GetEntropyDensityDistribution();
0160   double dx = ini->GetXStep();
0161   double dz = ini->GetZStep();
0162   double z_max = ini->GetZMax();
0163   int nz = ini->GetZSize();
0164   if (pre_eq_ptr == nullptr) {
0165     JSWARN << "Missing the pre-equilibrium module ...";
0166   } else {
0167     music_hydro_ptr->initialize_hydro_from_jetscape_preequilibrium_vectors(
0168         dx, dz, z_max, nz, pre_eq_ptr->e_, pre_eq_ptr->P_,
0169         pre_eq_ptr->utau_, pre_eq_ptr->ux_, pre_eq_ptr->uy_, pre_eq_ptr->ueta_,
0170         pre_eq_ptr->pi00_, pre_eq_ptr->pi01_, pre_eq_ptr->pi02_,
0171         pre_eq_ptr->pi03_, pre_eq_ptr->pi11_, pre_eq_ptr->pi12_,
0172         pre_eq_ptr->pi13_, pre_eq_ptr->pi22_, pre_eq_ptr->pi23_,
0173         pre_eq_ptr->pi33_, pre_eq_ptr->bulk_Pi_);
0174   }
0175
0176   JSINFO << "initial density profile dx = " << dx << " fm";
0177   hydro_status = INITIALIZED;
0178   JSINFO << "number of source terms: "
0179          << hydro_source_terms_ptr->get_number_of_sources()
0180          << ", total E = " << hydro_source_terms_ptr->get_total_E_of_sources()
0181          << " GeV.";
0182
0183   has_source_terms = false;
0184   if (hydro_source_terms_ptr->get_number_of_sources() > 0) {
0185     has_source_terms = true;
0186   }
0187
0188   if (hydro_status == INITIALIZED) {
0189     JSINFO << "running MUSIC ...";
0190     music_hydro_ptr->run_hydro();
0191     hydro_status = FINISHED;
0192   }
0193
0194   if (flag_output_evo_to_file == 1) {
0195     if (!has_source_terms) {
0196       // only the first hydro without source term will be stored
0197       // in memory for jet energy loss calculations
0198       PassHydroEvolutionHistoryToFramework();
0199       JSINFO << "number of fluid cells received by the JETSCAPE: "
0200              << bulk_info.data.size();
0201     }
0202     music_hydro_ptr->clear_hydro_info_from_memory();
0203
0204     // add hydro_id to the hydro evolution filename
0205     std::ostringstream system_command;
0206     system_command << "mv evolution_all_xyeta.dat "
0207                    << "evolution_all_xyeta_" << GetId() << ".dat";
0208     system(system_command.str().c_str());
0209
0210     //std::vector<SurfaceCellInfo> surface_cells;
0211     //if (freezeout_temperature > 0.0) {
0212     //  FindAConstantTemperatureSurface(freezeout_temperature, surface_cells);
0213     //}
0214   }
0215
0216   collect_freeze_out_surface();
0217
0218   if (hydro_status == FINISHED && doCooperFrye == 1) {
0219     music_hydro_ptr->run_Cooper_Frye();
0220   }
0221 }
0222
0223 void MpiMusic::collect_freeze_out_surface() {
0224   system("rm surface.dat 2> /dev/null");
0225
0226   std::ostringstream surface_filename;
0227   surface_filename << "surface_" << GetId() << ".dat";
0228
0229   std::ostringstream system_command;
0230   system_command << "rm " << surface_filename.str() << " 2> /dev/null";
0231   system(system_command.str().c_str());
0232   system_command.str("");
0233   system_command.clear();
0234   system_command << "cat surface_eps* >> " << surface_filename.str();
0235   system(system_command.str().c_str());
0236   system_command.str("");
0237   system_command.clear();
0238
0239   system_command << "ln -s " << surface_filename.str() << " surface.dat";
0240   system(system_command.str().c_str());
0241   system_command.str("");
0242   system_command.clear();
0243   system("rm surface_eps* 2> /dev/null");
0244 }
0245
0246 void MpiMusic::SetHydroGridInfo() {
0247   bulk_info.neta = music_hydro_ptr->get_neta();
0248   bulk_info.ntau = music_hydro_ptr->get_ntau();
0249   bulk_info.nx = music_hydro_ptr->get_nx();
0250   bulk_info.ny = music_hydro_ptr->get_nx();
0251   bulk_info.tau_min = music_hydro_ptr->get_hydro_tau0();
0252   bulk_info.dtau = music_hydro_ptr->get_hydro_dtau();
0253   bulk_info.x_min = -music_hydro_ptr->get_hydro_x_max();
0254   bulk_info.dx = music_hydro_ptr->get_hydro_dx();
0255   bulk_info.y_min = -music_hydro_ptr->get_hydro_x_max();
0256   bulk_info.dy = music_hydro_ptr->get_hydro_dx();
0257   bulk_info.eta_min = -music_hydro_ptr->get_hydro_eta_max();
0258   bulk_info.deta = music_hydro_ptr->get_hydro_deta();
0259
0260   bulk_info.boost_invariant = music_hydro_ptr->is_boost_invariant();
0261 }
0262
0263 void MpiMusic::PassHydroEvolutionHistoryToFramework() {
0264   clear_up_evolution_data();
0265
0266   JSINFO << "Passing hydro evolution information to JETSCAPE ... ";
0267   auto number_of_cells = music_hydro_ptr->get_number_of_fluid_cells();
0268   JSINFO << "total number of fluid cells: " << number_of_cells;
0269
0270   SetHydroGridInfo();
0271
0272   fluidCell *fluidCell_ptr = new fluidCell;
0273   for (int i = 0; i < number_of_cells; i++) {
0274     std::unique_ptr<FluidCellInfo> fluid_cell_info_ptr(new FluidCellInfo);
0275     music_hydro_ptr->get_fluid_cell_with_index(i, fluidCell_ptr);
0276
0277     fluid_cell_info_ptr->energy_density = fluidCell_ptr->ed;
0278     fluid_cell_info_ptr->entropy_density = fluidCell_ptr->sd;
0279     fluid_cell_info_ptr->temperature = fluidCell_ptr->temperature;
0280     fluid_cell_info_ptr->pressure = fluidCell_ptr->pressure;
0281     fluid_cell_info_ptr->vx = fluidCell_ptr->vx;
0282     fluid_cell_info_ptr->vy = fluidCell_ptr->vy;
0283     fluid_cell_info_ptr->vz = fluidCell_ptr->vz;
0284     fluid_cell_info_ptr->mu_B = 0.0;
0285     fluid_cell_info_ptr->mu_C = 0.0;
0286     fluid_cell_info_ptr->mu_S = 0.0;
0287     fluid_cell_info_ptr->qgp_fraction = 0.0;
0288     for (int i = 0; i < 4; i++) {
0289       for (int j = 0; j < 4; j++) {
0290         fluid_cell_info_ptr->pi[i][j] = fluidCell_ptr->pi[i][j];
0291       }
0292     }
0293     fluid_cell_info_ptr->bulk_Pi = fluidCell_ptr->bulkPi;
0294     StoreHydroEvolutionHistory(fluid_cell_info_ptr);
0295   }
0296   delete fluidCell_ptr;
0297 }
0298
0299 void MpiMusic::GetHydroInfo(
0300     Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
0301     std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
0302   GetHydroInfo_JETSCAPE(t, x, y, z, fluid_cell_info_ptr);
0303   //GetHydroInfo_MUSIC(t, x, y, z, fluid_cell_info_ptr);
0304 }
0305
0306 void MpiMusic::GetHydroInfo_JETSCAPE(
0307     Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
0308     std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
0309   auto temp = bulk_info.get_tz(t, x, y, z);
0310   fluid_cell_info_ptr = std::unique_ptr<FluidCellInfo>(new FluidCellInfo(temp));
0311 }
0312
0313 void MpiMusic::GetHydroInfo_MUSIC(
0314     Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
0315     std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
0316   fluid_cell_info_ptr = Jetscape::make_unique<FluidCellInfo>();
0317   fluidCell *fluidCell_ptr = new fluidCell;
0318   music_hydro_ptr->get_hydro_info(x, y, z, t, fluidCell_ptr);
0319   fluid_cell_info_ptr->energy_density = fluidCell_ptr->ed;
0320   fluid_cell_info_ptr->entropy_density = fluidCell_ptr->sd;
0321   fluid_cell_info_ptr->temperature = fluidCell_ptr->temperature;
0322   fluid_cell_info_ptr->pressure = fluidCell_ptr->pressure;
0323   fluid_cell_info_ptr->vx = fluidCell_ptr->vx;
0324   fluid_cell_info_ptr->vy = fluidCell_ptr->vy;
0325   fluid_cell_info_ptr->vz = fluidCell_ptr->vz;
0326   fluid_cell_info_ptr->mu_B = 0.0;
0327   fluid_cell_info_ptr->mu_C = 0.0;
0328   fluid_cell_info_ptr->mu_S = 0.0;
0329   fluid_cell_info_ptr->qgp_fraction = 0.0;
0330
0331   for (int i = 0; i < 4; i++) {
0332     for (int j = 0; j < 4; j++) {
0333       fluid_cell_info_ptr->pi[i][j] = fluidCell_ptr->pi[i][j];
0334     }
0335   }
0336   fluid_cell_info_ptr->bulk_Pi = fluidCell_ptr->bulkPi;
0337   delete fluidCell_ptr;
0338 }