Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fld.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fld.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 <cmath>
23 #include <algorithm>
24 #include "inc.h"
25 #include "rmn.h"
26 #include "fld.h"
27 #include "cll.h"
28 #include "eos.h"
29 #include "trancoeff.h"
30 
31 #define OUTPI
32 
33 using namespace std;
34 
35 // returns the velocities in cartesian coordinates, fireball rest frame.
36 // Y=longitudinal rapidity of fluid
37 void Fluid::getCMFvariables(Cell *c, double tau, double &e, double &nb,
38  double &nq, double &ns, double &vx, double &vy,
39  double &Y) {
40  double p, vz;
41  c->getPrimVar(eos, tau, e, p, nb, nq, ns, vx, vy, vz);
42  double eta = getZ(c->getZ());
43  // Y = eta + TMath::ATanH(vz) ;
44  Y = eta + 1. / 2. * log((1. + vz) / (1. - vz));
45  vx = vx * cosh(Y - eta) / cosh(Y);
46  vy = vy * cosh(Y - eta) / cosh(Y);
47 }
48 
49 Fluid::Fluid(EoS *_eos, EoS *_eosH, TransportCoeff *_trcoeff, int _nx, int _ny,
50  int _nz, double _minx, double _maxx, double _miny, double _maxy,
51  double _minz, double _maxz, double _dt, double eCrit) {
52  eos = _eos;
53  eosH = _eosH;
54  trcoeff = _trcoeff;
55  nx = _nx;
56  ny = _ny;
57  nz = _nz;
58  minx = _minx;
59  maxx = _maxx;
60  miny = _miny;
61  maxy = _maxy;
62  minz = _minz;
63  maxz = _maxz;
64  dx = (maxx - minx) / (nx - 1);
65  dy = (maxy - miny) / (ny - 1);
66  dz = (maxz - minz) / (nz - 1);
67  dt = _dt;
68 
69  cell = new Cell[nx * ny * nz];
70 
71  cell0 = new Cell;
72  cell0->setPrimVar(eos, 1.0, 0., 0., 0., 0., 0., 0.,
73  0.); // tau is not important here, since *0
74  cell0->setAllM(0.);
75 
76  for (int ix = 0; ix < nx; ix++)
77  for (int iy = 0; iy < ny; iy++)
78  for (int iz = 0; iz < nz; iz++) {
79  getCell(ix, iy, iz)->setPrev(X_, getCell(ix - 1, iy, iz));
80  getCell(ix, iy, iz)->setNext(X_, getCell(ix + 1, iy, iz));
81  getCell(ix, iy, iz)->setPrev(Y_, getCell(ix, iy - 1, iz));
82  getCell(ix, iy, iz)->setNext(Y_, getCell(ix, iy + 1, iz));
83  getCell(ix, iy, iz)->setPrev(Z_, getCell(ix, iy, iz - 1));
84  getCell(ix, iy, iz)->setNext(Z_, getCell(ix, iy, iz + 1));
85  getCell(ix, iy, iz)->setPos(ix, iy, iz);
86  }
87 
88  output_nt = 0;
89  output_nx = 0;
90  output_ny = 0;
91 }
92 
93 Fluid::~Fluid() {
94  delete[] cell;
95  delete cell0;
96 }
97 
98 void Fluid::initOutput(char *dir, int maxstep, double tau0, int cmpr2dOut) {
99  // directory = dir ;
100  compress2dOut = cmpr2dOut;
101  cout << "maxstep=" << maxstep << endl;
102  char command[255];
103  sprintf(command, "mkdir -p %s", dir);
104  int return_mkdir = system(command);
105  cout << "mkdir returns: " << return_mkdir << endl;
106  string outx = dir;
107  outx.append("/outx.dat");
108  string outxvisc = dir;
109  outxvisc.append("/outx.visc.dat");
110  string outyvisc = dir;
111  outyvisc.append("/outy.visc.dat");
112  string outdiagvisc = dir;
113  outdiagvisc.append("/diag.visc.dat");
114  string outy = dir;
115  outy.append("/outy.dat");
116  string outdiag = dir;
117  outdiag.append("/outdiag.dat");
118  string outz = dir;
119  outz.append("/outz.dat");
120  string outaniz = dir;
121  outaniz.append("/out.aniz.dat");
122  string out2d = dir;
123  out2d.append("/out2D.dat");
124  string outfreeze = dir;
125  outfreeze.append("/freezeout.dat");
126  foutx.open(outx.c_str());
127  fouty.open(outy.c_str());
128  foutz.open(outz.c_str());
129  foutdiag.open(outdiag.c_str());
130  fout2d.open(out2d.c_str());
131  foutxvisc.open(outxvisc.c_str());
132  foutyvisc.open(outyvisc.c_str());
133  foutdiagvisc.open(outdiagvisc.c_str());
134  fout_aniz.open(outaniz.c_str());
135  ffreeze.open(outfreeze.c_str());
136  //################################################################
137  // important remark. for correct diagonal output, nx=ny must hold.
138  //################################################################
139  foutxvisc << maxstep + 1 << " " << getNX() << endl;
140  foutyvisc << maxstep + 1 << " " << getNY() << endl;
141  foutdiagvisc << maxstep + 1 << " " << getNX() << endl;
142  foutx << "# " << maxstep + 1 << " " << getNX() << endl;
143  fouty << "# " << maxstep + 1 << " " << getNY() << endl;
144  foutdiag << "# " << maxstep + 1 << " " << getNX() << endl;
145  foutz << "# " << maxstep + 1 << " " << getNZ() << endl;
146  fout2d << " " << maxstep + 1 << " " << (getNX() - 5) + 1 << " "
147  << (getNY() - 5) + 1 << endl;
148  fout2d << tau0 << " " << tau0 + 0.05 * maxstep << " " << getX(2) << " "
149  << getX(getNX() - 3) << " " << getY(2) << " " << getY(getNY() - 3)
150  << endl;
151  outputGnuplot(tau0);
152  fout_aniz << "# tau <<v_T>> e_p e'_p (to compare with SongHeinz)\n";
153 }
154 
155 void Fluid::correctImagCells(void) {
156  double Q[7];
157  // Z
158  for (int ix = 0; ix < nx; ix++)
159  for (int iy = 0; iy < ny; iy++) {
160  // left boundary
161  getCell(ix, iy, 2)->getQ(Q);
162  getCell(ix, iy, 1)->setQ(Q);
163  getCell(ix, iy, 0)->setQ(Q);
164  // right boundary
165  getCell(ix, iy, nz - 3)->getQ(Q);
166  getCell(ix, iy, nz - 2)->setQ(Q);
167  getCell(ix, iy, nz - 1)->setQ(Q);
168  }
169  // Y
170  for (int ix = 0; ix < nx; ix++)
171  for (int iz = 0; iz < nz; iz++) {
172  // left boundary
173  getCell(ix, 2, iz)->getQ(Q);
174  getCell(ix, 1, iz)->setQ(Q);
175  getCell(ix, 0, iz)->setQ(Q);
176  // right boundary
177  getCell(ix, ny - 3, iz)->getQ(Q);
178  getCell(ix, ny - 2, iz)->setQ(Q);
179  getCell(ix, ny - 1, iz)->setQ(Q);
180  }
181  // X
182  for (int iy = 0; iy < ny; iy++)
183  for (int iz = 0; iz < nz; iz++) {
184  // left boundary
185  getCell(2, iy, iz)->getQ(Q);
186  getCell(1, iy, iz)->setQ(Q);
187  getCell(0, iy, iz)->setQ(Q);
188  // right boundary
189  getCell(nx - 3, iy, iz)->getQ(Q);
190  getCell(nx - 2, iy, iz)->setQ(Q);
191  getCell(nx - 1, iy, iz)->setQ(Q);
192  }
193 }
194 
195 void Fluid::correctImagCellsFull(void) {
196  double Q[7], _pi[4][4], _Pi;
197  // Z
198  for (int ix = 0; ix < nx; ix++)
199  for (int iy = 0; iy < ny; iy++) {
200  // left boundary
201  getCell(ix, iy, 2)->getQ(Q);
202  getCell(ix, iy, 1)->setQ(Q);
203  getCell(ix, iy, 0)->setQ(Q);
204  for (int i = 0; i < 4; i++)
205  for (int j = 0; j < 4; j++) _pi[i][j] = getCell(ix, iy, 2)->getpi(i, j);
206  _Pi = getCell(ix, iy, 2)->getPi();
207 
208  for (int i = 0; i < 4; i++)
209  for (int j = 0; j <= i; j++) {
210  getCell(ix, iy, 0)->setpi(i, j, _pi[i][j]);
211  getCell(ix, iy, 1)->setpi(i, j, _pi[i][j]);
212  }
213  getCell(ix, iy, 0)->setPi(_Pi);
214  getCell(ix, iy, 1)->setPi(_Pi);
215  // right boundary
216  getCell(ix, iy, nz - 3)->getQ(Q);
217  getCell(ix, iy, nz - 2)->setQ(Q);
218  getCell(ix, iy, nz - 1)->setQ(Q);
219  for (int i = 0; i < 4; i++)
220  for (int j = 0; j < 4; j++)
221  _pi[i][j] = getCell(ix, iy, nz - 3)->getpi(i, j);
222  _Pi = getCell(ix, iy, nz - 3)->getPi();
223 
224  for (int i = 0; i < 4; i++)
225  for (int j = 0; j <= i; j++) {
226  getCell(ix, iy, nz - 2)->setpi(i, j, _pi[i][j]);
227  getCell(ix, iy, nz - 1)->setpi(i, j, _pi[i][j]);
228  }
229  getCell(ix, iy, nz - 2)->setPi(_Pi);
230  getCell(ix, iy, nz - 1)->setPi(_Pi);
231  }
232  // Y
233  for (int ix = 0; ix < nx; ix++)
234  for (int iz = 0; iz < nz; iz++) {
235  // left boundary
236  getCell(ix, 2, iz)->getQ(Q);
237  getCell(ix, 1, iz)->setQ(Q);
238  getCell(ix, 0, iz)->setQ(Q);
239  for (int i = 0; i < 4; i++)
240  for (int j = 0; j < 4; j++) _pi[i][j] = getCell(ix, 2, iz)->getpi(i, j);
241  _Pi = getCell(ix, 2, iz)->getPi();
242 
243  for (int i = 0; i < 4; i++)
244  for (int j = 0; j <= i; j++) {
245  getCell(ix, 0, iz)->setpi(i, j, _pi[i][j]);
246  getCell(ix, 1, iz)->setpi(i, j, _pi[i][j]);
247  }
248  getCell(ix, 0, iz)->setPi(_Pi);
249  getCell(ix, 1, iz)->setPi(_Pi);
250  // right boundary
251  getCell(ix, ny - 3, iz)->getQ(Q);
252  getCell(ix, ny - 2, iz)->setQ(Q);
253  getCell(ix, ny - 1, iz)->setQ(Q);
254  for (int i = 0; i < 4; i++)
255  for (int j = 0; j < 4; j++)
256  _pi[i][j] = getCell(ix, ny - 3, iz)->getpi(i, j);
257  _Pi = getCell(ix, ny - 3, iz)->getPi();
258 
259  for (int i = 0; i < 4; i++)
260  for (int j = 0; j <= i; j++) {
261  getCell(ix, ny - 2, iz)->setpi(i, j, _pi[i][j]);
262  getCell(ix, ny - 1, iz)->setpi(i, j, _pi[i][j]);
263  }
264  getCell(ix, ny - 2, iz)->setPi(_Pi);
265  getCell(ix, ny - 1, iz)->setPi(_Pi);
266  }
267  // X
268  for (int iy = 0; iy < ny; iy++)
269  for (int iz = 0; iz < nz; iz++) {
270  // left boundary
271  getCell(2, iy, iz)->getQ(Q);
272  getCell(1, iy, iz)->setQ(Q);
273  getCell(0, iy, iz)->setQ(Q);
274  for (int i = 0; i < 4; i++)
275  for (int j = 0; j < 4; j++) _pi[i][j] = getCell(2, iy, iz)->getpi(i, j);
276  _Pi = getCell(2, iy, iz)->getPi();
277 
278  for (int i = 0; i < 4; i++)
279  for (int j = 0; j <= i; j++) {
280  getCell(0, iy, iz)->setpi(i, j, _pi[i][j]);
281  getCell(1, iy, iz)->setpi(i, j, _pi[i][j]);
282  }
283  getCell(0, iy, iz)->setPi(_Pi);
284  getCell(1, iy, iz)->setPi(_Pi);
285  // right boundary
286  getCell(nx - 3, iy, iz)->getQ(Q);
287  getCell(nx - 2, iy, iz)->setQ(Q);
288  getCell(nx - 1, iy, iz)->setQ(Q);
289  for (int i = 0; i < 4; i++)
290  for (int j = 0; j < 4; j++)
291  _pi[i][j] = getCell(nx - 3, iy, iz)->getpi(i, j);
292  _Pi = getCell(nx - 3, iy, iz)->getPi();
293 
294  for (int i = 0; i < 4; i++)
295  for (int j = 0; j <= i; j++) {
296  getCell(nx - 2, iy, iz)->setpi(i, j, _pi[i][j]);
297  getCell(nx - 1, iy, iz)->setpi(i, j, _pi[i][j]);
298  }
299  getCell(nx - 2, iy, iz)->setPi(_Pi);
300  getCell(nx - 1, iy, iz)->setPi(_Pi);
301  }
302 }
303 
304 void Fluid::updateM(double tau, double dt) {
305  for (int ix = 0; ix < getNX(); ix++)
306  for (int iy = 0; iy < getNY(); iy++)
307  for (int iz = 0; iz < getNZ(); iz++) {
308  Cell *c = getCell(ix, iy, iz);
309  c->setDM(X_, 0.);
310  c->setDM(Y_, 0.);
311  c->setDM(Z_, 0.);
312  if (getCell(ix, iy, iz)->getMaxM() < 1.) {
313  if (getCell(ix + 1, iy, iz)->getM(X_) >= 1. ||
314  getCell(ix - 1, iy, iz)->getM(X_) >= 1.)
315  c->setDM(X_, dt / dx);
316  if (getCell(ix, iy + 1, iz)->getM(Y_) >= 1. ||
317  getCell(ix, iy - 1, iz)->getM(Y_) >= 1.)
318  c->setDM(Y_, dt / dy);
319  if (getCell(ix, iy, iz + 1)->getM(Z_) >= 1. ||
320  getCell(ix, iy, iz - 1)->getM(Z_) >= 1.)
321  c->setDM(Z_, dt / dz / tau);
322 
323  if (c->getDM(X_) == 0. && c->getDM(Y_) == 0.) {
324  if (getCell(ix + 1, iy + 1, iz)->getMaxM() >= 1. ||
325  getCell(ix + 1, iy - 1, iz)->getMaxM() >= 1. ||
326  getCell(ix - 1, iy + 1, iz)->getMaxM() >= 1. ||
327  getCell(ix - 1, iy - 1, iz)->getMaxM() >= 1.) {
328  c->setDM(X_, 0.707 * dt / dx);
329  c->setDM(Y_, 0.707 * dt / dy);
330  }
331  }
332  } // if
333  }
334 
335  for (int ix = 0; ix < getNX(); ix++)
336  for (int iy = 0; iy < getNY(); iy++)
337  for (int iz = 0; iz < getNZ(); iz++) {
338  Cell *c = getCell(ix, iy, iz);
339  c->addM(X_, c->getDM(X_));
340  c->addM(Y_, c->getDM(Y_));
341  c->addM(Z_, c->getDM(Z_));
342  }
343 }
344 
345 
346 void Fluid::outputGnuplot(double tau) {
347  double e, p, nb, nq, ns, t, mub, muq, mus, vx, vy, vz;
348 
349  // X direction
350  for (int ix = 0; ix < nx; ix++) {
351  double x = getX(ix);
352  Cell *c = getCell(ix, ny / 2, nz / 2);
353  getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
354  eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
355  foutx << setw(14) << tau << setw(14) << x << setw(14) << vx << setw(14)
356  << vy << setw(14) << e << setw(14) << nb << setw(14) << t << setw(14)
357  << mub;
358  foutx << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
359  << setw(14) << c->getpi(0, 2);
360  foutx << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
361  << setw(14) << c->getpi(1, 2);
362  foutx << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
363  << setw(14) << c->getpi(2, 3);
364  foutx << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
365  << c->getViscCorrCutFlag() << endl;
366  }
367  foutx << endl;
368 
369  // Y direction
370  for (int iy = 0; iy < ny; iy++) {
371  double y = getY(iy);
372  Cell *c = getCell(nx / 2, iy, nz / 2);
373  getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
374  eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
375  fouty << setw(14) << tau << setw(14) << y << setw(14) << vy << setw(14)
376  << vx << setw(14) << e << setw(14) << nb << setw(14) << t << setw(14)
377  << mub;
378  fouty << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
379  << setw(14) << c->getpi(0, 2);
380  fouty << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
381  << setw(14) << c->getpi(1, 2);
382  fouty << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
383  << setw(14) << c->getpi(2, 3);
384  fouty << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
385  << c->getViscCorrCutFlag() << endl;
386  }
387  fouty << endl;
388 
389  // diagonal
390  for (int ix = 0; ix < nx; ix++) {
391  double x = getY(ix);
392  Cell *c = getCell(ix, ix, nz / 2);
393  getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
394  eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
395  foutdiag << setw(14) << tau << setw(14) << sqrt(2.) * x << setw(14) << vx
396  << setw(14) << vy << setw(14) << e << setw(14) << nb << setw(14)
397  << t << setw(14) << mub << endl;
398  foutdiag << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
399  << setw(14) << c->getpi(0, 2);
400  foutdiag << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
401  << setw(14) << c->getpi(1, 2);
402  foutdiag << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
403  << setw(14) << c->getpi(2, 3);
404  foutdiag << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
405  << c->getViscCorrCutFlag() << endl;
406  }
407  foutdiag << endl;
408 
409  // Z direction
410  for (int iz = 0; iz < nz; iz++) {
411  double z = getZ(iz);
412  Cell *c = getCell(nx / 2, ny / 2, iz);
413  getCMFvariables(getCell(nx / 2, ny / 2, iz), tau, e, nb, nq, ns, vx, vy,
414  vz);
415  eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
416  foutz << setw(14) << tau << setw(14) << z << setw(14) << vz << setw(14)
417  << vx << setw(14) << e << setw(14) << nb << setw(14) << t << setw(14)
418  << mub;
419  foutz << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
420  << setw(14) << c->getpi(0, 2);
421  foutz << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
422  << setw(14) << c->getpi(1, 2);
423  foutz << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
424  << setw(14) << c->getpi(2, 3);
425  foutz << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
426  << c->getViscCorrCutFlag() << endl;
427  }
428  foutz << endl;
429 
430  // averaged quantities for y=0
431  double T0xx = 0.0, T0yy = 0.0, Txx = 0.0, Tyy = 0.0, vtNum = 0.0, vtDen = 0.0;
432  for (int ix = 0; ix < nx; ix++)
433  for (int iy = 0; iy < ny; iy++) {
434  Cell *c = getCell(ix, iy, nz / 2);
435  getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
436  eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
437  double dTxx =
438  (e + p) * vx * vx / (1.0 - vx * vx - vy * vy - pow(tanh(vz), 2)) + p;
439  double dTyy =
440  (e + p) * vy * vy / (1.0 - vx * vx - vy * vy - pow(tanh(vz), 2)) + p;
441  T0xx += dTxx;
442  T0yy += dTyy;
443  Txx += dTxx + c->getpi(1, 1);
444  Tyy += dTyy + c->getpi(2, 2);
445  vtNum += sqrt(vx * vx + vy * vy) * e /
446  sqrt(1.0 - vx * vx - vy * vy - pow(tanh(vz), 2));
447  vtDen += e / sqrt(1.0 - vx * vx - vy * vy - pow(tanh(vz), 2));
448  }
449  fout_aniz << setw(14) << tau << setw(14) << vtNum / vtDen << setw(14)
450  << (T0xx - T0yy) / (T0xx + T0yy) << setw(14)
451  << (Txx - Tyy) / (Txx + Tyy) << endl;
452 }
453 
454 // unput: geom. rapidity + velocities in Bjorken frame, --> output: velocities
455 // in lab.frame
456 void transformToLab(double eta, double &vx, double &vy, double &vz) {
457  const double Y = eta + 1. / 2. * log((1. + vz) / (1. - vz));
458  vx = vx * cosh(Y - eta) / cosh(Y);
459  vy = vy * cosh(Y - eta) / cosh(Y);
460  vz = tanh(Y);
461 }
462 
463 void Fluid::calcTotals(double tau) {
464  double e, p, nb, nq, ns, t, mub, muq, mus, vx, vy, vz, Q[7];
465  double E = 0., Efull = 0., Px = 0., Nb1 = 0., Nb2 = 0., S = 0.;
466  double eta = 0;
467 
468  fout2d << endl;
469  for (int ix = 2; ix < nx - 2; ix++)
470  for (int iy = 2; iy < ny - 2; iy++)
471  for (int iz = 2; iz < nz - 2; iz++) {
472  Cell *c = getCell(ix, iy, iz);
473  getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
474  c->getQ(Q);
475  eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
476  double s = eos->s(e, nb, nq, ns);
477  eta = getZ(iz);
478  E += tau * (e + p) / (1. - vx * vx - vy * vy - tanh(vz) * tanh(vz)) *
479  (cosh(eta) - tanh(vz) * sinh(eta)) -
480  tau * p * cosh(eta);
481  Nb1 += Q[NB_];
482  Nb2 += tau * nb * (cosh(eta) - tanh(vz) * sinh(eta)) /
483  sqrt(1. - vx * vx - vy * vy - tanh(vz) * tanh(vz));
484  // -- noneq. corrections to entropy flux
485  const double gmumu[4] = {1., -1., -1., -1.};
486  double deltas = 0.;
487  for (int i = 0; i < 4; i++)
488  for (int j = 0; j < 4; j++)
489  deltas += pow(c->getpi(i, j), 2) * gmumu[i] * gmumu[j];
490  if (t > 0.05) s += 1.5 * deltas / ((e + p) * t);
491  S += tau * s * (cosh(eta) - tanh(vz) * sinh(eta)) /
492  sqrt(1. - vx * vx - vy * vy - tanh(vz) * tanh(vz));
493  Efull +=
494  tau * (e + p) / (1. - vx * vx - vy * vy - tanh(vz) * tanh(vz)) *
495  (cosh(eta) - tanh(vz) * sinh(eta)) -
496  tau * p * cosh(eta);
497  if (trcoeff->isViscous())
498  Efull += tau * c->getpi(0, 0) * cosh(eta) +
499  tau * c->getpi(0, 3) * sinh(eta);
500  }
501  E = E * dx * dy * dz;
502  Efull = Efull * dx * dy * dz;
503  Nb1 *= dx * dy * dz;
504  Nb2 *= dx * dy * dz;
505  S *= dx * dy * dz;
506  cout << setw(16) << "calcTotals: E = " << setw(14) << E
507  << " Efull = " << setw(14) << Efull << endl;
508  cout << setw(16) << "Px = " << setw(14) << Px << " S = " << setw(14) << S
509  << endl;
510 }