vhlle/src/ctrl.cpp

0001 #include "ctrl.h"
0002
0003 #include "inc.h"
0004 #include <ctime>
0005 #include <cmath>
0006 #include <iomanip>
0007 //#include <TMath.h>
0008 #include "hdo.h"
0009 #include "ickw.h"
0010 #include "conv.h"
0011 #include "eo3.h"
0012 #include "eo1.h"
0013 #include "trancoeff.h"
0014 #include <unistd.h>
0015
0016 #include <execinfo.h>
0017 #include <signal.h>
0018
0019 ofstream ofile ; // cout + cerr output
0020
0021 EoS *eos ;
0022 TransportCoeff *trcoeff ;
0023 Fluid *f ;
0024 IC_KW *ic_kw ;
0025 Hydro* h ;
0026 time_t start, end ;
0027
0028 double tau0, tau_max, dtau ;
0029 int maxstep ;
0030 //string directory, outkw ;
0031
0032 void handler(int sig) {
0033   void *array[10];
0034   size_t size;
0035
0036   // get void*'s for all entries on the stack
0037   size = backtrace(array, 10);
0038
0039   // print out all the frames to stderr
0040   cout<<"Error: signal "<<sig<<endl ;
0041   backtrace_symbols_fd(array, size, 2);
0042   exit(1);
0043 }
0044
0045 void redirecterrorshlle_(char* filename)
0046 {
0047     signal(SIGSEGV, handler); // handler for segfault
0048     ofile.open(filename) ;
0049     cout.rdbuf(ofile.rdbuf()) ;
0050     cerr.rdbuf(ofile.rdbuf()) ;
0051 }
0052
0053 void initeoshlle_(char *filename, int* ncols)
0054 {
0055     eos = new EoSs(filename,*ncols) ; // << CE ; set BAG_ASYMPT !!
0056 }
0057
0058
0059 void initeoshlle3f_(char *filename, double *B, double *volex0, double *delta0, double *aaa, double *bbb)
0060 {
0061     eos = new EoS3f(filename, *B, *volex0, *delta0, *aaa, *bbb) ;
0062 }
0063
0064 void initeoshlle1f_(char *filename)
0065 {
0066     eos = new EoS1f(filename) ;
0067 }
0068
0069
0070 void inittrcoeff_(double *etaS, double *zetaS)
0071 {
0072  trcoeff = new TransportCoeff(*etaS, *zetaS, eos) ;
0073 }
0074
0075 void eoshlle_(double *e, double *nb, double *nq, double *ns, double *T, double *mub, double *muq, double *mus, double *p)
0076 {
0077     eos->eos(*e, *nb, *nq, *ns, *T, *mub, *muq, *mus, *p) ;
0078 }
0079
0080
0081 void eosrangeshlle_(double *emax, double *e0, double *nmax, double *n0, int *ne, int *nn)
0082 {
0083     eos->eosranges(*emax, *e0, *nmax, *n0, *ne, *nn) ;
0084 }
0085
0086
0087 void initfluidhlle_(int* nx, int* ny, int* nz, double* minx, double* maxx, double* miny, double* maxy, double* minz, double* maxz)
0088 {
0089     f = new Fluid(eos, trcoeff, *nx, *ny, *nz, *minx, *maxx, *miny, *maxy, *minz, *maxz) ;
0090 }
0091
0092
0093
0094 void initichlle_(char *filename, double *tau0)
0095 {
0096     ic_kw = new IC_KW(filename) ; //---------for Klaus
0097 //    ic_kw->setIC(f, eos, *tau0) ; //!!! z-symmetry is included in ic_kw.cpp
0098 }
0099
0100
0101 void icgethlle3f_(double* x, double* y, double* eta, double* e, double* nB, double* nQ, double* nS, double* vx, double* vy, double* vz)
0102 {
0103     double nu, nd, ns ;
0104     ic_kw->getICs(*x, *y, *eta, *e, *vx, *vy, *vz, nu, nd, ns) ;
0105 //  if(*e<1.) *e = 0. ;
0106     *nB = 1./3.*(nu + nd + ns) ;
0107     *nQ = 1./3.*(2.*nu - nd - ns) ;
0108     *nS = - ns ;
0109 }
0110
0111
0112 void icgethlle_(double* x, double* y, double* eta, double* e, double* nB, double* nQ, double* nS, double* vx, double* vy, double* vz)
0113 {
0114     double nu, nd, ns ;
0115     ic_kw->getICs(*x, *y, *eta, *e, *vx, *vy, *vz, nu, nd, ns) ;
0116     *nB = *nQ = *nS = 0. ;
0117 }
0118
0119
0120 void icsethlle_(int* ix, int* iy, int* iz, double* tau0, double* e, double* nb, double* nq, double* ns, double* vx, double* vy, double* vz)
0121 {
0122 Cell *c = f->getCell(*ix-1,*iy-1,*iz-1) ;
0123 if((*vx)*(*vx)+(*vy)*(*vy)+(*vz)*(*vz)>1.){
0124 //  cerr << "setIC : " << ix << "  " << iy << "  " << iz << "  e = " << e << "  v^2 = " << (*vx)*(*vx)+(*vy)*(*vy)+vz*vz << endl ;
0125     double factor = sqrt((*vx)*(*vx)+(*vy)*(*vy)+(*vz)*(*vz)) ;
0126     (*vx) = (*vx)*0.99/factor ;
0127     (*vy) = (*vy)*0.99/factor ;
0128     (*vz) = (*vz)*0.99/factor ;
0129   }
0130   if(isinf(*e) || isnan(*e) || isinf(*nb) || isnan(*nb) || isinf(*nq) || isnan(*nq) || isinf(*ns) || isnan(*ns) ||
0131   isinf(*vx) || isnan(*vx) || isinf(*vy) || isnan(*vy) || isinf(*vz) || isnan(*vz))
0132   {cout<<"fluid: bad init,"<<setw(14)<<"e"<<setw(14)<<"nb"<<setw(14)<<"nq"<<setw(14)<<"ns"<<setw(14)<<"vx"<<setw(14)<<"vy"<<setw(14)<<"vz"<<endl ;
0133    cout<<"----------------"<<setw(14)<<e<<setw(14)<<nb<<setw(14)<<nq<<setw(14)<<ns<<setw(14)<<vx<<setw(14)<<vy<<setw(14)<<vz<<endl ;
0134    cout<<"cell  "<<ix<<"  "<<iy<<"  "<<iz<<endl ;
0135    exit(1) ; }
0136   c->setPrimVar(eos, *tau0, *e, *nb, *nq, *ns, (*vx), (*vy), (*vz)) ;
0137   c->saveQprev() ;
0138   if(*e>0.) c->setAllM(1.) ;
0139 }
0140
0141
0142 void inithydrohlle_(double* _tau0, double* _tau_max, double* _dtau)
0143 {
0144     tau0 = *_tau0 ;
0145     tau_max = *_tau_max ;
0146     dtau = *_dtau ;
0147     h = new Hydro(f, eos, trcoeff, tau0, dtau) ;
0148     maxstep = ceil((tau_max-tau0)/dtau) ;
0149
0150     start = 0;
0151     time(&start);
0152     h->setNSvalues() ; // initialize viscous terms
0153 }
0154
0155
0156 void dtauhlle_(double* dtau)
0157 {
0158   h->setDtau(*dtau) ;
0159 }
0160
0161
0162 void initoutputhlle_(char* dir)
0163 {
0164   f->initOutput(dir, maxstep, tau0, 2) ;
0165   f->outputPDirections(h->getTau());
0166   f->outputTransverseAverages(h->getTau()) ;
0167 }
0168
0169
0170 int getmaxstephlle_(void)
0171 {
0172     return maxstep ;
0173 }
0174
0175
0176 void makestephlle_(int *i)
0177 {
0178   h->performStep() ;
0179     f->outputPDirections(h->getTau());
0180   f->outputTransverseAverages(h->getTau()) ;
0181 }
0182
0183
0184 void getvalueshlle_(int* ix, int* iy, int* iz, double* e, double *p, double *nb, double *nq, double *ns, double* vx, double* vy, double* vz, double* viscCorrCutFlag)
0185 {
0186     Cell *c = f->getCell(*ix-1, *iy-1, *iz-1) ;
0187     c->getPrimVar(eos, h->getTau(), *e, *p, *nb, *nq, *ns, *vx, *vy, *vz) ;
0188     *viscCorrCutFlag = c->getViscCorrCutFlag() ;
0189 }
0190
0191
0192 void getvflaghlle_(int* ix, int* iy, int* iz, double* viscCorrCutFlag)
0193 {
0194   *viscCorrCutFlag = f->getCell(*ix-1, *iy-1, *iz-1)->getViscCorrCutFlag() ;
0195 }
0196
0197
0198 void getvischlle_(int* ix, int* iy, int* iz, double *pi, double *Pi)
0199 {
0200  Cell* c=f->getCell(*ix-1, *iy-1, *iz-1) ;
0201  // fortran and C have reverse array alignment, (a,b) --> [b,a]
0202  // but since pi[mu][nu] is symmetric, this is not important
0203  for(int i=0; i<4; i++)
0204  for(int j=0; j<4; j++){
0205   *(pi+4*i+j) = c->getpi(i,j) ;
0206  }
0207  *Pi=c->getPi();
0208  //if(c->Pi!=0.0) cout <<"nonzero Pi: " << c->Pi << endl ;
0209 }
0210
0211
0212 double getxhlle_(int *ix)
0213 {
0214     return f->getX(*ix-1) ;
0215 }
0216
0217
0218 double getyhlle_(int *iy)
0219 {
0220     return f->getY(*iy-1) ;
0221 }
0222
0223
0224 double getzhlle_(int *iz)
0225 {
0226     return f->getZ(*iz-1) ;
0227 }
0228
0229
0230 void destroyhlle_(void)
0231 {
0232 //    delete ic ; // do we need IC instance?
0233     delete f ;
0234     delete h ;
0235     end=0 ;
0236     time(&end); float diff = difftime(end, start);
0237 //    cout<<"Execution time = "<<diff<< " [sec]" << endl;
0238     ofile.close() ; // close output and error file
0239 }
0240
0241
0242 void destroyeoshlle_(void)
0243 {
0244     delete eos ;
0245 }
0246
0247
0248 double gettimehlle_(void)
0249 {
0250     end=0 ;
0251     time(&end); float diff = difftime(end, start);
0252     return diff ;
0253 }
0254
0255 double getenergyhlle_(void)
0256 {
0257     double ene = 0., e, p, nb, nq, ns, vx, vy, vz, _vx, _vy, _vz ;
0258     double tau = h->getTau() ;
0259     for(int ix=0; ix<f->getNX(); ix++)
0260     for(int iy=0; iy<f->getNY(); iy++)
0261     for(int iz=0; iz<f->getNZ(); iz++){
0262     f->getCell(ix, iy, iz)->getPrimVar(eos, tau, e, p, nb, nq, ns, _vx, _vy, _vz) ;
0263     double eta = f->getZ(iz) ;
0264     vx = _vx/(cosh(eta) + _vz*sinh(eta)) ;
0265     vy = _vy/(cosh(eta) + _vz*sinh(eta)) ;
0266     vz = (_vz*cosh(eta)+sinh(eta))/(cosh(eta) + _vz*sinh(eta)) ;
0267     ene += (e+p)/(1.-vx*vx-vy*vy-vz*vz) - p ;
0268     }
0269     return e ;
0270 }