File indexing completed on 2025-08-05 08:19:16
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 "icGubser.h"
0011 #include "inc.h"
0012 #include "s95p.h"
0013
0014 using namespace std;
0015
0016 ICGubser::ICGubser() {
0017 }
0018
0019 ICGubser::~ICGubser(void) {}
0020
0021
0022 void ICGubser::setIC(Fluid *f, EoS *eos, double tau) {
0023 double e, nb, nq, vx = 0., vy = 0., vz = 0.;
0024 Cell *c;
0025
0026 double avv_num = 0., avv_den = 0.;
0027 double Etotal = 0.0;
0028
0029 for (int ix = 0; ix < f->getNX(); ix++)
0030 for (int iy = 0; iy < f->getNY(); iy++)
0031 for (int iz = 0; iz < f->getNZ(); iz++) {
0032 c = f->getCell(ix, iy, iz);
0033 double x = f->getX(ix);
0034 double y = f->getY(iy);
0035 double eta = f->getZ(iz);
0036
0037 const double _t = 1.0;
0038 const double q = 1.0;
0039 const double r = sqrt(x * x + y * y);
0040 const double _k =
0041 2. * q * q * _t * r / (1. + q * q * _t * _t + q * q * r * r);
0042 e = 4. * q * q / (1. + 2. * q * q * (_t * _t + r * r) +
0043 pow(q * q * (_t * _t - r * r), 2)) /
0044 _t;
0045 if (e < 0.) e = 0.;
0046 e = pow(e, 4. / 3.);
0047 vx = x / (r + 1e-50) * _k;
0048 vy = y / (r + 1e-50) * _k;
0049 nb = nq = 0.0;
0050 vz = 0.0;
0051
0052 avv_num += sqrt(vx * vx + vy * vy) * e;
0053 avv_den += e;
0054
0055 c->setPrimVar(eos, tau, e, nb, nq, 0., vx, vy, vz);
0056 double _p = eos->p(e, nb, nq, 0.);
0057 const double gamma2 = 1.0 / (1.0 - vx * vx - vy * vy - vz * vz);
0058 Etotal +=
0059 ((e + _p) * gamma2 * (cosh(eta) + vz * sinh(eta)) - _p * cosh(eta));
0060 c->saveQprev();
0061
0062 if (e > 0.) c->setAllM(1.);
0063 }
0064 cout << "average initial flow = " << avv_num / avv_den << endl;
0065 cout << "total energy = " << Etotal *f->getDx() * f->getDy() * f->getDz() *
0066 tau << endl;
0067 }
