Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hdo.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hdo.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 <algorithm>
23 #include <unistd.h>
24 #include "hdo.h"
25 #include "inc.h"
26 #include "rmn.h"
27 #include "fld.h"
28 #include "eos.h"
29 #include "cll.h"
30 #include "trancoeff.h"
31 
32 using namespace std;
33 
34 extern bool debugRiemann;
35 
36 double sign(double x) {
37  if (x > 0)
38  return 1.;
39  else if (x < 0.)
40  return -1.;
41  else
42  return 0.;
43 }
44 
45 // this version contains NO PRE-ADVECTION for the IS solution
46 
47 // enable this to use formal solution for the relaxation part of
48 // Israel-Stewart equations (not recommended)
49 //#define FORMAL_SOLUTION
50 // else: use 1st order finite difference update
51 
52 Hydro::Hydro(Fluid *_f, EoS *_eos, TransportCoeff *_trcoeff, double _t0,
53  double _dt) {
54  eos = _eos;
55  trcoeff = _trcoeff;
56  f = _f;
57  dt = _dt;
58  tau = _t0;
59 }
60 
61 Hydro::~Hydro() {}
62 
63 void Hydro::setDtau(double deltaTau) {
64  dt = deltaTau;
65  if (dt > f->getDx() / 2. ||
66  dt > f->getDy() / 2. /*|| dt>tau*f->getDz()/2. */) {
67  cout << "too big delta_tau " << dt << " " << f->getDx() << " "
68  << f->getDy() << " " << tau * f->getDz() << endl;
69  exit(1);
70  }
71 }
72 
73 void Hydro::hlle_flux(Cell *left, Cell *right, int direction, int mode) {
74  // for all variables, suffix "l" = left state, "r" = right state
75  // with respect to the cell boundary
76  double el, er, pl, pr, nbl, nql, nsl, nbr, nqr, nsr, vxl, vxr, vyl, vyr, vzl,
77  vzr, bl = 0., br = 0., csb, vb, El, Er, dx = 0.;
78  double Ftl = 0., Fxl = 0., Fyl = 0., Fzl = 0., Fbl = 0., Fql = 0., Fsl = 0.,
79  Ftr = 0., Fxr = 0., Fyr = 0., Fzr = 0., Fbr = 0., Fqr = 0., Fsr = 0.;
80  double U1l, U2l, U3l, U4l, Ubl, Uql, Usl, U1r, U2r, U3r, U4r, Ubr, Uqr, Usr;
81  double flux[7];
82  const double dta = mode == 0 ? dt / 2. : dt;
83  double tauFactor; // fluxes are also multiplied by tau
84  if (mode == PREDICT) {
85  // get primitive quantities from Q_{i+} at previous timestep
86  left->getPrimVarRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
87  direction);
88  // ... and Q_{(i+1)-}
89  right->getPrimVarLeft(eos, tau, er, pr, nbr, nqr, nsr, vxr, vyr, vzr,
90  direction);
91  El = (el + pl) / (1 - vxl * vxl - vyl * vyl - vzl * vzl);
92  Er = (er + pr) / (1 - vxr * vxr - vyr * vyr - vzr * vzr);
93  tauFactor = tau;
94  } else {
95  // use half-step updated Q's for corrector step
96  left->getPrimVarHRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
97  direction);
98  right->getPrimVarHLeft(eos, tau, er, pr, nbr, nqr, nsr, vxr, vyr, vzr,
99  direction);
100  El = (el + pl) / (1 - vxl * vxl - vyl * vyl - vzl * vzl);
101  Er = (er + pr) / (1 - vxr * vxr - vyr * vyr - vzr * vzr);
102  tauFactor = tau + dt / 2.;
103  }
104 
105  if (el < 0.) {
106  el = 0.;
107  pl = 0.;
108  }
109  if (er < 0.) {
110  er = 0.;
111  pr = 0.;
112  }
113 
114  if (el > 1e10) {
115  cout << "e>1e10; debug info below:\n";
116  left->Dump(tau);
117  debugRiemann = true;
118  if (mode == PREDICT)
119  left->getPrimVarRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
120  direction);
121  else
122  left->getPrimVarHRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
123  direction);
124  debugRiemann = false;
125  exit(0);
126  }
127 
128  // skip the procedure for two empty cells
129  if (el == 0. && er == 0.) return;
130  if (pr < 0.) {
131  cout << "Negative pressure" << endl;
132  left->getPrimVarRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
133  direction);
134  right->getPrimVarLeft(eos, tau, er, pr, nbr, nqr, nsr, vxr, vyr, vzr,
135  direction);
136  }
137 
138  // skip the procedure for two partially vacuum cells
139  if (left->getM(direction) < 1. && right->getM(direction) < 1.) return;
140 
141  double gammal = 1. / sqrt(1 - vxl * vxl - vyl * vyl - vzl * vzl);
142  double gammar = 1. / sqrt(1 - vxr * vxr - vyr * vyr - vzr * vzr);
143  U1l = gammal * gammal * (el + pl) * vxl;
144  U2l = gammal * gammal * (el + pl) * vyl;
145  U3l = gammal * gammal * (el + pl) * vzl;
146  U4l = gammal * gammal * (el + pl) - pl;
147  Ubl = gammal * nbl;
148  Uql = gammal * nql;
149  Usl = gammal * nsl;
150 
151  U1r = gammar * gammar * (er + pr) * vxr;
152  U2r = gammar * gammar * (er + pr) * vyr;
153  U3r = gammar * gammar * (er + pr) * vzr;
154  U4r = gammar * gammar * (er + pr) - pr;
155  Ubr = gammar * nbr;
156  Uqr = gammar * nqr;
157  Usr = gammar * nsr;
158 
159  if (direction == X_) {
160  Ftl = U4l * vxl + pl * vxl;
161  Fxl = U1l * vxl + pl;
162  Fyl = U2l * vxl;
163  Fzl = U3l * vxl;
164  Fbl = Ubl * vxl;
165  Fql = Uql * vxl;
166  Fsl = Usl * vxl;
167 
168  Ftr = U4r * vxr + pr * vxr;
169  Fxr = U1r * vxr + pr;
170  Fyr = U2r * vxr;
171  Fzr = U3r * vxr;
172  Fbr = Ubr * vxr;
173  Fqr = Uqr * vxr;
174  Fsr = Usr * vxr;
175 
176  // for the case of constant c_s only
177  csb = sqrt(eos->cs2() + 0.5 * sqrt(El * Er) / pow(sqrt(El) + sqrt(Er), 2) *
178  pow(vxl - vxr, 2));
179  vb = (sqrt(El) * vxl + sqrt(Er) * vxr) / (sqrt(El) + sqrt(Er));
180  bl = min(0., min((vb - csb) / (1 - vb * csb),
181  (vxl - eos->cs()) / (1 - vxl * eos->cs())));
182  br = max(0., max((vb + csb) / (1 + vb * csb),
183  (vxr + eos->cs()) / (1 + vxr * eos->cs())));
184 
185  dx = f->getDx();
186 
187  // bl or br in the case of boundary with vacuum
188  if (el == 0.) bl = -1.;
189  if (er == 0.) br = 1.;
190  }
191  if (direction == Y_) {
192  Ftl = U4l * vyl + pl * vyl;
193  Fxl = U1l * vyl;
194  Fyl = U2l * vyl + pl;
195  Fzl = U3l * vyl;
196  Fbl = Ubl * vyl;
197  Fql = Uql * vyl;
198  Fsl = Usl * vyl;
199 
200  Ftr = U4r * vyr + pr * vyr;
201  Fxr = U1r * vyr;
202  Fyr = U2r * vyr + pr;
203  Fzr = U3r * vyr;
204  Fbr = Ubr * vyr;
205  Fqr = Uqr * vyr;
206  Fsr = Usr * vyr;
207 
208  // for the case of constant c_s only
209  csb = sqrt(eos->cs2() + 0.5 * sqrt(El * Er) / pow(sqrt(El) + sqrt(Er), 2) *
210  pow(vyl - vyr, 2));
211  vb = (sqrt(El) * vyl + sqrt(Er) * vyr) / (sqrt(El) + sqrt(Er));
212  bl = min(0., min((vb - csb) / (1 - vb * csb),
213  (vyl - eos->cs()) / (1 - vyl * eos->cs())));
214  br = max(0., max((vb + csb) / (1 + vb * csb),
215  (vyr + eos->cs()) / (1 + vyr * eos->cs())));
216 
217  dx = f->getDy();
218 
219  // bl or br in the case of boundary with vacuum
220  if (el == 0.) bl = -1.;
221  if (er == 0.) br = 1.;
222  }
223  if (direction == Z_) {
224  double tau1 = tau_z;
225  Ftl = U4l * vzl / tau1 + pl * vzl / tau1;
226  Fxl = U1l * vzl / tau1;
227  Fyl = U2l * vzl / tau1;
228  Fzl = U3l * vzl / tau1 + pl / tau1;
229  Fbl = Ubl * vzl / tau1;
230  Fql = Uql * vzl / tau1;
231  Fsl = Usl * vzl / tau1;
232 
233  Ftr = U4r * vzr / tau1 + pr * vzr / tau1;
234  Fxr = U1r * vzr / tau1;
235  Fyr = U2r * vzr / tau1;
236  Fzr = U3r * vzr / tau1 + pr / tau1;
237  Fbr = Ubr * vzr / tau1;
238  Fqr = Uqr * vzr / tau1;
239  Fsr = Usr * vzr / tau1;
240 
241  // for the case of constant c_s only
242  // factor 1/tau accounts for eta-coordinate
243 
244  // different estimate
245  csb = sqrt(eos->cs2() + 0.5 * sqrt(El * Er) / pow(sqrt(El) + sqrt(Er), 2) *
246  pow(vzl - vzr, 2));
247  vb = (sqrt(El) * vzl + sqrt(Er) * vzr) / (sqrt(El) + sqrt(Er));
248  bl = 1. / tau * min(0., min((vb - csb) / (1 - vb * csb),
249  (vzl - eos->cs()) / (1 - vzl * eos->cs())));
250  br = 1. / tau * max(0., max((vb + csb) / (1 + vb * csb),
251  (vzr + eos->cs()) / (1 + vzr * eos->cs())));
252 
253  dx = f->getDz();
254 
255  // bl or br in the case of boundary with vacuum
256  if (el == 0.) bl = -1. / tau;
257  if (er == 0.) br = 1. / tau;
258  }
259 
260  if (bl == 0. && br == 0.) return;
261 
262  // finally, HLLE formula for the fluxes
263  flux[T_] = tauFactor * dta / dx *
264  (-bl * br * (U4l - U4r) + br * Ftl - bl * Ftr) / (-bl + br);
265  flux[X_] = tauFactor * dta / dx *
266  (-bl * br * (U1l - U1r) + br * Fxl - bl * Fxr) / (-bl + br);
267  flux[Y_] = tauFactor * dta / dx *
268  (-bl * br * (U2l - U2r) + br * Fyl - bl * Fyr) / (-bl + br);
269  flux[Z_] = tauFactor * dta / dx *
270  (-bl * br * (U3l - U3r) + br * Fzl - bl * Fzr) / (-bl + br);
271  flux[NB_] = tauFactor * dta / dx *
272  (-bl * br * (Ubl - Ubr) + br * Fbl - bl * Fbr) / (-bl + br);
273  flux[NQ_] = tauFactor * dta / dx *
274  (-bl * br * (Uql - Uqr) + br * Fql - bl * Fqr) / (-bl + br);
275  flux[NS_] = tauFactor * dta / dx *
276  (-bl * br * (Usl - Usr) + br * Fsl - bl * Fsr) / (-bl + br);
277 
278  if (flux[NB_] != flux[NB_]) { // if things failed
279  cout << "---- error in hlle_flux: f_nb undefined!\n";
280  cout << setw(12) << U4l << setw(12) << U1l << setw(12) << U2l << setw(12)
281  << U3l << endl;
282  cout << setw(12) << U4r << setw(12) << U1r << setw(12) << U2r << setw(12)
283  << U3r << endl;
284  cout << setw(12) << Ubl << setw(12) << Uql << setw(12) << Usl << endl;
285  cout << setw(12) << Ubr << setw(12) << Uqr << setw(12) << Usr << endl;
286  cout << setw(12) << Ftl << setw(12) << Fxl << setw(12) << Fyl << setw(12)
287  << Fzl << endl;
288  cout << setw(12) << Ftr << setw(12) << Fxr << setw(12) << Fyr << setw(12)
289  << Fzr << endl;
290  exit(1);
291  }
292 
293  // update the cumulative fluxes in both neighbouring cells
294  left->addFlux(-flux[T_], -flux[X_], -flux[Y_], -flux[Z_], -flux[NB_],
295  -flux[NQ_], -flux[NS_]);
296  right->addFlux(flux[T_], flux[X_], flux[Y_], flux[Z_], flux[NB_], flux[NQ_],
297  flux[NS_]);
298 }
299 
300 void Hydro::source(double tau1, double x, double y, double z, double e,
301  double p, double nb, double nq, double ns, double vx,
302  double vy, double vz, double S[7]) {
303  double Q[7];
304  transformCV(e, p, nb, nq, ns, vx, vy, vz, Q);
305  S[T_] = -Q[T_] * vz * vz - p * (1. + vz * vz);
306  S[X_] = 0.;
307  S[Y_] = 0.;
308  S[Z_] = -Q[Z_];
309  S[NB_] = 0.;
310  S[NQ_] = 0.;
311  S[NS_] = 0.;
312 }
313 
314 void Hydro::source_step(int ix, int iy, int iz, int mode) {
315  double _dt;
316  if (mode == PREDICT)
317  _dt = dt / 2.;
318  else
319  _dt = dt;
320 
321  double e, p, nb, nq, ns, vx, vy, vz;
322  double k[7];
323 
324  double x = f->getX(ix), y = f->getY(iy), z = f->getZ(iz);
325  Cell *c = f->getCell(ix, iy, iz);
326 
327  if (mode == PREDICT) {
328  c->getPrimVar(eos, tau, e, p, nb, nq, ns, vx, vy, vz);
329  } else {
330  c->getPrimVarHCenter(eos, tau + dt / 2., e, p, nb, nq, ns, vx, vy, vz);
331  }
332 
333  source(tau, x, y, z, e, p, nb, nq, ns, vx, vy, vz, k); // setting k1
334  for (int i = 0; i < 7; i++) k[i] *= _dt;
335 
336  if (k[NB_] != k[NB_]) { // something failed
337  cout << "---- error in source_step: k_nb undefined!\n";
338  cout << setw(12) << k[0] << setw(12) << k[1] << setw(12) << k[2] << setw(12)
339  << k[3] << endl;
340  cout << setw(12) << k[4] << setw(12) << k[5] << setw(12) << k[6] << endl;
341  exit(1);
342  }
343  c->addFlux(k[T_], k[X_], k[Y_], k[Z_], k[NB_], k[NQ_], k[NS_]);
344 }
345 
346 void Hydro::visc_source_step(int ix, int iy, int iz) {
347  double e, p, nb, nq, ns, vx, vy, vz;
348  double uuu[4];
349  double k[7];
350 
351  Cell *c = f->getCell(ix, iy, iz);
352 
353  c->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vx, vy, vz);
354  if (e <= 0.) return;
355  uuu[0] = 1. / sqrt(1. - vx * vx - vy * vy - vz * vz);
356  uuu[1] = uuu[0] * vx;
357  uuu[2] = uuu[0] * vy;
358  uuu[3] = uuu[0] * vz;
359 
360  k[T_] = - c->getpiH(3, 3) +
361  c->getPiH() * ( - 1.0 - uuu[3] * uuu[3]);
362  k[X_] = 0.;
363  k[Y_] = 0.;
364  k[Z_] = - (c->getpiH(0, 3) + c->getPiH() * uuu[0] * uuu[3]);
365  for (int i = 0; i < 4; i++) k[i] *= dt;
366 
367  c->addFlux(k[T_], k[X_], k[Y_], k[Z_], 0., 0., 0.);
368 }
369 
370 // for the procedure below, the following approximations are used:
371 // dv/d(tau) = v^{t+dt}_ideal - v^{t}
372 // dv/dx_i ~ v^{x+dx}-v{x-dx},
373 // which makes sense after non-viscous step
374 void Hydro::NSquant(int ix, int iy, int iz, double pi[4][4], double &Pi,
375  double dmu[4][4], double &du) {
376  const double VMIN = 1e-2;
377  const double UDIFF = 3.0;
378  double e0, e1, p, nb, nq, ns, vx1, vy1, vz1, vx0, vy0, vz0, vxH, vyH, vzH;
379  double ut0, ux0, uy0, uz0, ut1, ux1, uy1, uz1;
380  // double dmu [4][4] ; // \partial_\mu u^\nu matrix
381  // coordinates: 0=tau, 1=x, 2=y, 3=eta
382  double Z[4][4][4][4]; // Z[mu][nu][lambda][rho]
383  double uuu[4]; // the 4-velocity
384  double gmunu[4][4] = {{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0},
385  {0, 0, 0, -1}}; // omit 1/tau^2 in g^{eta,eta}
386  Cell *c = f->getCell(ix, iy, iz);
387  double dx = f->getDx(), dy = f->getDy(), dz = f->getDz();
388  // check if the cell is next to vacuum from +-x, +-y side:
389  if (c->getNext(X_)->getMaxM() <= 0.9 || c->getNext(Y_)->getMaxM() <= 0.9 ||
390  c->getPrev(X_)->getMaxM() <= 0.9 || c->getPrev(Y_)->getMaxM() <= 0.9 ||
391  f->getCell(ix + 1, iy + 1, iz)->getMaxM() <= 0.9 ||
392  f->getCell(ix + 1, iy - 1, iz)->getMaxM() <= 0.9 ||
393  f->getCell(ix - 1, iy + 1, iz)->getMaxM() <= 0.9 ||
394  f->getCell(ix - 1, iy - 1, iz)->getMaxM() <= 0.9) {
395  for (int i = 0; i < 4; i++)
396  for (int j = 0; j < 4; j++) {
397  pi[i][j] = 0.;
398  dmu[i][j] = 0.;
399  }
400  Pi = du = 0.;
401  return;
402  }
403  // calculation of \partial_\mu u^\nu matrix
404  // mu=first index, nu=second index
405  // centered differences with respect to the values at (it+1/2, ix, iy, iz)
406  // d_tau u^\mu
407  c->getPrimVarPrev(eos, tau - dt, e0, p, nb, nq, ns, vx0, vy0, vz0);
408  c->getPrimVar(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
409  c->getPrimVarHCenter(eos, tau - 0.5 * dt, e1, p, nb, nq, ns, vxH, vyH, vzH);
410  //############## get transport coefficients
411  double T, mub, muq, mus;
412  double etaS, zetaS;
413  double s = eos->s(e1, nb, nq, ns); // entropy density in the current cell
414  eos->eos(e1, nb, nq, ns, T, mub, muq, mus, p);
415  trcoeff->getEta(e1, T, etaS, zetaS);
416  //##############
417  // if(e1<0.00004) s=0. ; // negative pressure due to pi^zz for small e
418  ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
419  ux0 = ut0 * vx0;
420  uy0 = ut0 * vy0;
421  uz0 = ut0 * vz0;
422  ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
423  ux1 = ut1 * vx1;
424  uy1 = ut1 * vy1;
425  uz1 = ut1 * vz1;
426  uuu[0] = 1.0 / sqrt(1.0 - vxH * vxH - vyH * vyH - vzH * vzH);
427  uuu[1] = uuu[0] * vxH;
428  uuu[2] = uuu[0] * vyH;
429  uuu[3] = uuu[0] * vzH;
430 
431  dmu[0][0] = (ut1 * ut1 - ut0 * ut0) / 2. / uuu[0] / dt;
432  dmu[0][1] = (ux1 * ux1 - ux0 * ux0) / 2. / uuu[1] / dt;
433  dmu[0][2] = (uy1 * uy1 - uy0 * uy0) / 2. / uuu[2] / dt;
434  dmu[0][3] = (uz1 * uz1 - uz0 * uz0) / 2. / uuu[3] / dt;
435  if (fabs(0.5 * (ut1 + ut0) / ut1) > UDIFF) dmu[0][0] = (ut1 - ut0) / dt;
436  if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / ux1) > UDIFF)
437  dmu[0][1] = (ux1 - ux0) / dt;
438  if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uy1) > UDIFF)
439  dmu[0][2] = (uy1 - uy0) / dt;
440  if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uz1) > UDIFF)
441  dmu[0][3] = (uz1 - uz0) / dt;
442  if (e1 <= 0. || e0 <= 0.) { // matter-vacuum
443  dmu[0][0] = dmu[0][1] = dmu[0][2] = dmu[0][3] = 0.;
444  }
445  // d_x u^\mu
446  f->getCell(ix + 1, iy, iz)
447  ->getPrimVarHCenter(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
448  f->getCell(ix - 1, iy, iz)
449  ->getPrimVarHCenter(eos, tau, e0, p, nb, nq, ns, vx0, vy0, vz0);
450  if (e1 > 0. && e0 > 0.) {
451  ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
452  ux0 = ut0 * vx0;
453  uy0 = ut0 * vy0;
454  uz0 = ut0 * vz0;
455  ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
456  ux1 = ut1 * vx1;
457  uy1 = ut1 * vy1;
458  uz1 = ut1 * vz1;
459  dmu[1][0] = 0.25 * (ut1 * ut1 - ut0 * ut0) / uuu[0] / dx;
460  dmu[1][1] = 0.25 * (ux1 * ux1 - ux0 * ux0) / uuu[1] / dx;
461  dmu[1][2] = 0.25 * (uy1 * uy1 - uy0 * uy0) / uuu[2] / dx;
462  dmu[1][3] = 0.25 * (uz1 * uz1 - uz0 * uz0) / uuu[3] / dx;
463  if (fabs(0.5 * (ut1 + ut0) / uuu[0]) > UDIFF)
464  dmu[1][0] = 0.5 * (ut1 - ut0) / dx;
465  if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / uuu[1]) > UDIFF)
466  dmu[1][1] = 0.5 * (ux1 - ux0) / dx;
467  if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uuu[2]) > UDIFF)
468  dmu[1][2] = 0.5 * (uy1 - uy0) / dx;
469  if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uuu[3]) > UDIFF)
470  dmu[1][3] = 0.5 * (uz1 - uz0) / dx;
471  } else { // matter-vacuum
472  dmu[1][0] = dmu[1][1] = dmu[1][2] = dmu[1][3] = 0.;
473  }
474  if (fabs(dmu[1][3]) > 1e+10)
475  cout << "dmu[1][3]: " << uz1 << " " << uz0 << " " << uuu[3] << endl;
476  // d_y u^\mu
477  f->getCell(ix, iy + 1, iz)
478  ->getPrimVarHCenter(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
479  f->getCell(ix, iy - 1, iz)
480  ->getPrimVarHCenter(eos, tau, e0, p, nb, nq, ns, vx0, vy0, vz0);
481  if (e1 > 0. && e0 > 0.) {
482  ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
483  ux0 = ut0 * vx0;
484  uy0 = ut0 * vy0;
485  uz0 = ut0 * vz0;
486  ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
487  ux1 = ut1 * vx1;
488  uy1 = ut1 * vy1;
489  uz1 = ut1 * vz1;
490  dmu[2][0] = 0.25 * (ut1 * ut1 - ut0 * ut0) / uuu[0] / dy;
491  dmu[2][1] = 0.25 * (ux1 * ux1 - ux0 * ux0) / uuu[1] / dy;
492  dmu[2][2] = 0.25 * (uy1 * uy1 - uy0 * uy0) / uuu[2] / dy;
493  dmu[2][3] = 0.25 * (uz1 * uz1 - uz0 * uz0) / uuu[3] / dy;
494  if (fabs(0.5 * (ut1 + ut0) / uuu[0]) > UDIFF)
495  dmu[2][0] = 0.5 * (ut1 - ut0) / dy;
496  if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / uuu[1]) > UDIFF)
497  dmu[2][1] = 0.5 * (ux1 - ux0) / dy;
498  if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uuu[2]) > UDIFF)
499  dmu[2][2] = 0.5 * (uy1 - uy0) / dy;
500  if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uuu[3]) > UDIFF)
501  dmu[2][3] = 0.5 * (uz1 - uz0) / dy;
502  } else { // matter-vacuum
503  dmu[2][0] = dmu[2][1] = dmu[2][2] = dmu[2][3] = 0.;
504  }
505  // d_z u^\mu
506  f->getCell(ix, iy, iz + 1)
507  ->getPrimVarHCenter(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
508  f->getCell(ix, iy, iz - 1)
509  ->getPrimVarHCenter(eos, tau, e0, p, nb, nq, ns, vx0, vy0, vz0);
510  if (e1 > 0. && e0 > 0.) {
511  ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
512  ux0 = ut0 * vx0;
513  uy0 = ut0 * vy0;
514  uz0 = ut0 * vz0;
515  ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
516  ux1 = ut1 * vx1;
517  uy1 = ut1 * vy1;
518  uz1 = ut1 * vz1;
519  dmu[3][0] = 0.25 * (ut1 * ut1 - ut0 * ut0) / uuu[0] / dz / (tau + 0.5 * dt);
520  dmu[3][1] = 0.25 * (ux1 * ux1 - ux0 * ux0) / uuu[1] / dz / (tau + 0.5 * dt);
521  dmu[3][2] = 0.25 * (uy1 * uy1 - uy0 * uy0) / uuu[2] / dz / (tau + 0.5 * dt);
522  dmu[3][3] = 0.25 * (uz1 * uz1 - uz0 * uz0) / uuu[3] / dz / (tau + 0.5 * dt);
523  if (fabs(0.5 * (ut1 + ut0) / uuu[0]) > UDIFF)
524  dmu[3][0] = 0.5 * (ut1 - ut0) / dz / (tau + 0.5 * dt);
525  if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / uuu[1]) > UDIFF)
526  dmu[3][1] = 0.5 * (ux1 - ux0) / dz / (tau + 0.5 * dt);
527  if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uuu[2]) > UDIFF)
528  dmu[3][2] = 0.5 * (uy1 - uy0) / dz / (tau + 0.5 * dt);
529  if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uuu[3]) > UDIFF)
530  dmu[3][3] = 0.5 * (uz1 - uz0) / dz / (tau + 0.5 * dt);
531  } else { // matter-vacuum
532  dmu[3][0] = dmu[3][1] = dmu[3][2] = dmu[3][3] = 0.;
533  }
534  // additional terms from Christoffel symbols :)
535  dmu[3][0] += uuu[3] / (tau - 0.5 * dt);
536  dmu[3][3] += uuu[0] / (tau - 0.5 * dt);
537  // calculation of Z[mu][nu][lambda][rho]
538  for (int i = 0; i < 4; i++)
539  for (int j = 0; j < 4; j++)
540  for (int k = 0; k < 4; k++)
541  for (int l = 0; l < 4; l++) Z[i][j][k][l] = 0.0;
542  // filling Z matrix
543  for (int mu = 0; mu < 4; mu++)
544  for (int nu = 0; nu < 4; nu++)
545  for (int lam = 0; lam < 4; lam++)
546  for (int rho = 0; rho < 4; rho++) {
547  if (nu == rho)
548  Z[mu][nu][lam][rho] += 0.5 * (gmunu[mu][lam] - uuu[mu] * uuu[lam]);
549  if (mu == rho)
550  Z[mu][nu][lam][rho] += 0.5 * (gmunu[nu][lam] - uuu[nu] * uuu[lam]);
551  if (lam == rho)
552  Z[mu][nu][lam][rho] -= (gmunu[mu][nu] - uuu[mu] * uuu[nu]) / 3.0;
553  }
554  // calculating sigma[mu][nu]
555  for (int i = 0; i < 4; i++)
556  for (int j = 0; j < 4; j++) {
557  pi[i][j] = 0.0;
558  for (int k = 0; k < 4; k++)
559  for (int l = 0; l < 4; l++) {
560  pi[i][j] += Z[i][j][k][l] * dmu[k][l] * 2.0 * etaS * s / 5.068;
561  }
562  }
563  Pi = -zetaS * s * (dmu[0][0] + dmu[1][1] + dmu[2][2] + dmu[3][3]) /
564  5.068; // fm^{-4} --> GeV/fm^3
565  du = dmu[0][0] + dmu[1][1] + dmu[2][2] + dmu[3][3];
566  //--------- debug part: NaN/inf check, trace check, diag check, transversality
567  //check
568  for (int i = 0; i < 4; i++)
569  for (int j = 0; j < 4; j++) {
570  if (pi[i][j] != 0. && fabs(1.0 - pi[j][i] / pi[i][j]) > 1e-10)
571  cout << "non-diag: " << pi[i][j] << " " << pi[j][i] << endl;
572  if (isinf(pi[i][j]) || isnan(pi[i][j])) {
573  cout << "hydro:NSquant: inf/nan i " << i << " j " << j << endl;
574  exit(1);
575  }
576  }
577 }
578 
579 void Hydro::setNSvalues() {
580  double e, p, nb, nq, ns, vx, vy, vz, piNS[4][4], PiNS;
581  for (int ix = 0; ix < f->getNX(); ix++)
582  for (int iy = 0; iy < f->getNY(); iy++)
583  for (int iz = 0; iz < f->getNZ(); iz++) {
584  Cell *c = f->getCell(ix, iy, iz);
585  c->getPrimVar(eos, tau, e, p, nb, nq, ns, vx, vy, vz);
586  if (e <= 0.) continue;
587  // NSquant(ix, iy, iz, piNS, PiNS, dmu, du) ;
588  //############## set NS values assuming initial zero flow + Bjorken z
589  //flow
590  double T, mub, muq, mus;
591  double etaS, zetaS;
592  double s =
593  eos->s(e, nb, nq, ns); // entropy density in the current cell
594  eos->eos(e, nb, nq, ns, T, mub, muq, mus, p);
595  trcoeff->getEta(e, T, etaS, zetaS);
596  for (int i = 0; i < 4; i++)
597  for (int j = 0; j < 4; j++) piNS[i][j] = 0.0; // reset piNS
598  piNS[1][1] = piNS[2][2] = 2.0 / 3.0 * etaS * s / tau / 5.068;
599  piNS[3][3] = -2.0 * piNS[1][1];
600  PiNS = 0.0;
601  for (int i = 0; i < 4; i++)
602  for (int j = 0; j <= i; j++) c->setpi(i, j, piNS[i][j]);
603  c->setPi(PiNS);
604  }
605  cout << "setNS done\n";
606 }
607 
608 void Hydro::ISformal() {
609  double e, p, nb, nq, ns, vx, vy, vz, T, mub, muq, mus;
610  double piNS[4][4], PiNS, dmu[4][4], du, pi[4][4], piH[4][4], Pi, PiH;
611  const double gmumu[4] = {1., -1., -1., -1.};
612 
613  // loop #1 (relaxation+source terms)
614  for (int ix = 0; ix < f->getNX(); ix++)
615  for (int iy = 0; iy < f->getNY(); iy++)
616  for (int iz = 0; iz < f->getNZ(); iz++) {
617  Cell *c = f->getCell(ix, iy, iz);
618  c->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vx, vy,
619  vz); // instead of getPrimVar()
620  if (e <= 0.) { // empty cell?
621  for (int i = 0; i < 4; i++)
622  for (int j = 0; j <= i; j++) {
623  c->setpiH0(i, j, 0.0);
624  c->setpi0(i, j, 0.0);
625  }
626  c->setPiH0(0.0);
627  c->setPi0(0.0);
628  } else { // non-empty cell
629  // 1) relaxation(pi)+source(pi) terms for half-step
630  double gamma = 1.0 / sqrt(1.0 - vx * vx - vy * vy - vz * vz);
631  double u[4];
632  u[0] = gamma;
633  u[1] = u[0] * vx;
634  u[2] = u[0] * vy;
635  u[3] = u[0] * vz;
636  // source term + tau*delta_Q_i/delta_tau
637  double flux[4];
638  for(int i=0; i<4; i++)
639  flux[i] = (tau-dt)*(c->getpi(0,i) + c->getPi()*u[0]*u[i]);
640  flux[0] += - (tau-dt)*c->getPi();
641  c->addFlux(flux[0], flux[1], flux[2], flux[3], 0., 0., 0.);
642  // now calculating viscous terms in NS limit
643  NSquant(ix, iy, iz, piNS, PiNS, dmu, du);
644  eos->eos(e, nb, nq, ns, T, mub, muq, mus, p);
645  //############# get relaxation times
646  double taupi, tauPi;
647  trcoeff->getTau(T, taupi, tauPi);
648  //#############
649  // relaxation term, piH,PiH-->half-step
650  for (int i = 0; i < 4; i++)
651  for (int j = 0; j <= i; j++) {
652 #ifdef FORMAL_SOLUTION
653  c->setpiH0(i, j, (c->getpi(i, j) - piNS[i][j]) *
654  exp(-dt / 2.0 / gamma / taupi) +
655  piNS[i][j]);
656 #else
657  c->setpiH0(i, j, c->getpi(i, j) - (c->getpi(i, j) - piNS[i][j]) *
658  dt / 2.0 / gamma / taupi);
659 #endif
660  }
661 #ifdef FORMAL_SOLUTION
662  c->setPiH0((c->getPi() - PiNS) * exp(-dt / 2.0 / gamma / tauPi) +
663  PiNS);
664 #else
665  c->setPiH0(c->getPi() -
666  (c->getPi() - PiNS) * dt / 2.0 / gamma / tauPi);
667 #endif
668  // sources from Christoffel symbols from \dot pi_munu
669  double tau1 = tau - dt * 0.75;
670  c->addpiH0(0, 0, -2. * vz * c->getpi(0, 3) / tau1 * dt /
671  2.); // *gamma/gamma
672  c->addpiH0(3, 3, -(2. * vz / tau1 * c->getpi(0, 3)) * dt / 2.);
673  c->addpiH0(
674  3, 0, -(vz / tau1 * c->getpi(0, 0) + vz / tau1 * c->getpi(3, 3)) *
675  dt / 2.);
676  c->addpiH0(1, 0, -vz / tau1 * c->getpi(1, 3) * dt / 2.);
677  c->addpiH0(2, 0, -vz / tau1 * c->getpi(2, 3) * dt / 2.);
678  c->addpiH0(3, 1, -(vz / tau1 * c->getpi(0, 1)) * dt / 2.);
679  c->addpiH0(3, 2, -(vz / tau1 * c->getpi(0, 2)) * dt / 2.);
680  // source from full IS equations (see draft for the description)
681  for (int i = 0; i < 4; i++)
682  for (int j = 0; j <= i; j++) {
683  c->addpiH0(i, j,
684  -4. / 3. * c->getpi(i, j) * du / gamma * dt / 2.);
685  for (int k = 0; k < 4;
686  k++) // terms to achieve better transverality to u^mu
687  for (int l = 0; l < 4; l++)
688  c->addpiH0(i, j,
689  -c->getpi(i, k) * u[j] * u[l] * dmu[l][k] *
690  gmumu[k] / gamma * dt / 2. -
691  c->getpi(j, k) * u[i] * u[l] * dmu[l][k] *
692  gmumu[k] / gamma * dt / 2.);
693  }
694  c->addPiH0(-4. / 3. * c->getPi() * du / gamma * dt / 2.);
695  // 1) relaxation(piH)+source(piH) terms for full-step
696  for (int i = 0; i < 4; i++)
697  for (int j = 0; j <= i; j++) {
698 #ifdef FORMAL_SOLUTION
699  c->setpi0(i, j, (c->getpi(i, j) - piNS[i][j]) *
700  exp(-dt / gamma / taupi) +
701  piNS[i][j]);
702 #else
703  c->setpi0(i, j, c->getpi(i, j) - (c->getpiH0(i, j) - piNS[i][j]) *
704  dt / gamma / taupi);
705 #endif
706  }
707 #ifdef FORMAL_SOLUTION
708  c->setPi0((c->getPi() - PiNS) * exp(-dt / gamma / tauPi) + PiNS);
709 #else
710  c->setPi0(c->getPi() - (c->getPiH0() - PiNS) * dt / gamma / tauPi);
711 #endif
712  tau1 = tau - dt * 0.5;
713  c->addpi0(0, 0,
714  -2. * vz / tau1 * c->getpiH0(0, 3) * dt); // *gamma/gamma
715  c->addpi0(3, 3, -(2. * vz / tau1 * c->getpiH0(0, 3)) * dt);
716  c->addpi0(3, 0, -(vz / tau1 * c->getpiH0(0, 0) +
717  vz / tau1 * c->getpiH0(3, 3)) *
718  dt);
719  c->addpi0(1, 0, -vz / tau1 * c->getpiH0(1, 3) * dt);
720  c->addpi0(2, 0, -vz / tau1 * c->getpiH0(2, 3) * dt);
721  c->addpi0(3, 1, -(vz / tau1 * c->getpiH0(0, 1)) * dt);
722  c->addpi0(3, 2, -(vz / tau1 * c->getpiH0(0, 2)) * dt);
723  // source from full IS equations (see draft for the description)
724  for (int i = 0; i < 4; i++)
725  for (int j = 0; j <= i; j++) {
726  c->addpi0(i, j, -4. / 3. * c->getpiH0(i, j) * du / gamma * dt);
727  for (int k = 0; k < 4;
728  k++) // terms to achieve better transverality to u^mu
729  for (int l = 0; l < 4; l++)
730  c->addpi0(i, j, -c->getpiH0(i, k) * u[j] * u[l] * dmu[l][k] *
731  gmumu[k] / gamma * dt -
732  c->getpiH0(j, k) * u[i] * u[l] *
733  dmu[l][k] * gmumu[k] / gamma * dt);
734  }
735  c->addPi0(-4. / 3. * c->getPiH0() * du / gamma * dt);
736  } // end non-empty cell
737  } // end loop #1
738 
739  // 3) -- advection ---
740  for (int ix = 0; ix < f->getNX(); ix++)
741  for (int iy = 0; iy < f->getNY(); iy++)
742  for (int iz = 0; iz < f->getNZ(); iz++) {
743  Cell *c = f->getCell(ix, iy, iz);
744  c->getPrimVarHCenter(eos, tau - 0.5 * dt, e, p, nb, nq, ns, vx, vy,
745  vz); // getPrimVar() before
746  if (e <= 0.) continue;
747  double xm = -vx * dt / f->getDx();
748  double ym = -vy * dt / f->getDy();
749  double zm = -vz * dt / f->getDz() / (tau - 0.5 * dt);
750  double xmH = -vx * dt / f->getDx() / 2.0;
751  double ymH = -vy * dt / f->getDy() / 2.0;
752  double zmH = -vz * dt / f->getDz() / (tau - 0.5 * dt) / 2.0;
753  double wx[2] = {(1. - fabs(xm)), fabs(xm)};
754  double wy[2] = {(1. - fabs(ym)), fabs(ym)};
755  double wz[2] = {(1. - fabs(zm)), fabs(zm)};
756  double wxH[2] = {(1. - fabs(xmH)), fabs(xmH)};
757  double wyH[2] = {(1. - fabs(ymH)), fabs(ymH)};
758  double wzH[2] = {(1. - fabs(zmH)), fabs(zmH)};
759  for (int i = 0; i < 4; i++)
760  for (int j = 0; j < 4; j++) {
761  pi[i][j] = piH[i][j] = 0.0;
762  }
763  Pi = PiH = 0.0;
764  for (int jx = 0; jx < 2; jx++)
765  for (int jy = 0; jy < 2; jy++)
766  for (int jz = 0; jz < 2; jz++) {
767  // pi,Pi-->full step, piH,PiH-->half-step
768  Cell *c1 = f->getCell(ix + jx * sign(xm), iy + jy * sign(ym),
769  iz + jz * sign(zm));
770  for (int i = 0; i < 4; i++)
771  for (int j = 0; j < 4; j++) {
772  pi[i][j] += wx[jx] * wy[jy] * wz[jz] * c1->getpi0(i, j);
773  piH[i][j] += wxH[jx] * wyH[jy] * wzH[jz] * c1->getpiH0(i, j);
774  }
775  Pi += wx[jx] * wy[jy] * wz[jz] * c1->getPi0();
776  PiH += wxH[jx] * wyH[jy] * wzH[jz] * c1->getPiH0();
777  }
778 
779  //--------- debug part: trace check, diag check, transversality check
780  for (int i = 0; i < 4; i++)
781  for (int j = 0; j < 4; j++) {
782  if (pi[i][j] != 0. && fabs(1.0 - pi[j][i] / pi[i][j]) > 1e-10)
783  cout << "non-diag: " << pi[i][j] << " " << pi[j][i] << endl;
784  }
785  //------ end debug
786  //======= hydro applicability check (viscous corrections limiter):
787  // double maxT0 = max(fabs((e+p)*vx*vx/(1.-vx*vx-vy*vy-vz*vz)+p),
788  // fabs((e+p)*vy*vy/(1.-vx*vx-vy*vy-vz*vz)+p)) ;
789  double maxT0 = max((e + p) / (1. - vx * vx - vy * vy - vz * vz) - p,
790  (e + p) * (vx * vx + vy * vy + vz * vz) /
791  (1. - vx * vx - vy * vy - vz * vz) +
792  p);
793  // double maxpi = max(fabs(pi[1][1]),fabs(pi[2][2])) ;
794  double maxpi = 0.;
795  for (int i = 0; i < 4; i++)
796  for (int j = 0; j < 4; j++)
797  if (fabs(pi[i][j]) > maxpi) maxpi = fabs(pi[i][j]);
798  bool rescaled = false;
799  if (maxT0 / maxpi < 1.0) {
800  for (int i = 0; i < 4; i++)
801  for (int j = 0; j < 4; j++) {
802  pi[i][j] = 0.1 * pi[i][j] * maxT0 / maxpi;
803  piH[i][j] = 0.1 * piH[i][j] * maxT0 / maxpi;
804  }
805  rescaled = true;
806  }
807  if (fabs(Pi) > p) {
808  if (Pi != 0.) Pi = 0.1 * Pi / fabs(Pi) * p;
809  if (PiH != 0.) PiH = 0.1 * PiH / fabs(PiH) * p;
810  rescaled = true;
811  }
812  if (rescaled)
813  c->setViscCorrCutFlag(maxT0 / maxpi);
814  else
815  c->setViscCorrCutFlag(1.);
816  // updating to the new values
817  for (int i = 0; i < 4; i++)
818  for (int j = 0; j <= i; j++) {
819  c->setpi(i, j, pi[i][j]);
820  c->setpiH(i, j, piH[i][j]);
821  }
822  c->setPi(Pi);
823  c->setPiH(PiH);
824  // source term - (tau+dt)*delta_Q_(i+1)/delta_tau
825  double gamma = 1.0 / sqrt(1.0 - vx * vx - vy * vy - vz * vz);
826  double u[4];
827  u[0] = gamma;
828  u[1] = u[0] * vx;
829  u[2] = u[0] * vy;
830  u[3] = u[0] * vz;
831  double flux[4];
832  for(int i=0; i<4; i++)
833  flux[i] = - tau*(c->getpi(0,i) + c->getPi()*u[0]*u[i]);
834  flux[0] += tau*c->getPi();
835  c->addFlux(flux[0], flux[1], flux[2], flux[3], 0., 0., 0.);
836  } // advection loop (all cells)
837 }
838 
839 
840 // this procedure explicitly uses T_==0, X_==1, Y_==2, Z_==3
841 void Hydro::visc_flux(Cell *left, Cell *right, int direction) {
842  double flux[4];
843  int ind2 = 0;
844  double dxa = 0.;
845  // exit if noth cells are not full with matter
846  if (left->getM(direction) < 1. && right->getM(direction) < 1.) return;
847 
848  if (direction == X_)
849  dxa = f->getDx();
850  else if (direction == Y_)
851  dxa = f->getDy();
852  else if (direction == Z_)
853  dxa = f->getDz() * (tau + 0.5 * dt);
854  double e, p, nb, nq, ns, vxl, vyl, vzl, vxr, vyr, vzr;
855  // we need to know the velocities at both cell centers at (n+1/2) in order to
856  // interpolate to
857  // get the value at the interface
858  left->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vxl, vyl, vzl);
859  right->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vxr, vyr, vzr);
860  vxl = 0.5 * (vxl + vxr);
861  vyl = 0.5 * (vyl + vyr);
862  vzl = 0.5 * (vzl + vzr);
863  double v = sqrt(vxl * vxl + vyl * vyl + vzl * vzl);
864  if (v > 1.) {
865  vxl = 0.99 * vxl / v;
866  vyl = 0.99 * vyl / v;
867  vzl = 0.99 * vzl / v;
868  }
869  double gamma = 1. / sqrt(1. - v * v);
870  double uuu[4] = {gamma, gamma * vxl, gamma * vyl, gamma * vzl};
871  double gmumu[4] = {1., -1., -1., -1.};
872  if (direction == X_)
873  ind2 = 1;
874  else if (direction == Y_)
875  ind2 = 2;
876  else if (direction == Z_)
877  ind2 = 3;
878  for (int ind1 = 0; ind1 < 4; ind1++) {
879  flux[ind1] = 0.5 * (left->getpiH(ind1, ind2) + right->getpiH(ind1, ind2));
880  if (ind1 == ind2)
881  flux[ind1] += -0.5 * (left->getPiH() + right->getPiH()) *
882  gmumu[ind1]; // gmunu is diagonal
883  flux[ind1] +=
884  0.5 * (left->getPiH() + right->getPiH()) * uuu[ind1] * uuu[ind2];
885  }
886  for (int i = 0; i < 4; i++) flux[i] = flux[i] * (tau - 0.5*dt) * dt / dxa;
887  left->addFlux(-flux[T_], -flux[X_], -flux[Y_], -flux[Z_], 0., 0., 0.);
888  right->addFlux(flux[T_], flux[X_], flux[Y_], flux[Z_], 0., 0., 0.);
889 }
890 
891 void Hydro::performStep(void) {
892  debugRiemann = false; // turn off debug output
893 
894  f->updateM(tau, dt);
895 
896  tau_z = dt / 2. / log(1 + dt / 2. / tau);
897 
898  //-----PREDICTOR-ideal
899  for (int iy = 0; iy < f->getNY(); iy++)
900  for (int iz = 0; iz < f->getNZ(); iz++)
901  for (int ix = 0; ix < f->getNX(); ix++) {
902  Cell *c = f->getCell(ix, iy, iz);
903  c->saveQprev();
904  c->clearFlux();
905  }
906  // X dir
907  for (int iy = 0; iy < f->getNY(); iy++)
908  for (int iz = 0; iz < f->getNZ(); iz++)
909  for (int ix = 0; ix < f->getNX() - 1; ix++) {
910  hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix + 1, iy, iz), X_,
911  PREDICT);
912  }
913  // cout << "predictor X done\n" ;
914  // Y dir
915  for (int iz = 0; iz < f->getNZ(); iz++)
916  for (int ix = 0; ix < f->getNX(); ix++)
917  for (int iy = 0; iy < f->getNY() - 1; iy++) {
918  hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy + 1, iz), Y_,
919  PREDICT);
920  }
921  // cout << "predictor Y done\n" ;
922  // Z dir
923  for (int ix = 0; ix < f->getNX(); ix++)
924  for (int iy = 0; iy < f->getNY(); iy++)
925  for (int iz = 0; iz < f->getNZ() - 1; iz++) {
926  hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy, iz + 1), Z_,
927  PREDICT);
928  }
929  // cout << "predictor Z done\n" ;
930 
931  for (int iy = 0; iy < f->getNY(); iy++)
932  for (int iz = 0; iz < f->getNZ(); iz++)
933  for (int ix = 0; ix < f->getNX(); ix++) {
934  Cell *c = f->getCell(ix, iy, iz);
935  source_step(ix, iy, iz, PREDICT);
936  c->updateQtoQhByFlux();
937  c->clearFlux();
938  }
939 
940  //----CORRECTOR-ideal
941 
942  tau_z = dt / log(1 + dt / tau);
943  // X dir
944  for (int iy = 0; iy < f->getNY(); iy++)
945  for (int iz = 0; iz < f->getNZ(); iz++)
946  for (int ix = 0; ix < f->getNX() - 1; ix++) {
947  hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix + 1, iy, iz), X_,
948  CORRECT);
949  }
950  // cout << "corrector X done\n" ;
951  // Y dir
952  for (int iz = 0; iz < f->getNZ(); iz++)
953  for (int ix = 0; ix < f->getNX(); ix++)
954  for (int iy = 0; iy < f->getNY() - 1; iy++) {
955  hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy + 1, iz), Y_,
956  CORRECT);
957  }
958  // cout << "corrector Y done\n" ;
959  // Z dir
960  for (int ix = 0; ix < f->getNX(); ix++)
961  for (int iy = 0; iy < f->getNY(); iy++)
962  for (int iz = 0; iz < f->getNZ() - 1; iz++) {
963  hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy, iz + 1), Z_,
964  CORRECT);
965  }
966  // cout << "corrector Z done\n" ;
967 
968  for (int iy = 0; iy < f->getNY(); iy++)
969  for (int iz = 0; iz < f->getNZ(); iz++)
970  for (int ix = 0; ix < f->getNX(); ix++) {
971  Cell *c = f->getCell(ix, iy, iz);
972  source_step(ix, iy, iz, CORRECT);
973  c->updateByFlux();
974  c->clearFlux();
975  }
976  tau += dt;
977  f->correctImagCells();
978 
979  //===== viscous part ======
980  if (trcoeff->isViscous()) {
981  ISformal(); // evolution of viscous quantities according to IS equations
982 
983  // X dir
984  for (int iy = 0; iy < f->getNY(); iy++)
985  for (int iz = 0; iz < f->getNZ(); iz++)
986  for (int ix = 0; ix < f->getNX() - 1; ix++) {
987  visc_flux(f->getCell(ix, iy, iz), f->getCell(ix + 1, iy, iz), X_);
988  }
989  // cout << "visc_flux X done\n" ;
990  // Y dir
991  for (int iz = 0; iz < f->getNZ(); iz++)
992  for (int ix = 0; ix < f->getNX(); ix++)
993  for (int iy = 0; iy < f->getNY() - 1; iy++) {
994  visc_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy + 1, iz), Y_);
995  }
996  // cout << "visc_flux Y done\n" ;
997  // Z dir
998  for (int ix = 0; ix < f->getNX(); ix++)
999  for (int iy = 0; iy < f->getNY(); iy++)
1000  for (int iz = 0; iz < f->getNZ() - 1; iz++) {
1001  visc_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy, iz + 1), Z_);
1002  }
1003 
1004  for (int iy = 0; iy < f->getNY(); iy++)
1005  for (int iz = 0; iz < f->getNZ(); iz++)
1006  for (int ix = 0; ix < f->getNX(); ix++) {
1007  visc_source_step(ix, iy, iz);
1008  f->getCell(ix, iy, iz)->updateByFlux();
1009  f->getCell(ix, iy, iz)->clearFlux();
1010  }
1011  } else { // end viscous part
1012  }
1013  //==== finishing work ====
1014  f->correctImagCellsFull();
1015 }