Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MusicWrapper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MusicWrapper.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 <string>
21 #include <sstream>
22 #include <vector>
23 #include <memory>
24 
25 #include "JetScapeLogger.h"
26 #include "MusicWrapper.h"
27 
28 using namespace Jetscape;
29 
30 // Register the module with the base class
31 RegisterJetScapeModule<MpiMusic> MpiMusic::reg("MUSIC");
32 
33 MpiMusic::MpiMusic() {
34  hydro_status = NOT_START;
35  freezeout_temperature = 0.0;
36  doCooperFrye = 0;
37  flag_output_evo_to_file = 0;
38  has_source_terms = false;
39  SetId("MUSIC");
40  hydro_source_terms_ptr =
41  std::shared_ptr<HydroSourceJETSCAPE>(new HydroSourceJETSCAPE());
42 }
43 
44 MpiMusic::~MpiMusic() {}
45 
46 void MpiMusic::InitializeHydro(Parameter parameter_list) {
47  JSINFO << "Initialize MUSIC ...";
48  VERBOSE(8);
49 
50  string input_file = GetXMLElementText({"Hydro", "MUSIC", "MUSIC_input_file"});
51  doCooperFrye =
52  GetXMLElementInt({"Hydro", "MUSIC", "Perform_CooperFrye_Feezeout"});
53 
54  music_hydro_ptr = std::unique_ptr<MUSIC>(new MUSIC(input_file));
55 
56  // overwrite input options
57  int echoLevel = GetXMLElementInt({"vlevel"});
58  music_hydro_ptr->set_parameter("JSechoLevel", echoLevel);
59 
60  flag_output_evo_to_file = (
61  GetXMLElementInt({"Hydro", "MUSIC", "output_evolution_to_file"}));
62  if (flag_output_evo_to_file == 1) {
63  music_hydro_ptr->set_parameter("output_evolution_to_file", 2);
64  }
65  double tau_hydro = (
66  GetXMLElementDouble({"Hydro", "MUSIC", "Initial_time_tau_0"}));
67  music_hydro_ptr->set_parameter("Initial_time_tau_0", tau_hydro);
68  int beastMode = GetXMLElementInt({"Hydro", "MUSIC", "beastMode"});
69  music_hydro_ptr->set_parameter("beastMode", beastMode);
70 
71  double eta_over_s =
72  GetXMLElementDouble({"Hydro", "MUSIC", "shear_viscosity_eta_over_s"});
73  if (eta_over_s > 1e-6) {
74  music_hydro_ptr->set_parameter("Viscosity_Flag_Yes_1_No_0", 1);
75  music_hydro_ptr->set_parameter("Include_Shear_Visc_Yes_1_No_0", 1);
76  music_hydro_ptr->set_parameter("Shear_to_S_ratio", eta_over_s);
77  } else if (eta_over_s >= 0.) {
78  music_hydro_ptr->set_parameter("Viscosity_Flag_Yes_1_No_0", 0);
79  music_hydro_ptr->set_parameter("Include_Shear_Visc_Yes_1_No_0", 0);
80  } else {
81  JSWARN << "The input shear viscosity is negative! eta/s = " << eta_over_s;
82  exit(1);
83  }
84 
85  int flag_shear_Tdep = (
86  GetXMLElementInt({"Hydro", "MUSIC", "T_dependent_Shear_to_S_ratio"}));
87  if (flag_shear_Tdep > 0) {
88  music_hydro_ptr->set_parameter("Viscosity_Flag_Yes_1_No_0", 1);
89  music_hydro_ptr->set_parameter("T_dependent_Shear_to_S_ratio",
90  flag_shear_Tdep);
91  if (flag_shear_Tdep == 3) {
92  double shear_kinkT = (
93  GetXMLElementDouble({"Hydro", "MUSIC", "eta_over_s_T_kink_in_GeV"}));
94  music_hydro_ptr->set_parameter("shear_viscosity_3_T_kink_in_GeV",
95  shear_kinkT);
96  double shear_lowTslope = (
97  GetXMLElementDouble({"Hydro", "MUSIC",
98  "eta_over_s_low_T_slope_in_GeV"}));
99  music_hydro_ptr->set_parameter("shear_viscosity_3_low_T_slope_in_GeV",
100  shear_lowTslope);
101  double shear_highTslope = (
102  GetXMLElementDouble({"Hydro", "MUSIC",
103  "eta_over_s_high_T_slope_in_GeV"}));
104  music_hydro_ptr->set_parameter("shear_viscosity_3_high_T_slope_in_GeV",
105  shear_highTslope);
106  double shear_kink = (
107  GetXMLElementDouble({"Hydro", "MUSIC", "eta_over_s_at_kink"}));
108  music_hydro_ptr->set_parameter("shear_viscosity_3_at_kink", shear_kink);
109  }
110  }
111 
112  int flag_bulkvis = GetXMLElementInt(
113  {"Hydro", "MUSIC", "temperature_dependent_bulk_viscosity"});
114  if (flag_bulkvis != 0) {
115  music_hydro_ptr->set_parameter("Include_Bulk_Visc_Yes_1_No_0", 1);
116  music_hydro_ptr->set_parameter("T_dependent_Bulk_to_S_ratio",
117  flag_bulkvis);
118  if (flag_bulkvis == 3) {
119  double bulk_max = GetXMLElementDouble(
120  {"Hydro", "MUSIC", "zeta_over_s_max"});
121  music_hydro_ptr->set_parameter("bulk_viscosity_3_max", bulk_max);
122  double bulk_peakT = GetXMLElementDouble(
123  {"Hydro", "MUSIC", "zeta_over_s_T_peak_in_GeV"});
124  music_hydro_ptr->set_parameter("bulk_viscosity_3_T_peak_in_GeV",
125  bulk_peakT);
126  double bulk_width = GetXMLElementDouble(
127  {"Hydro", "MUSIC", "zeta_over_s_width_in_GeV"});
128  music_hydro_ptr->set_parameter("bulk_viscosity_3_width_in_GeV",
129  bulk_width);
130  double bulk_asy = GetXMLElementDouble(
131  {"Hydro", "MUSIC", "zeta_over_s_lambda_asymm"});
132  music_hydro_ptr->set_parameter("bulk_viscosity_3_lambda_asymm",
133  bulk_asy);
134  }
135  }
136 
137  int flag_secondorderTerms = GetXMLElementInt(
138  {"Hydro", "MUSIC", "Include_second_order_terms"});
139  if (flag_secondorderTerms == 1) {
140  music_hydro_ptr->set_parameter("Include_second_order_terms", 1);
141  }
142 
143  freezeout_temperature =
144  GetXMLElementDouble({"Hydro", "MUSIC", "freezeout_temperature"});
145  if (freezeout_temperature > 0.05) {
146  music_hydro_ptr->set_parameter("T_freeze", freezeout_temperature);
147  } else {
148  JSWARN << "The input freeze-out temperature is too low! T_frez = "
149  << freezeout_temperature << " GeV!";
150  exit(1);
151  }
152 
153  music_hydro_ptr->add_hydro_source_terms(hydro_source_terms_ptr);
154 }
155 
156 void MpiMusic::EvolveHydro() {
157  VERBOSE(8);
158  JSINFO << "Initialize density profiles in MUSIC ...";
159  std::vector<double> entropy_density = ini->GetEntropyDensityDistribution();
160  double dx = ini->GetXStep();
161  double dz = ini->GetZStep();
162  double z_max = ini->GetZMax();
163  int nz = ini->GetZSize();
164  if (pre_eq_ptr == nullptr) {
165  JSWARN << "Missing the pre-equilibrium module ...";
166  } else {
167  music_hydro_ptr->initialize_hydro_from_jetscape_preequilibrium_vectors(
168  dx, dz, z_max, nz, pre_eq_ptr->e_, pre_eq_ptr->P_,
169  pre_eq_ptr->utau_, pre_eq_ptr->ux_, pre_eq_ptr->uy_, pre_eq_ptr->ueta_,
170  pre_eq_ptr->pi00_, pre_eq_ptr->pi01_, pre_eq_ptr->pi02_,
171  pre_eq_ptr->pi03_, pre_eq_ptr->pi11_, pre_eq_ptr->pi12_,
172  pre_eq_ptr->pi13_, pre_eq_ptr->pi22_, pre_eq_ptr->pi23_,
173  pre_eq_ptr->pi33_, pre_eq_ptr->bulk_Pi_);
174  }
175 
176  JSINFO << "initial density profile dx = " << dx << " fm";
177  hydro_status = INITIALIZED;
178  JSINFO << "number of source terms: "
179  << hydro_source_terms_ptr->get_number_of_sources()
180  << ", total E = " << hydro_source_terms_ptr->get_total_E_of_sources()
181  << " GeV.";
182 
183  has_source_terms = false;
184  if (hydro_source_terms_ptr->get_number_of_sources() > 0) {
185  has_source_terms = true;
186  }
187 
188  if (hydro_status == INITIALIZED) {
189  JSINFO << "running MUSIC ...";
190  music_hydro_ptr->run_hydro();
191  hydro_status = FINISHED;
192  }
193 
194  if (flag_output_evo_to_file == 1) {
195  if (!has_source_terms) {
196  // only the first hydro without source term will be stored
197  // in memory for jet energy loss calculations
198  PassHydroEvolutionHistoryToFramework();
199  JSINFO << "number of fluid cells received by the JETSCAPE: "
200  << bulk_info.data.size();
201  }
202  music_hydro_ptr->clear_hydro_info_from_memory();
203 
204  // add hydro_id to the hydro evolution filename
205  std::ostringstream system_command;
206  system_command << "mv evolution_all_xyeta.dat "
207  << "evolution_all_xyeta_" << GetId() << ".dat";
208  system(system_command.str().c_str());
209 
210  //std::vector<SurfaceCellInfo> surface_cells;
211  //if (freezeout_temperature > 0.0) {
212  // FindAConstantTemperatureSurface(freezeout_temperature, surface_cells);
213  //}
214  }
215 
216  collect_freeze_out_surface();
217 
218  if (hydro_status == FINISHED && doCooperFrye == 1) {
219  music_hydro_ptr->run_Cooper_Frye();
220  }
221 }
222 
223 void MpiMusic::collect_freeze_out_surface() {
224  system("rm surface.dat 2> /dev/null");
225 
226  std::ostringstream surface_filename;
227  surface_filename << "surface_" << GetId() << ".dat";
228 
229  std::ostringstream system_command;
230  system_command << "rm " << surface_filename.str() << " 2> /dev/null";
231  system(system_command.str().c_str());
232  system_command.str("");
233  system_command.clear();
234  system_command << "cat surface_eps* >> " << surface_filename.str();
235  system(system_command.str().c_str());
236  system_command.str("");
237  system_command.clear();
238 
239  system_command << "ln -s " << surface_filename.str() << " surface.dat";
240  system(system_command.str().c_str());
241  system_command.str("");
242  system_command.clear();
243  system("rm surface_eps* 2> /dev/null");
244 }
245 
246 void MpiMusic::SetHydroGridInfo() {
247  bulk_info.neta = music_hydro_ptr->get_neta();
248  bulk_info.ntau = music_hydro_ptr->get_ntau();
249  bulk_info.nx = music_hydro_ptr->get_nx();
250  bulk_info.ny = music_hydro_ptr->get_nx();
251  bulk_info.tau_min = music_hydro_ptr->get_hydro_tau0();
252  bulk_info.dtau = music_hydro_ptr->get_hydro_dtau();
253  bulk_info.x_min = -music_hydro_ptr->get_hydro_x_max();
254  bulk_info.dx = music_hydro_ptr->get_hydro_dx();
255  bulk_info.y_min = -music_hydro_ptr->get_hydro_x_max();
256  bulk_info.dy = music_hydro_ptr->get_hydro_dx();
257  bulk_info.eta_min = -music_hydro_ptr->get_hydro_eta_max();
258  bulk_info.deta = music_hydro_ptr->get_hydro_deta();
259 
260  bulk_info.boost_invariant = music_hydro_ptr->is_boost_invariant();
261 }
262 
263 void MpiMusic::PassHydroEvolutionHistoryToFramework() {
264  clear_up_evolution_data();
265 
266  JSINFO << "Passing hydro evolution information to JETSCAPE ... ";
267  auto number_of_cells = music_hydro_ptr->get_number_of_fluid_cells();
268  JSINFO << "total number of fluid cells: " << number_of_cells;
269 
270  SetHydroGridInfo();
271 
272  fluidCell *fluidCell_ptr = new fluidCell;
273  for (int i = 0; i < number_of_cells; i++) {
274  std::unique_ptr<FluidCellInfo> fluid_cell_info_ptr(new FluidCellInfo);
275  music_hydro_ptr->get_fluid_cell_with_index(i, fluidCell_ptr);
276 
277  fluid_cell_info_ptr->energy_density = fluidCell_ptr->ed;
278  fluid_cell_info_ptr->entropy_density = fluidCell_ptr->sd;
279  fluid_cell_info_ptr->temperature = fluidCell_ptr->temperature;
280  fluid_cell_info_ptr->pressure = fluidCell_ptr->pressure;
281  fluid_cell_info_ptr->vx = fluidCell_ptr->vx;
282  fluid_cell_info_ptr->vy = fluidCell_ptr->vy;
283  fluid_cell_info_ptr->vz = fluidCell_ptr->vz;
284  fluid_cell_info_ptr->mu_B = 0.0;
285  fluid_cell_info_ptr->mu_C = 0.0;
286  fluid_cell_info_ptr->mu_S = 0.0;
287  fluid_cell_info_ptr->qgp_fraction = 0.0;
288  for (int i = 0; i < 4; i++) {
289  for (int j = 0; j < 4; j++) {
290  fluid_cell_info_ptr->pi[i][j] = fluidCell_ptr->pi[i][j];
291  }
292  }
293  fluid_cell_info_ptr->bulk_Pi = fluidCell_ptr->bulkPi;
294  StoreHydroEvolutionHistory(fluid_cell_info_ptr);
295  }
296  delete fluidCell_ptr;
297 }
298 
299 void MpiMusic::GetHydroInfo(
300  Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
301  std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
302  GetHydroInfo_JETSCAPE(t, x, y, z, fluid_cell_info_ptr);
303  //GetHydroInfo_MUSIC(t, x, y, z, fluid_cell_info_ptr);
304 }
305 
306 void MpiMusic::GetHydroInfo_JETSCAPE(
307  Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
308  std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
309  auto temp = bulk_info.get_tz(t, x, y, z);
310  fluid_cell_info_ptr = std::unique_ptr<FluidCellInfo>(new FluidCellInfo(temp));
311 }
312 
313 void MpiMusic::GetHydroInfo_MUSIC(
314  Jetscape::real t, Jetscape::real x, Jetscape::real y, Jetscape::real z,
315  std::unique_ptr<FluidCellInfo> &fluid_cell_info_ptr) {
316  fluid_cell_info_ptr = Jetscape::make_unique<FluidCellInfo>();
317  fluidCell *fluidCell_ptr = new fluidCell;
318  music_hydro_ptr->get_hydro_info(x, y, z, t, fluidCell_ptr);
319  fluid_cell_info_ptr->energy_density = fluidCell_ptr->ed;
320  fluid_cell_info_ptr->entropy_density = fluidCell_ptr->sd;
321  fluid_cell_info_ptr->temperature = fluidCell_ptr->temperature;
322  fluid_cell_info_ptr->pressure = fluidCell_ptr->pressure;
323  fluid_cell_info_ptr->vx = fluidCell_ptr->vx;
324  fluid_cell_info_ptr->vy = fluidCell_ptr->vy;
325  fluid_cell_info_ptr->vz = fluidCell_ptr->vz;
326  fluid_cell_info_ptr->mu_B = 0.0;
327  fluid_cell_info_ptr->mu_C = 0.0;
328  fluid_cell_info_ptr->mu_S = 0.0;
329  fluid_cell_info_ptr->qgp_fraction = 0.0;
330 
331  for (int i = 0; i < 4; i++) {
332  for (int j = 0; j < 4; j++) {
333  fluid_cell_info_ptr->pi[i][j] = fluidCell_ptr->pi[i][j];
334  }
335  }
336  fluid_cell_info_ptr->bulk_Pi = fluidCell_ptr->bulkPi;
337  delete fluidCell_ptr;
338 }