vhlle/src/ic.cpp

0001 #include <fstream>
0002 #include <iomanip>
0003 #include <iomanip>
0004 #include <TF1.h>
0005 #include <TF2.h>
0006 #include <TGraph.h>
0007
0008 #include "fld.h"
0009 #include "eos.h"
0010 #include "ic.h"
0011 #include "inc.h"
0012 #include "s95p.h"
0013
0014 using namespace std;
0015
0016
0017 IC::IC(char* icInputFile, double s0ScaleFactor)
0018 {
0019  s95p::loadSongIC(icInputFile, s0ScaleFactor);
0020 }
0021
0022 IC::~IC(void) {}
0023
0024
0025 void IC::setIC(Fluid *f, EoS *eos, double tau) {
0026   double e, nb, nq, vx = 0., vy = 0., vz = 0.;
0027   Cell *c;
0028
0029   double avv_num = 0., avv_den = 0.;
0030   double Etotal = 0.0;
0031
0032   for (int ix = 0; ix < f->getNX(); ix++)
0033     for (int iy = 0; iy < f->getNY(); iy++)
0034       for (int iz = 0; iz < f->getNZ(); iz++) {
0035         c = f->getCell(ix, iy, iz);
0036         double x = f->getX(ix);
0037         double y = f->getY(iy);
0038         double eta = f->getZ(iz);
0039
0040         double eta1 = fabs(eta) < 1.3 ? 0.0 : fabs(eta) - 1.3;
0041         double etaFactor =
0042               exp(-eta1 * eta1 / 2.1 / 2.1) * (fabs(eta) < 5.3 ? 1.0 : 0.0);
0043         e = s95p::getSongEps(x, y) * etaFactor;
0044         if (e < 0.5) e = 0.0;
0045         vx = vy = 0.0;
0046         nb = nq = 0.0;
0047         vz = 0.0;
0048
0049         avv_num += sqrt(vx * vx + vy * vy) * e;
0050         avv_den += e;
0051
0052         c->setPrimVar(eos, tau, e, nb, nq, 0., vx, vy, vz);
0053         double _p = eos->p(e, nb, nq, 0.);
0054         const double gamma2 = 1.0 / (1.0 - vx * vx - vy * vy - vz * vz);
0055         Etotal +=
0056             ((e + _p) * gamma2 * (cosh(eta) + vz * sinh(eta)) - _p * cosh(eta));
0057         c->saveQprev();
0058
0059         if (e > 0.) c->setAllM(1.);
0060       }
0061   cout << "average initial flow = " << avv_num / avv_den << endl;
0062   cout << "total energy = " << Etotal *f->getDx() * f->getDy() * f->getDz() *
0063                                    tau << endl;
0064 }