external_packages/hydro_from_external_file/main.cpp

0001 /////////////////////////////////////////////////////////////////////////
0002 //                      hydrodynamics analysis
0003 //
0004 //              author: Chun Shen <chunshen@physics.mcgill.ca>
0005 //              Copyright: Chun Shen 2014
0006 //
0007 //  This program load hydrodynamic evolution files and perform various
0008 //  kinds of analysis
0009 //
0010 //
0011 //  To do in the future:
0012 /////////////////////////////////////////////////////////////////////////
0013
0014 #include <iostream>
0015 #include <sstream>
0016 #include <fstream>
0017 #include <cmath>
0018 #include <iomanip>
0019 #include <cstdlib>
0020
0021 #ifdef USE_HDF5
0022 #include "./Hydroinfo_h5.h"
0023 #endif
0024
0025 #include "./Hydroinfo_MUSIC.h"
0026 #include "./Stopwatch.h"
0027 #include "./FluidcellStatistic.h"
0028 #include "./ParameterReader.h"
0029 #include "./SurfaceFinder.h"
0030
0031 using namespace std;
0032
0033 int main(int argc, char *argv[]) {
0034     ParameterReader *paraRdr = new ParameterReader();
0035     paraRdr->readFromFile("parameters.dat");
0036     paraRdr->readFromArguments(argc, argv);
0037     paraRdr->echo();
0038
0039     int load_viscous = paraRdr->getVal("load_viscous_info");
0040     int hydro_type = paraRdr->getVal("hydro_type");
0041
0042     void* hydroinfo_ptr_in = NULL;
0043
0044     string input_file = "results/music_input";
0045     string hydro_ideal_file = "results/evolution_xyeta.dat";
0046     string hydro_shear_file =
0047                     "results/evolution_Wmunu_over_epsilon_plus_P_xyeta.dat";
0048     string hydro_bulk_file = "results/evolution_bulk_pressure_xyeta.dat";
0049
0050     Stopwatch sw;
0051     sw.tic();
0052     // hydro data file pointer
0053     if (hydro_type == 1) {
0054 #ifdef USE_HDF5
0055         HydroinfoH5* hydroinfo_ptr = new HydroinfoH5("JetData.h5", 500,
0056                                                      load_viscous);
0057         hydroinfo_ptr_in = hydroinfo_ptr;
0058 #endif
0059     } else if (hydro_type == 2) {
0060         Hydroinfo_MUSIC* hydroinfo_ptr = new Hydroinfo_MUSIC();
0061         int hydro_mode = 8;
0062         int nskip_tau = paraRdr->getVal("hydro_nskip_tau");
0063         hydroinfo_ptr->readHydroData(hydro_mode, nskip_tau,
0064             input_file, hydro_ideal_file, hydro_shear_file, hydro_bulk_file);
0065         hydroinfo_ptr_in = hydroinfo_ptr;
0066     } else if (hydro_type == 3) {
0067         Hydroinfo_MUSIC* hydroinfo_ptr = new Hydroinfo_MUSIC();
0068         int hydro_mode = 9;
0069         int nskip_tau = paraRdr->getVal("hydro_nskip_tau");
0070         hydroinfo_ptr->readHydroData(hydro_mode, nskip_tau,
0071             input_file, hydro_ideal_file, hydro_shear_file, hydro_bulk_file);
0072         hydroinfo_ptr_in = hydroinfo_ptr;
0073     } else if (hydro_type == 4) {
0074         Hydroinfo_MUSIC* hydroinfo_ptr = new Hydroinfo_MUSIC();
0075         int hydro_mode = 10;
0076         int nskip_tau = 1;
0077         hydroinfo_ptr->readHydroData(hydro_mode, nskip_tau,
0078             input_file, hydro_ideal_file, hydro_shear_file, hydro_bulk_file);
0079         hydroinfo_ptr_in = hydroinfo_ptr;
0080     } else {
0081         cout << "main: unrecognized hydro_type = " << hydro_type << endl;
0082         exit(1);
0083     }
0084
0085     FluidcellStatistic fluidcellanalysis(hydroinfo_ptr_in, paraRdr);
0086     double T_cut = paraRdr->getVal("T_cut");
0087     fluidcellanalysis.outputTempasTauvsX();
0088     fluidcellanalysis.outputKnudersonNumberasTauvsX();
0089     fluidcellanalysis.outputinverseReynoldsNumberasTauvsX();
0090     fluidcellanalysis.analysis_hydro_volume_for_photon(T_cut);
0091     fluidcellanalysis.output_temperature_vs_avg_utau();
0092     fluidcellanalysis.output_flowvelocity_vs_tau();
0093
0094     // construct freeze-out hyper-surface
0095     // SurfaceFinder* surface_ptr = new SurfaceFinder(hydroinfo_ptr, paraRdr);
0096     // surface_ptr->Find_full_hypersurface();
0097
0098     sw.toc();
0099     cout << "totally takes : " << sw.takeTime() << " seconds." << endl;
0100
0101     return(0);
0102 }
0103
0104