Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
main.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file main.cpp
1 /******************************************************************************
2 * *
3 * vHLLE : a 3D viscous hydrodynamic code *
4 * version 1.0, November 2013 *
5 * by Iurii Karpenko *
6 * contact: yu.karpenko@gmail.com *
7 * For the detailed description please refer to: *
8 * http://arxiv.org/abs/1312.4160 *
9 * *
10 * This code can be freely used and redistributed, provided that this *
11 * copyright appear in all the copies. If you decide to make modifications *
12 * to the code, please contact the authors, especially if you plan to publish *
13 * the results obtained with such modified code. Any publication of results *
14 * obtained using this code must include the reference to *
15 * arXiv:1312.4160 [nucl-th] or the published version of it, when available. *
16 * *
17 *******************************************************************************/
18 
19 #include <iostream>
20 #include <fstream>
21 #include <iomanip>
22 #include <cstring>
23 #include <ctime>
24 #include <sstream>
25 #include "fld.h"
26 #include "hdo.h"
27 #include "ic.h"
28 #include "icGlauber.h"
29 #include "icGubser.h"
30 #include "eos.h"
31 #include "trancoeff.h"
32 
33 using namespace std;
34 
35 // program parameters, to be read from file
36 int nx, ny, nz, eosType;
37 double xmin, xmax, ymin, ymax, etamin, etamax, tau0, tauMax, dtau;
38 char outputDir[255];
39 char icInputFile[255];
40 double etaS, zetaS, eCrit;
41 int icModel,
42  glauberVariable =
43  1; // icModel=1 for pure Glauber, 2 for table input (Glissando etc)
44 double epsilon0, alpha, impactPar, s0ScaleFactor;
45 
46 void readParameters(char *parFile) {
47  char parName[255], parValue[255];
48  ifstream fin(parFile);
49  if (!fin.is_open()) {
50  cout << "cannot open parameters file " << parFile << endl;
51  exit(1);
52  }
53  while (fin.good()) {
54  string line;
55  getline(fin, line);
56  istringstream sline(line);
57  sline >> parName >> parValue;
58  if (strcmp(parName, "outputDir") == 0)
59  strcpy(outputDir, parValue);
60  else if (strcmp(parName, "eosType") == 0)
61  eosType = atoi(parValue);
62  else if (strcmp(parName, "icInputFile") == 0)
63  strcpy(icInputFile, parValue);
64  else if (strcmp(parName, "nx") == 0)
65  nx = atoi(parValue);
66  else if (strcmp(parName, "ny") == 0)
67  ny = atoi(parValue);
68  else if (strcmp(parName, "nz") == 0)
69  nz = atoi(parValue);
70  else if (strcmp(parName, "icModel") == 0)
71  icModel = atoi(parValue);
72  else if (strcmp(parName, "glauberVar") == 0)
73  glauberVariable = atoi(parValue);
74  else if (strcmp(parName, "xmin") == 0)
75  xmin = atof(parValue);
76  else if (strcmp(parName, "xmax") == 0)
77  xmax = atof(parValue);
78  else if (strcmp(parName, "ymin") == 0)
79  ymin = atof(parValue);
80  else if (strcmp(parName, "ymax") == 0)
81  ymax = atof(parValue);
82  else if (strcmp(parName, "etamin") == 0)
83  etamin = atof(parValue);
84  else if (strcmp(parName, "etamax") == 0)
85  etamax = atof(parValue);
86  else if (strcmp(parName, "tau0") == 0)
87  tau0 = atof(parValue);
88  else if (strcmp(parName, "tauMax") == 0)
89  tauMax = atof(parValue);
90  else if (strcmp(parName, "dtau") == 0)
91  dtau = atof(parValue);
92  else if (strcmp(parName, "e_crit") == 0)
93  eCrit = atof(parValue);
94  else if (strcmp(parName, "etaS") == 0)
95  etaS = atof(parValue);
96  else if (strcmp(parName, "zetaS") == 0)
97  zetaS = atof(parValue);
98  else if (strcmp(parName, "epsilon0") == 0)
99  epsilon0 = atof(parValue);
100  else if (strcmp(parName, "alpha") == 0)
101  alpha = atof(parValue);
102  else if (strcmp(parName, "impactPar") == 0)
103  impactPar = atof(parValue);
104  else if (strcmp(parName, "s0ScaleFactor") == 0)
105  s0ScaleFactor = atof(parValue);
106  else if (parName[0] == '!')
107  cout << "CCC " << sline.str() << endl;
108  else
109  cout << "UUU " << sline.str() << endl;
110  }
111 }
112 
113 void printParameters() {
114  cout << "====== parameters ======\n";
115  cout << "outputDir = " << outputDir << endl;
116  cout << "eosType = " << eosType << endl;
117  cout << "nx = " << nx << endl;
118  cout << "ny = " << ny << endl;
119  cout << "nz = " << nz << endl;
120  cout << "icModel = " << icModel << endl;
121  cout << "glauberVar = " << glauberVariable << " ! 0=epsilon,1=entropy"
122  << endl;
123  cout << "xmin = " << xmin << endl;
124  cout << "xmax = " << xmax << endl;
125  cout << "ymin = " << ymin << endl;
126  cout << "ymax = " << ymax << endl;
127  cout << "etamin = " << etamin << endl;
128  cout << "etamax = " << etamax << endl;
129  cout << "tau0 = " << tau0 << endl;
130  cout << "tauMax = " << tauMax << endl;
131  cout << "dtau = " << dtau << endl;
132  cout << "e_crit = " << eCrit << endl;
133  cout << "eta/s = " << etaS << endl;
134  cout << "zeta/s = " << zetaS << endl;
135  cout << "epsilon0 = " << epsilon0 << endl;
136  cout << "alpha = " << alpha << endl;
137  cout << "impactPar = " << impactPar << endl;
138  cout << "s0ScaleFactor = " << s0ScaleFactor << endl;
139  cout << "======= end parameters =======\n";
140 }
141 
142 // program parameters, to be read from file
143 // int nx, ny, nz, eosType ;
144 // double xmin, xmax, ymin, ymax, zmin, zmax, tau0, tauMax, dtau ;
145 // char outputDir[255], eosFile[255], chiBfile[255], chiSfile[255] ;
146 // char icInputFile [255] ;
147 // double T_ch, mu_b, mu_q, mu_s, gammaS, gammaFactor, exclVolume ;
148 // int icModel, NPART, glauberVariable=1 ;
149 // double epsilon0, alpha, impactPar, s0ScaleFactor ;
150 
151 int main(int argc, char **argv) {
152  // pointers to all the main objects
153  EoS *eos;
154  TransportCoeff *trcoeff;
155  Fluid *f;
156  Hydro *h;
157  time_t start = 0, end;
158 
159  time(&start);
160 
161  // read parameters from file
162  char *parFile;
163  if (argc == 1) {
164  cout << "NO PARAMETERS, exiting\n";
165  exit(1);
166  } else {
167  parFile = argv[1];
168  }
169  readParameters(parFile);
170  printParameters();
171 
172  // EoS
173  char eosfile[] = "eos/Laine_nf3.dat";
174  int ncols = 3;
175  eos = new EoSs(eosfile, ncols);
176 
177  // transport coefficients
178  trcoeff = new TransportCoeff(etaS, zetaS, eos);
179 
180  f = new Fluid(eos, eos, trcoeff, nx, ny, nz, xmin, xmax, ymin, ymax, etamin,
181  etamax, dtau, eCrit);
182  cout << "fluid allocation done\n";
183 
184  // initilal conditions
185  if(icModel==1){ // optical Glauber
186  ICGlauber *ic = new ICGlauber(epsilon0, impactPar, tau0);
187  ic->setIC(f, eos);
188  delete ic;
189  }else if(icModel==2){ // Glauber_table + parametrized rapidity dependence
190  IC *ic = new IC(icInputFile, s0ScaleFactor);
191  ic->setIC(f, eos, tau0);
192  delete ic;
193  }else if(icModel==4){ // analytical Gubser solution
194  ICGubser *ic = new ICGubser();
195  ic->setIC(f, eos, tau0);
196  delete ic;
197  }else{
198  cout << "icModel = " << icModel << " not implemented\n";
199  }
200  cout << "IC done\n";
201 
202  time_t tinit = 0;
203  time(&tinit);
204  float diff = difftime(tinit, start);
205  cout << "Init time = " << diff << " [sec]" << endl;
206 
207  // hydro init
208  h = new Hydro(f, eos, trcoeff, tau0, dtau);
209  int maxstep = ceil((tauMax - tau0) / dtau);
210  start = 0;
211  time(&start);
212  h->setNSvalues(); // initialize viscous terms
213 
214  f->initOutput(outputDir, maxstep, tau0, 2);
215  f->calcTotals(h->getTau());
216 
217  for (int istep = 0; istep < maxstep; istep++) {
218  // decrease timestep automatically, but use fixed dtau for output
219  int nSubSteps = 1;
220  while (dtau / nSubSteps > 0.5 * (tau0 + dtau * istep))
221  nSubSteps *= 2; // 0.02 in "old" coordinates
222  h->setDtau(dtau / nSubSteps);
223  for (int j = 0; j < nSubSteps; j++) {
224  h->performStep();
225  }
226  cout << "step= " << istep << " dtau= " << dtau / nSubSteps << "\n"
227  << endl; // "\r" << flush
228  f->outputGnuplot(h->getTau());
229  f->calcTotals(h->getTau());
230  }
231 
232  end = 0;
233  time(&end);
234  float diff2 = difftime(end, start);
235  cout << "Execution time = " << diff2 << " [sec]" << endl;
236 
237  delete f;
238  delete h;
239  delete eos;
240 }