vhlle/src/main.cpp

0001 /******************************************************************************
0002 *                                                                             *
0003 *            vHLLE : a 3D viscous hydrodynamic code                           *
0004 *            version 1.0,            November 2013                            *
0005 *            by Iurii Karpenko                                                *
0006 *  contact:  yu.karpenko@gmail.com                                            *
0007 *  For the detailed description please refer to:                              *
0008 *  http://arxiv.org/abs/1312.4160                                             *
0009 *                                                                             *
0010 *  This code can be freely used and redistributed, provided that this         *
0011 *  copyright appear in all the copies. If you decide to make modifications    *
0012 *  to the code, please contact the authors, especially if you plan to publish *
0013 * the results obtained with such modified code. Any publication of results    *
0014 * obtained using this code must include the reference to                      *
0015 * arXiv:1312.4160 [nucl-th] or the published version of it, when available.   *
0016 *                                                                             *
0017 *******************************************************************************/
0018
0019 #include <iostream>
0020 #include <fstream>
0021 #include <iomanip>
0022 #include <cstring>
0023 #include <ctime>
0024 #include <sstream>
0025 #include "fld.h"
0026 #include "hdo.h"
0027 #include "ic.h"
0028 #include "icGlauber.h"
0029 #include "icGubser.h"
0030 #include "eos.h"
0031 #include "trancoeff.h"
0032
0033 using namespace std;
0034
0035 // program parameters, to be read from file
0036 int nx, ny, nz, eosType;
0037 double xmin, xmax, ymin, ymax, etamin, etamax, tau0, tauMax, dtau;
0038 char outputDir[255];
0039 char icInputFile[255];
0040 double etaS, zetaS, eCrit;
0041 int icModel,
0042     glauberVariable =
0043         1;  // icModel=1 for pure Glauber, 2 for table input (Glissando etc)
0044 double epsilon0, alpha, impactPar, s0ScaleFactor;
0045
0046 void readParameters(char *parFile) {
0047   char parName[255], parValue[255];
0048   ifstream fin(parFile);
0049   if (!fin.is_open()) {
0050     cout << "cannot open parameters file " << parFile << endl;
0051     exit(1);
0052   }
0053   while (fin.good()) {
0054     string line;
0055     getline(fin, line);
0056     istringstream sline(line);
0057     sline >> parName >> parValue;
0058     if (strcmp(parName, "outputDir") == 0)
0059       strcpy(outputDir, parValue);
0060     else if (strcmp(parName, "eosType") == 0)
0061       eosType = atoi(parValue);
0062     else if (strcmp(parName, "icInputFile") == 0)
0063       strcpy(icInputFile, parValue);
0064     else if (strcmp(parName, "nx") == 0)
0065       nx = atoi(parValue);
0066     else if (strcmp(parName, "ny") == 0)
0067       ny = atoi(parValue);
0068     else if (strcmp(parName, "nz") == 0)
0069       nz = atoi(parValue);
0070     else if (strcmp(parName, "icModel") == 0)
0071       icModel = atoi(parValue);
0072     else if (strcmp(parName, "glauberVar") == 0)
0073       glauberVariable = atoi(parValue);
0074     else if (strcmp(parName, "xmin") == 0)
0075       xmin = atof(parValue);
0076     else if (strcmp(parName, "xmax") == 0)
0077       xmax = atof(parValue);
0078     else if (strcmp(parName, "ymin") == 0)
0079       ymin = atof(parValue);
0080     else if (strcmp(parName, "ymax") == 0)
0081       ymax = atof(parValue);
0082     else if (strcmp(parName, "etamin") == 0)
0083       etamin = atof(parValue);
0084     else if (strcmp(parName, "etamax") == 0)
0085       etamax = atof(parValue);
0086     else if (strcmp(parName, "tau0") == 0)
0087       tau0 = atof(parValue);
0088     else if (strcmp(parName, "tauMax") == 0)
0089       tauMax = atof(parValue);
0090     else if (strcmp(parName, "dtau") == 0)
0091       dtau = atof(parValue);
0092     else if (strcmp(parName, "e_crit") == 0)
0093       eCrit = atof(parValue);
0094     else if (strcmp(parName, "etaS") == 0)
0095       etaS = atof(parValue);
0096     else if (strcmp(parName, "zetaS") == 0)
0097       zetaS = atof(parValue);
0098     else if (strcmp(parName, "epsilon0") == 0)
0099       epsilon0 = atof(parValue);
0100     else if (strcmp(parName, "alpha") == 0)
0101       alpha = atof(parValue);
0102     else if (strcmp(parName, "impactPar") == 0)
0103       impactPar = atof(parValue);
0104     else if (strcmp(parName, "s0ScaleFactor") == 0)
0105       s0ScaleFactor = atof(parValue);
0106     else if (parName[0] == '!')
0107       cout << "CCC " << sline.str() << endl;
0108     else
0109       cout << "UUU " << sline.str() << endl;
0110   }
0111 }
0112
0113 void printParameters() {
0114   cout << "====== parameters ======\n";
0115   cout << "outputDir = " << outputDir << endl;
0116   cout << "eosType = " << eosType << endl;
0117   cout << "nx = " << nx << endl;
0118   cout << "ny = " << ny << endl;
0119   cout << "nz = " << nz << endl;
0120   cout << "icModel = " << icModel << endl;
0121   cout << "glauberVar = " << glauberVariable << "   ! 0=epsilon,1=entropy"
0122        << endl;
0123   cout << "xmin = " << xmin << endl;
0124   cout << "xmax = " << xmax << endl;
0125   cout << "ymin = " << ymin << endl;
0126   cout << "ymax = " << ymax << endl;
0127   cout << "etamin = " << etamin << endl;
0128   cout << "etamax = " << etamax << endl;
0129   cout << "tau0 = " << tau0 << endl;
0130   cout << "tauMax = " << tauMax << endl;
0131   cout << "dtau = " << dtau << endl;
0132   cout << "e_crit = " << eCrit << endl;
0133   cout << "eta/s = " << etaS << endl;
0134   cout << "zeta/s = " << zetaS << endl;
0135   cout << "epsilon0 = " << epsilon0 << endl;
0136   cout << "alpha = " << alpha << endl;
0137   cout << "impactPar = " << impactPar << endl;
0138   cout << "s0ScaleFactor = " << s0ScaleFactor << endl;
0139   cout << "======= end parameters =======\n";
0140 }
0141
0142 // program parameters, to be read from file
0143 // int nx, ny, nz, eosType ;
0144 // double xmin, xmax, ymin, ymax, zmin, zmax, tau0, tauMax, dtau ;
0145 // char outputDir[255], eosFile[255], chiBfile[255], chiSfile[255] ;
0146 // char icInputFile [255] ;
0147 // double T_ch, mu_b, mu_q, mu_s, gammaS, gammaFactor, exclVolume ;
0148 // int icModel, NPART, glauberVariable=1 ;
0149 // double epsilon0, alpha, impactPar, s0ScaleFactor ;
0150
0151 int main(int argc, char **argv) {
0152   // pointers to all the main objects
0153   EoS *eos;
0154   TransportCoeff *trcoeff;
0155   Fluid *f;
0156   Hydro *h;
0157   time_t start = 0, end;
0158
0159   time(&start);
0160
0161   // read parameters from file
0162   char *parFile;
0163   if (argc == 1) {
0164     cout << "NO PARAMETERS, exiting\n";
0165     exit(1);
0166   } else {
0167     parFile = argv[1];
0168   }
0169   readParameters(parFile);
0170   printParameters();
0171
0172   // EoS
0173   char eosfile[] = "eos/Laine_nf3.dat";
0174   int ncols = 3;
0175   eos = new EoSs(eosfile, ncols);
0176
0177   // transport coefficients
0178   trcoeff = new TransportCoeff(etaS, zetaS, eos);
0179
0180   f = new Fluid(eos, eos, trcoeff, nx, ny, nz, xmin, xmax, ymin, ymax, etamin,
0181                 etamax, dtau, eCrit);
0182   cout << "fluid allocation done\n";
0183
0184   // initilal conditions
0185   if(icModel==1){ // optical Glauber
0186    ICGlauber *ic = new ICGlauber(epsilon0, impactPar, tau0);
0187    ic->setIC(f, eos);
0188    delete ic;
0189   }else if(icModel==2){ // Glauber_table + parametrized rapidity dependence
0190    IC *ic = new IC(icInputFile, s0ScaleFactor);
0191    ic->setIC(f, eos, tau0);
0192    delete ic;
0193   }else if(icModel==4){ // analytical Gubser solution
0194    ICGubser *ic = new ICGubser();
0195    ic->setIC(f, eos, tau0);
0196    delete ic;
0197   }else{
0198    cout << "icModel = " << icModel << " not implemented\n";
0199   }
0200   cout << "IC done\n";
0201
0202   time_t tinit = 0;
0203   time(&tinit);
0204   float diff = difftime(tinit, start);
0205   cout << "Init time = " << diff << " [sec]" << endl;
0206
0207   // hydro init
0208   h = new Hydro(f, eos, trcoeff, tau0, dtau);
0209   int maxstep = ceil((tauMax - tau0) / dtau);
0210   start = 0;
0211   time(&start);
0212   h->setNSvalues();  // initialize viscous terms
0213
0214   f->initOutput(outputDir, maxstep, tau0, 2);
0215   f->calcTotals(h->getTau());
0216
0217   for (int istep = 0; istep < maxstep; istep++) {
0218     // decrease timestep automatically, but use fixed dtau for output
0219     int nSubSteps = 1;
0220     while (dtau / nSubSteps > 0.5 * (tau0 + dtau * istep))
0221       nSubSteps *= 2;  // 0.02 in "old" coordinates
0222     h->setDtau(dtau / nSubSteps);
0223     for (int j = 0; j < nSubSteps; j++) {
0224       h->performStep();
0225     }
0226     cout << "step= " << istep << "  dtau= " << dtau / nSubSteps << "\n"
0227          << endl;  // "\r" << flush
0228     f->outputGnuplot(h->getTau());
0229     f->calcTotals(h->getTau());
0230   }
0231
0232   end = 0;
0233   time(&end);
0234   float diff2 = difftime(end, start);
0235   cout << "Execution time = " << diff2 << " [sec]" << endl;
0236
0237   delete f;
0238   delete h;
0239   delete eos;
0240 }