Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:19:16 

0001 /******************************************************************************
0002 *                                                                             *
0003 *            vHLLE : a 3D viscous hydrodynamic code                           *
0004 *            version 1.0,            November 2013                            *
0005 *            by Iurii Karpenko                                                *
0006 *  contact:  yu.karpenko@gmail.com                                            *
0007 *  For the detailed description please refer to:                              *
0008 *  http://arxiv.org/abs/1312.4160                                             *
0009 *                                                                             *
0010 *  This code can be freely used and redistributed, provided that this         *
0011 *  copyright appear in all the copies. If you decide to make modifications    *
0012 *  to the code, please contact the authors, especially if you plan to publish *
0013 * the results obtained with such modified code. Any publication of results    *
0014 * obtained using this code must include the reference to                      *
0015 * arXiv:1312.4160 [nucl-th] or the published version of it, when available.   *
0016 *                                                                             *
0017 *******************************************************************************/
0018 
0019 #include <iostream>
0020 #include <fstream>
0021 #include <iomanip>
0022 #include <algorithm>
0023 #include <unistd.h>
0024 #include "hdo.h"
0025 #include "inc.h"
0026 #include "rmn.h"
0027 #include "fld.h"
0028 #include "eos.h"
0029 #include "cll.h"
0030 #include "trancoeff.h"
0031 
0032 using namespace std;
0033 
0034 extern bool debugRiemann;
0035 
0036 double sign(double x) {
0037   if (x > 0)
0038     return 1.;
0039   else if (x < 0.)
0040     return -1.;
0041   else
0042     return 0.;
0043 }
0044 
0045 // this version contains NO PRE-ADVECTION for the IS solution
0046 
0047 // enable this to use formal solution for the relaxation part of
0048 // Israel-Stewart equations (not recommended)
0049 //#define FORMAL_SOLUTION
0050 // else: use 1st order finite difference update
0051 
0052 Hydro::Hydro(Fluid *_f, EoS *_eos, TransportCoeff *_trcoeff, double _t0,
0053              double _dt) {
0054   eos = _eos;
0055   trcoeff = _trcoeff;
0056   f = _f;
0057   dt = _dt;
0058   tau = _t0;
0059 }
0060 
0061 Hydro::~Hydro() {}
0062 
0063 void Hydro::setDtau(double deltaTau) {
0064   dt = deltaTau;
0065   if (dt > f->getDx() / 2. ||
0066       dt > f->getDy() / 2. /*|| dt>tau*f->getDz()/2. */) {
0067     cout << "too big delta_tau " << dt << "  " << f->getDx() << "  "
0068          << f->getDy() << "  " << tau * f->getDz() << endl;
0069     exit(1);
0070   }
0071 }
0072 
0073 void Hydro::hlle_flux(Cell *left, Cell *right, int direction, int mode) {
0074   // for all variables, suffix "l" = left state, "r" = right state
0075   // with respect to the cell boundary
0076   double el, er, pl, pr, nbl, nql, nsl, nbr, nqr, nsr, vxl, vxr, vyl, vyr, vzl,
0077       vzr, bl = 0., br = 0., csb, vb, El, Er, dx = 0.;
0078   double Ftl = 0., Fxl = 0., Fyl = 0., Fzl = 0., Fbl = 0., Fql = 0., Fsl = 0.,
0079          Ftr = 0., Fxr = 0., Fyr = 0., Fzr = 0., Fbr = 0., Fqr = 0., Fsr = 0.;
0080   double U1l, U2l, U3l, U4l, Ubl, Uql, Usl, U1r, U2r, U3r, U4r, Ubr, Uqr, Usr;
0081   double flux[7];
0082   const double dta = mode == 0 ? dt / 2. : dt;
0083   double tauFactor;  // fluxes are also multiplied by tau
0084   if (mode == PREDICT) {
0085     // get primitive quantities from Q_{i+} at previous timestep
0086     left->getPrimVarRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
0087                           direction);
0088     // ... and Q_{(i+1)-}
0089     right->getPrimVarLeft(eos, tau, er, pr, nbr, nqr, nsr, vxr, vyr, vzr,
0090                           direction);
0091     El = (el + pl) / (1 - vxl * vxl - vyl * vyl - vzl * vzl);
0092     Er = (er + pr) / (1 - vxr * vxr - vyr * vyr - vzr * vzr);
0093     tauFactor = tau;
0094   } else {
0095     // use half-step updated Q's for corrector step
0096     left->getPrimVarHRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
0097                            direction);
0098     right->getPrimVarHLeft(eos, tau, er, pr, nbr, nqr, nsr, vxr, vyr, vzr,
0099                            direction);
0100     El = (el + pl) / (1 - vxl * vxl - vyl * vyl - vzl * vzl);
0101     Er = (er + pr) / (1 - vxr * vxr - vyr * vyr - vzr * vzr);
0102     tauFactor = tau + dt / 2.;
0103   }
0104 
0105   if (el < 0.) {
0106     el = 0.;
0107     pl = 0.;
0108   }
0109   if (er < 0.) {
0110     er = 0.;
0111     pr = 0.;
0112   }
0113 
0114   if (el > 1e10) {
0115     cout << "e>1e10; debug info below:\n";
0116     left->Dump(tau);
0117     debugRiemann = true;
0118     if (mode == PREDICT)
0119       left->getPrimVarRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
0120                             direction);
0121     else
0122       left->getPrimVarHRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
0123                              direction);
0124     debugRiemann = false;
0125     exit(0);
0126   }
0127 
0128   // skip the procedure for two empty cells
0129   if (el == 0. && er == 0.) return;
0130   if (pr < 0.) {
0131     cout << "Negative pressure" << endl;
0132     left->getPrimVarRight(eos, tau, el, pl, nbl, nql, nsl, vxl, vyl, vzl,
0133                           direction);
0134     right->getPrimVarLeft(eos, tau, er, pr, nbr, nqr, nsr, vxr, vyr, vzr,
0135                           direction);
0136   }
0137 
0138   // skip the procedure for two partially vacuum cells
0139   if (left->getM(direction) < 1. && right->getM(direction) < 1.) return;
0140 
0141   double gammal = 1. / sqrt(1 - vxl * vxl - vyl * vyl - vzl * vzl);
0142   double gammar = 1. / sqrt(1 - vxr * vxr - vyr * vyr - vzr * vzr);
0143   U1l = gammal * gammal * (el + pl) * vxl;
0144   U2l = gammal * gammal * (el + pl) * vyl;
0145   U3l = gammal * gammal * (el + pl) * vzl;
0146   U4l = gammal * gammal * (el + pl) - pl;
0147   Ubl = gammal * nbl;
0148   Uql = gammal * nql;
0149   Usl = gammal * nsl;
0150 
0151   U1r = gammar * gammar * (er + pr) * vxr;
0152   U2r = gammar * gammar * (er + pr) * vyr;
0153   U3r = gammar * gammar * (er + pr) * vzr;
0154   U4r = gammar * gammar * (er + pr) - pr;
0155   Ubr = gammar * nbr;
0156   Uqr = gammar * nqr;
0157   Usr = gammar * nsr;
0158 
0159   if (direction == X_) {
0160     Ftl = U4l * vxl + pl * vxl;
0161     Fxl = U1l * vxl + pl;
0162     Fyl = U2l * vxl;
0163     Fzl = U3l * vxl;
0164     Fbl = Ubl * vxl;
0165     Fql = Uql * vxl;
0166     Fsl = Usl * vxl;
0167 
0168     Ftr = U4r * vxr + pr * vxr;
0169     Fxr = U1r * vxr + pr;
0170     Fyr = U2r * vxr;
0171     Fzr = U3r * vxr;
0172     Fbr = Ubr * vxr;
0173     Fqr = Uqr * vxr;
0174     Fsr = Usr * vxr;
0175 
0176     // for the case of constant c_s only
0177     csb = sqrt(eos->cs2() + 0.5 * sqrt(El * Er) / pow(sqrt(El) + sqrt(Er), 2) *
0178                                 pow(vxl - vxr, 2));
0179     vb = (sqrt(El) * vxl + sqrt(Er) * vxr) / (sqrt(El) + sqrt(Er));
0180     bl = min(0., min((vb - csb) / (1 - vb * csb),
0181                      (vxl - eos->cs()) / (1 - vxl * eos->cs())));
0182     br = max(0., max((vb + csb) / (1 + vb * csb),
0183                      (vxr + eos->cs()) / (1 + vxr * eos->cs())));
0184 
0185     dx = f->getDx();
0186 
0187     // bl or br in the case of boundary with vacuum
0188     if (el == 0.) bl = -1.;
0189     if (er == 0.) br = 1.;
0190   }
0191   if (direction == Y_) {
0192     Ftl = U4l * vyl + pl * vyl;
0193     Fxl = U1l * vyl;
0194     Fyl = U2l * vyl + pl;
0195     Fzl = U3l * vyl;
0196     Fbl = Ubl * vyl;
0197     Fql = Uql * vyl;
0198     Fsl = Usl * vyl;
0199 
0200     Ftr = U4r * vyr + pr * vyr;
0201     Fxr = U1r * vyr;
0202     Fyr = U2r * vyr + pr;
0203     Fzr = U3r * vyr;
0204     Fbr = Ubr * vyr;
0205     Fqr = Uqr * vyr;
0206     Fsr = Usr * vyr;
0207 
0208     // for the case of constant c_s only
0209     csb = sqrt(eos->cs2() + 0.5 * sqrt(El * Er) / pow(sqrt(El) + sqrt(Er), 2) *
0210                                 pow(vyl - vyr, 2));
0211     vb = (sqrt(El) * vyl + sqrt(Er) * vyr) / (sqrt(El) + sqrt(Er));
0212     bl = min(0., min((vb - csb) / (1 - vb * csb),
0213                      (vyl - eos->cs()) / (1 - vyl * eos->cs())));
0214     br = max(0., max((vb + csb) / (1 + vb * csb),
0215                      (vyr + eos->cs()) / (1 + vyr * eos->cs())));
0216 
0217     dx = f->getDy();
0218 
0219     // bl or br in the case of boundary with vacuum
0220     if (el == 0.) bl = -1.;
0221     if (er == 0.) br = 1.;
0222   }
0223   if (direction == Z_) {
0224     double tau1 = tau_z;
0225     Ftl = U4l * vzl / tau1 + pl * vzl / tau1;
0226     Fxl = U1l * vzl / tau1;
0227     Fyl = U2l * vzl / tau1;
0228     Fzl = U3l * vzl / tau1 + pl / tau1;
0229     Fbl = Ubl * vzl / tau1;
0230     Fql = Uql * vzl / tau1;
0231     Fsl = Usl * vzl / tau1;
0232 
0233     Ftr = U4r * vzr / tau1 + pr * vzr / tau1;
0234     Fxr = U1r * vzr / tau1;
0235     Fyr = U2r * vzr / tau1;
0236     Fzr = U3r * vzr / tau1 + pr / tau1;
0237     Fbr = Ubr * vzr / tau1;
0238     Fqr = Uqr * vzr / tau1;
0239     Fsr = Usr * vzr / tau1;
0240 
0241     // for the case of constant c_s only
0242     // factor 1/tau accounts for eta-coordinate
0243 
0244     // different estimate
0245     csb = sqrt(eos->cs2() + 0.5 * sqrt(El * Er) / pow(sqrt(El) + sqrt(Er), 2) *
0246                                 pow(vzl - vzr, 2));
0247     vb = (sqrt(El) * vzl + sqrt(Er) * vzr) / (sqrt(El) + sqrt(Er));
0248     bl = 1. / tau * min(0., min((vb - csb) / (1 - vb * csb),
0249                                 (vzl - eos->cs()) / (1 - vzl * eos->cs())));
0250     br = 1. / tau * max(0., max((vb + csb) / (1 + vb * csb),
0251                                 (vzr + eos->cs()) / (1 + vzr * eos->cs())));
0252 
0253     dx = f->getDz();
0254 
0255     // bl or br in the case of boundary with vacuum
0256     if (el == 0.) bl = -1. / tau;
0257     if (er == 0.) br = 1. / tau;
0258   }
0259 
0260   if (bl == 0. && br == 0.) return;
0261 
0262   // finally, HLLE formula for the fluxes
0263   flux[T_] = tauFactor * dta / dx *
0264              (-bl * br * (U4l - U4r) + br * Ftl - bl * Ftr) / (-bl + br);
0265   flux[X_] = tauFactor * dta / dx *
0266              (-bl * br * (U1l - U1r) + br * Fxl - bl * Fxr) / (-bl + br);
0267   flux[Y_] = tauFactor * dta / dx *
0268              (-bl * br * (U2l - U2r) + br * Fyl - bl * Fyr) / (-bl + br);
0269   flux[Z_] = tauFactor * dta / dx *
0270              (-bl * br * (U3l - U3r) + br * Fzl - bl * Fzr) / (-bl + br);
0271   flux[NB_] = tauFactor * dta / dx *
0272               (-bl * br * (Ubl - Ubr) + br * Fbl - bl * Fbr) / (-bl + br);
0273   flux[NQ_] = tauFactor * dta / dx *
0274               (-bl * br * (Uql - Uqr) + br * Fql - bl * Fqr) / (-bl + br);
0275   flux[NS_] = tauFactor * dta / dx *
0276               (-bl * br * (Usl - Usr) + br * Fsl - bl * Fsr) / (-bl + br);
0277 
0278   if (flux[NB_] != flux[NB_]) {  // if things failed
0279     cout << "---- error in hlle_flux: f_nb undefined!\n";
0280     cout << setw(12) << U4l << setw(12) << U1l << setw(12) << U2l << setw(12)
0281          << U3l << endl;
0282     cout << setw(12) << U4r << setw(12) << U1r << setw(12) << U2r << setw(12)
0283          << U3r << endl;
0284     cout << setw(12) << Ubl << setw(12) << Uql << setw(12) << Usl << endl;
0285     cout << setw(12) << Ubr << setw(12) << Uqr << setw(12) << Usr << endl;
0286     cout << setw(12) << Ftl << setw(12) << Fxl << setw(12) << Fyl << setw(12)
0287          << Fzl << endl;
0288     cout << setw(12) << Ftr << setw(12) << Fxr << setw(12) << Fyr << setw(12)
0289          << Fzr << endl;
0290     exit(1);
0291   }
0292 
0293   // update the cumulative fluxes in both neighbouring cells
0294   left->addFlux(-flux[T_], -flux[X_], -flux[Y_], -flux[Z_], -flux[NB_],
0295                 -flux[NQ_], -flux[NS_]);
0296   right->addFlux(flux[T_], flux[X_], flux[Y_], flux[Z_], flux[NB_], flux[NQ_],
0297                  flux[NS_]);
0298 }
0299 
0300 void Hydro::source(double tau1, double x, double y, double z, double e,
0301                    double p, double nb, double nq, double ns, double vx,
0302                    double vy, double vz, double S[7]) {
0303   double Q[7];
0304   transformCV(e, p, nb, nq, ns, vx, vy, vz, Q);
0305   S[T_] = -Q[T_] * vz * vz - p * (1. + vz * vz);
0306   S[X_] = 0.;
0307   S[Y_] = 0.;
0308   S[Z_] = -Q[Z_];
0309   S[NB_] = 0.;
0310   S[NQ_] = 0.;
0311   S[NS_] = 0.;
0312 }
0313 
0314 void Hydro::source_step(int ix, int iy, int iz, int mode) {
0315   double _dt;
0316   if (mode == PREDICT)
0317     _dt = dt / 2.;
0318   else
0319     _dt = dt;
0320 
0321   double e, p, nb, nq, ns, vx, vy, vz;
0322   double k[7];
0323 
0324   double x = f->getX(ix), y = f->getY(iy), z = f->getZ(iz);
0325   Cell *c = f->getCell(ix, iy, iz);
0326 
0327   if (mode == PREDICT) {
0328     c->getPrimVar(eos, tau, e, p, nb, nq, ns, vx, vy, vz);
0329   } else {
0330     c->getPrimVarHCenter(eos, tau + dt / 2., e, p, nb, nq, ns, vx, vy, vz);
0331   }
0332 
0333   source(tau, x, y, z, e, p, nb, nq, ns, vx, vy, vz, k);  // setting k1
0334   for (int i = 0; i < 7; i++) k[i] *= _dt;
0335 
0336   if (k[NB_] != k[NB_]) {  // something failed
0337     cout << "---- error in source_step: k_nb undefined!\n";
0338     cout << setw(12) << k[0] << setw(12) << k[1] << setw(12) << k[2] << setw(12)
0339          << k[3] << endl;
0340     cout << setw(12) << k[4] << setw(12) << k[5] << setw(12) << k[6] << endl;
0341     exit(1);
0342   }
0343   c->addFlux(k[T_], k[X_], k[Y_], k[Z_], k[NB_], k[NQ_], k[NS_]);
0344 }
0345 
0346 void Hydro::visc_source_step(int ix, int iy, int iz) {
0347   double e, p, nb, nq, ns, vx, vy, vz;
0348   double uuu[4];
0349   double k[7];
0350 
0351   Cell *c = f->getCell(ix, iy, iz);
0352 
0353   c->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vx, vy, vz);
0354   if (e <= 0.) return;
0355   uuu[0] = 1. / sqrt(1. - vx * vx - vy * vy - vz * vz);
0356   uuu[1] = uuu[0] * vx;
0357   uuu[2] = uuu[0] * vy;
0358   uuu[3] = uuu[0] * vz;
0359 
0360   k[T_] = - c->getpiH(3, 3) +
0361             c->getPiH() * ( - 1.0 - uuu[3] * uuu[3]);
0362   k[X_] = 0.;
0363   k[Y_] = 0.;
0364   k[Z_] = - (c->getpiH(0, 3) + c->getPiH() * uuu[0] * uuu[3]);
0365   for (int i = 0; i < 4; i++) k[i] *= dt;
0366 
0367   c->addFlux(k[T_], k[X_], k[Y_], k[Z_], 0., 0., 0.);
0368 }
0369 
0370 // for the procedure below, the following approximations are used:
0371 // dv/d(tau) = v^{t+dt}_ideal - v^{t}
0372 // dv/dx_i ~ v^{x+dx}-v{x-dx},
0373 // which makes sense after non-viscous step
0374 void Hydro::NSquant(int ix, int iy, int iz, double pi[4][4], double &Pi,
0375                     double dmu[4][4], double &du) {
0376   const double VMIN = 1e-2;
0377   const double UDIFF = 3.0;
0378   double e0, e1, p, nb, nq, ns, vx1, vy1, vz1, vx0, vy0, vz0, vxH, vyH, vzH;
0379   double ut0, ux0, uy0, uz0, ut1, ux1, uy1, uz1;
0380   //    double dmu [4][4] ; // \partial_\mu u^\nu matrix
0381   // coordinates: 0=tau, 1=x, 2=y, 3=eta
0382   double Z[4][4][4][4];  // Z[mu][nu][lambda][rho]
0383   double uuu[4];         // the 4-velocity
0384   double gmunu[4][4] = {{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0},
0385                         {0, 0, 0, -1}};  // omit 1/tau^2 in g^{eta,eta}
0386   Cell *c = f->getCell(ix, iy, iz);
0387   double dx = f->getDx(), dy = f->getDy(), dz = f->getDz();
0388   // check if the cell is next to vacuum from +-x, +-y side:
0389   if (c->getNext(X_)->getMaxM() <= 0.9 || c->getNext(Y_)->getMaxM() <= 0.9 ||
0390       c->getPrev(X_)->getMaxM() <= 0.9 || c->getPrev(Y_)->getMaxM() <= 0.9 ||
0391       f->getCell(ix + 1, iy + 1, iz)->getMaxM() <= 0.9 ||
0392       f->getCell(ix + 1, iy - 1, iz)->getMaxM() <= 0.9 ||
0393       f->getCell(ix - 1, iy + 1, iz)->getMaxM() <= 0.9 ||
0394       f->getCell(ix - 1, iy - 1, iz)->getMaxM() <= 0.9) {
0395     for (int i = 0; i < 4; i++)
0396       for (int j = 0; j < 4; j++) {
0397         pi[i][j] = 0.;
0398         dmu[i][j] = 0.;
0399       }
0400     Pi = du = 0.;
0401     return;
0402   }
0403   // calculation of \partial_\mu u^\nu matrix
0404   // mu=first index, nu=second index
0405   // centered differences with respect to the values at (it+1/2, ix, iy, iz)
0406   // d_tau u^\mu
0407   c->getPrimVarPrev(eos, tau - dt, e0, p, nb, nq, ns, vx0, vy0, vz0);
0408   c->getPrimVar(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
0409   c->getPrimVarHCenter(eos, tau - 0.5 * dt, e1, p, nb, nq, ns, vxH, vyH, vzH);
0410   //############## get transport coefficients
0411   double T, mub, muq, mus;
0412   double etaS, zetaS;
0413   double s = eos->s(e1, nb, nq, ns);  // entropy density in the current cell
0414   eos->eos(e1, nb, nq, ns, T, mub, muq, mus, p);
0415   trcoeff->getEta(e1, T, etaS, zetaS);
0416   //##############
0417   // if(e1<0.00004) s=0. ; // negative pressure due to pi^zz for small e
0418   ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
0419   ux0 = ut0 * vx0;
0420   uy0 = ut0 * vy0;
0421   uz0 = ut0 * vz0;
0422   ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
0423   ux1 = ut1 * vx1;
0424   uy1 = ut1 * vy1;
0425   uz1 = ut1 * vz1;
0426   uuu[0] = 1.0 / sqrt(1.0 - vxH * vxH - vyH * vyH - vzH * vzH);
0427   uuu[1] = uuu[0] * vxH;
0428   uuu[2] = uuu[0] * vyH;
0429   uuu[3] = uuu[0] * vzH;
0430 
0431   dmu[0][0] = (ut1 * ut1 - ut0 * ut0) / 2. / uuu[0] / dt;
0432   dmu[0][1] = (ux1 * ux1 - ux0 * ux0) / 2. / uuu[1] / dt;
0433   dmu[0][2] = (uy1 * uy1 - uy0 * uy0) / 2. / uuu[2] / dt;
0434   dmu[0][3] = (uz1 * uz1 - uz0 * uz0) / 2. / uuu[3] / dt;
0435   if (fabs(0.5 * (ut1 + ut0) / ut1) > UDIFF) dmu[0][0] = (ut1 - ut0) / dt;
0436   if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / ux1) > UDIFF)
0437     dmu[0][1] = (ux1 - ux0) / dt;
0438   if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uy1) > UDIFF)
0439     dmu[0][2] = (uy1 - uy0) / dt;
0440   if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uz1) > UDIFF)
0441     dmu[0][3] = (uz1 - uz0) / dt;
0442   if (e1 <= 0. || e0 <= 0.) {  // matter-vacuum
0443     dmu[0][0] = dmu[0][1] = dmu[0][2] = dmu[0][3] = 0.;
0444   }
0445   // d_x u^\mu
0446   f->getCell(ix + 1, iy, iz)
0447       ->getPrimVarHCenter(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
0448   f->getCell(ix - 1, iy, iz)
0449       ->getPrimVarHCenter(eos, tau, e0, p, nb, nq, ns, vx0, vy0, vz0);
0450   if (e1 > 0. && e0 > 0.) {
0451     ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
0452     ux0 = ut0 * vx0;
0453     uy0 = ut0 * vy0;
0454     uz0 = ut0 * vz0;
0455     ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
0456     ux1 = ut1 * vx1;
0457     uy1 = ut1 * vy1;
0458     uz1 = ut1 * vz1;
0459     dmu[1][0] = 0.25 * (ut1 * ut1 - ut0 * ut0) / uuu[0] / dx;
0460     dmu[1][1] = 0.25 * (ux1 * ux1 - ux0 * ux0) / uuu[1] / dx;
0461     dmu[1][2] = 0.25 * (uy1 * uy1 - uy0 * uy0) / uuu[2] / dx;
0462     dmu[1][3] = 0.25 * (uz1 * uz1 - uz0 * uz0) / uuu[3] / dx;
0463     if (fabs(0.5 * (ut1 + ut0) / uuu[0]) > UDIFF)
0464       dmu[1][0] = 0.5 * (ut1 - ut0) / dx;
0465     if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / uuu[1]) > UDIFF)
0466       dmu[1][1] = 0.5 * (ux1 - ux0) / dx;
0467     if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uuu[2]) > UDIFF)
0468       dmu[1][2] = 0.5 * (uy1 - uy0) / dx;
0469     if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uuu[3]) > UDIFF)
0470       dmu[1][3] = 0.5 * (uz1 - uz0) / dx;
0471   } else {  // matter-vacuum
0472     dmu[1][0] = dmu[1][1] = dmu[1][2] = dmu[1][3] = 0.;
0473   }
0474   if (fabs(dmu[1][3]) > 1e+10)
0475     cout << "dmu[1][3]:  " << uz1 << "  " << uz0 << "  " << uuu[3] << endl;
0476   // d_y u^\mu
0477   f->getCell(ix, iy + 1, iz)
0478       ->getPrimVarHCenter(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
0479   f->getCell(ix, iy - 1, iz)
0480       ->getPrimVarHCenter(eos, tau, e0, p, nb, nq, ns, vx0, vy0, vz0);
0481   if (e1 > 0. && e0 > 0.) {
0482     ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
0483     ux0 = ut0 * vx0;
0484     uy0 = ut0 * vy0;
0485     uz0 = ut0 * vz0;
0486     ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
0487     ux1 = ut1 * vx1;
0488     uy1 = ut1 * vy1;
0489     uz1 = ut1 * vz1;
0490     dmu[2][0] = 0.25 * (ut1 * ut1 - ut0 * ut0) / uuu[0] / dy;
0491     dmu[2][1] = 0.25 * (ux1 * ux1 - ux0 * ux0) / uuu[1] / dy;
0492     dmu[2][2] = 0.25 * (uy1 * uy1 - uy0 * uy0) / uuu[2] / dy;
0493     dmu[2][3] = 0.25 * (uz1 * uz1 - uz0 * uz0) / uuu[3] / dy;
0494     if (fabs(0.5 * (ut1 + ut0) / uuu[0]) > UDIFF)
0495       dmu[2][0] = 0.5 * (ut1 - ut0) / dy;
0496     if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / uuu[1]) > UDIFF)
0497       dmu[2][1] = 0.5 * (ux1 - ux0) / dy;
0498     if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uuu[2]) > UDIFF)
0499       dmu[2][2] = 0.5 * (uy1 - uy0) / dy;
0500     if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uuu[3]) > UDIFF)
0501       dmu[2][3] = 0.5 * (uz1 - uz0) / dy;
0502   } else {  // matter-vacuum
0503     dmu[2][0] = dmu[2][1] = dmu[2][2] = dmu[2][3] = 0.;
0504   }
0505   // d_z u^\mu
0506   f->getCell(ix, iy, iz + 1)
0507       ->getPrimVarHCenter(eos, tau, e1, p, nb, nq, ns, vx1, vy1, vz1);
0508   f->getCell(ix, iy, iz - 1)
0509       ->getPrimVarHCenter(eos, tau, e0, p, nb, nq, ns, vx0, vy0, vz0);
0510   if (e1 > 0. && e0 > 0.) {
0511     ut0 = 1.0 / sqrt(1.0 - vx0 * vx0 - vy0 * vy0 - vz0 * vz0);
0512     ux0 = ut0 * vx0;
0513     uy0 = ut0 * vy0;
0514     uz0 = ut0 * vz0;
0515     ut1 = 1.0 / sqrt(1.0 - vx1 * vx1 - vy1 * vy1 - vz1 * vz1);
0516     ux1 = ut1 * vx1;
0517     uy1 = ut1 * vy1;
0518     uz1 = ut1 * vz1;
0519     dmu[3][0] = 0.25 * (ut1 * ut1 - ut0 * ut0) / uuu[0] / dz / (tau + 0.5 * dt);
0520     dmu[3][1] = 0.25 * (ux1 * ux1 - ux0 * ux0) / uuu[1] / dz / (tau + 0.5 * dt);
0521     dmu[3][2] = 0.25 * (uy1 * uy1 - uy0 * uy0) / uuu[2] / dz / (tau + 0.5 * dt);
0522     dmu[3][3] = 0.25 * (uz1 * uz1 - uz0 * uz0) / uuu[3] / dz / (tau + 0.5 * dt);
0523     if (fabs(0.5 * (ut1 + ut0) / uuu[0]) > UDIFF)
0524       dmu[3][0] = 0.5 * (ut1 - ut0) / dz / (tau + 0.5 * dt);
0525     if (fabs(uuu[1]) < VMIN || fabs(0.5 * (ux1 + ux0) / uuu[1]) > UDIFF)
0526       dmu[3][1] = 0.5 * (ux1 - ux0) / dz / (tau + 0.5 * dt);
0527     if (fabs(uuu[2]) < VMIN || fabs(0.5 * (uy1 + uy0) / uuu[2]) > UDIFF)
0528       dmu[3][2] = 0.5 * (uy1 - uy0) / dz / (tau + 0.5 * dt);
0529     if (fabs(uuu[3]) < VMIN || fabs(0.5 * (uz1 + uz0) / uuu[3]) > UDIFF)
0530       dmu[3][3] = 0.5 * (uz1 - uz0) / dz / (tau + 0.5 * dt);
0531   } else {  // matter-vacuum
0532     dmu[3][0] = dmu[3][1] = dmu[3][2] = dmu[3][3] = 0.;
0533   }
0534   // additional terms from Christoffel symbols :)
0535   dmu[3][0] += uuu[3] / (tau - 0.5 * dt);
0536   dmu[3][3] += uuu[0] / (tau - 0.5 * dt);
0537   // calculation of Z[mu][nu][lambda][rho]
0538   for (int i = 0; i < 4; i++)
0539     for (int j = 0; j < 4; j++)
0540       for (int k = 0; k < 4; k++)
0541         for (int l = 0; l < 4; l++) Z[i][j][k][l] = 0.0;
0542   // filling Z matrix
0543   for (int mu = 0; mu < 4; mu++)
0544     for (int nu = 0; nu < 4; nu++)
0545       for (int lam = 0; lam < 4; lam++)
0546         for (int rho = 0; rho < 4; rho++) {
0547           if (nu == rho)
0548             Z[mu][nu][lam][rho] += 0.5 * (gmunu[mu][lam] - uuu[mu] * uuu[lam]);
0549           if (mu == rho)
0550             Z[mu][nu][lam][rho] += 0.5 * (gmunu[nu][lam] - uuu[nu] * uuu[lam]);
0551           if (lam == rho)
0552             Z[mu][nu][lam][rho] -= (gmunu[mu][nu] - uuu[mu] * uuu[nu]) / 3.0;
0553         }
0554   // calculating sigma[mu][nu]
0555   for (int i = 0; i < 4; i++)
0556     for (int j = 0; j < 4; j++) {
0557       pi[i][j] = 0.0;
0558       for (int k = 0; k < 4; k++)
0559         for (int l = 0; l < 4; l++) {
0560           pi[i][j] += Z[i][j][k][l] * dmu[k][l] * 2.0 * etaS * s / 5.068;
0561         }
0562     }
0563   Pi = -zetaS * s * (dmu[0][0] + dmu[1][1] + dmu[2][2] + dmu[3][3]) /
0564        5.068;  // fm^{-4} --> GeV/fm^3
0565   du = dmu[0][0] + dmu[1][1] + dmu[2][2] + dmu[3][3];
0566   //--------- debug part: NaN/inf check, trace check, diag check, transversality
0567   //check
0568   for (int i = 0; i < 4; i++)
0569     for (int j = 0; j < 4; j++) {
0570       if (pi[i][j] != 0. && fabs(1.0 - pi[j][i] / pi[i][j]) > 1e-10)
0571         cout << "non-diag: " << pi[i][j] << "  " << pi[j][i] << endl;
0572       if (isinf(pi[i][j]) || isnan(pi[i][j])) {
0573         cout << "hydro:NSquant: inf/nan i " << i << " j " << j << endl;
0574         exit(1);
0575       }
0576     }
0577 }
0578 
0579 void Hydro::setNSvalues() {
0580   double e, p, nb, nq, ns, vx, vy, vz, piNS[4][4], PiNS;
0581   for (int ix = 0; ix < f->getNX(); ix++)
0582     for (int iy = 0; iy < f->getNY(); iy++)
0583       for (int iz = 0; iz < f->getNZ(); iz++) {
0584         Cell *c = f->getCell(ix, iy, iz);
0585         c->getPrimVar(eos, tau, e, p, nb, nq, ns, vx, vy, vz);
0586         if (e <= 0.) continue;
0587         // NSquant(ix, iy, iz, piNS, PiNS, dmu, du) ;
0588         //############## set NS values assuming initial zero flow + Bjorken z
0589         //flow
0590         double T, mub, muq, mus;
0591         double etaS, zetaS;
0592         double s =
0593             eos->s(e, nb, nq, ns);  // entropy density in the current cell
0594         eos->eos(e, nb, nq, ns, T, mub, muq, mus, p);
0595         trcoeff->getEta(e, T, etaS, zetaS);
0596         for (int i = 0; i < 4; i++)
0597           for (int j = 0; j < 4; j++) piNS[i][j] = 0.0;  // reset piNS
0598         piNS[1][1] = piNS[2][2] = 2.0 / 3.0 * etaS * s / tau / 5.068;
0599         piNS[3][3] = -2.0 * piNS[1][1];
0600         PiNS = 0.0;
0601         for (int i = 0; i < 4; i++)
0602           for (int j = 0; j <= i; j++) c->setpi(i, j, piNS[i][j]);
0603         c->setPi(PiNS);
0604       }
0605   cout << "setNS done\n";
0606 }
0607 
0608 void Hydro::ISformal() {
0609   double e, p, nb, nq, ns, vx, vy, vz, T, mub, muq, mus;
0610   double piNS[4][4], PiNS, dmu[4][4], du, pi[4][4], piH[4][4], Pi, PiH;
0611   const double gmumu[4] = {1., -1., -1., -1.};
0612 
0613   // loop #1 (relaxation+source terms)
0614   for (int ix = 0; ix < f->getNX(); ix++)
0615     for (int iy = 0; iy < f->getNY(); iy++)
0616       for (int iz = 0; iz < f->getNZ(); iz++) {
0617         Cell *c = f->getCell(ix, iy, iz);
0618         c->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vx, vy,
0619                              vz);  // instead of getPrimVar()
0620         if (e <= 0.) {             // empty cell?
0621           for (int i = 0; i < 4; i++)
0622             for (int j = 0; j <= i; j++) {
0623               c->setpiH0(i, j, 0.0);
0624               c->setpi0(i, j, 0.0);
0625             }
0626           c->setPiH0(0.0);
0627           c->setPi0(0.0);
0628         } else {  // non-empty cell
0629           // 1) relaxation(pi)+source(pi) terms for half-step
0630           double gamma = 1.0 / sqrt(1.0 - vx * vx - vy * vy - vz * vz);
0631           double u[4];
0632           u[0] = gamma;
0633           u[1] = u[0] * vx;
0634           u[2] = u[0] * vy;
0635           u[3] = u[0] * vz;
0636           // source term  + tau*delta_Q_i/delta_tau
0637           double flux[4];
0638           for(int i=0; i<4; i++)
0639            flux[i] = (tau-dt)*(c->getpi(0,i) + c->getPi()*u[0]*u[i]);
0640           flux[0] += - (tau-dt)*c->getPi();
0641           c->addFlux(flux[0], flux[1], flux[2], flux[3], 0., 0., 0.);
0642           // now calculating viscous terms in NS limit
0643           NSquant(ix, iy, iz, piNS, PiNS, dmu, du);
0644           eos->eos(e, nb, nq, ns, T, mub, muq, mus, p);
0645           //############# get relaxation times
0646           double taupi, tauPi;
0647           trcoeff->getTau(T, taupi, tauPi);
0648           //#############
0649           // relaxation term, piH,PiH-->half-step
0650           for (int i = 0; i < 4; i++)
0651             for (int j = 0; j <= i; j++) {
0652 #ifdef FORMAL_SOLUTION
0653               c->setpiH0(i, j, (c->getpi(i, j) - piNS[i][j]) *
0654                                        exp(-dt / 2.0 / gamma / taupi) +
0655                                    piNS[i][j]);
0656 #else
0657               c->setpiH0(i, j, c->getpi(i, j) - (c->getpi(i, j) - piNS[i][j]) *
0658                                                     dt / 2.0 / gamma / taupi);
0659 #endif
0660             }
0661 #ifdef FORMAL_SOLUTION
0662           c->setPiH0((c->getPi() - PiNS) * exp(-dt / 2.0 / gamma / tauPi) +
0663                      PiNS);
0664 #else
0665           c->setPiH0(c->getPi() -
0666                      (c->getPi() - PiNS) * dt / 2.0 / gamma / tauPi);
0667 #endif
0668           // sources from Christoffel symbols from \dot pi_munu
0669           double tau1 = tau - dt * 0.75;
0670           c->addpiH0(0, 0, -2. * vz * c->getpi(0, 3) / tau1 * dt /
0671                                2.);  // *gamma/gamma
0672           c->addpiH0(3, 3, -(2. * vz / tau1 * c->getpi(0, 3)) * dt / 2.);
0673           c->addpiH0(
0674               3, 0, -(vz / tau1 * c->getpi(0, 0) + vz / tau1 * c->getpi(3, 3)) *
0675                         dt / 2.);
0676           c->addpiH0(1, 0, -vz / tau1 * c->getpi(1, 3) * dt / 2.);
0677           c->addpiH0(2, 0, -vz / tau1 * c->getpi(2, 3) * dt / 2.);
0678           c->addpiH0(3, 1, -(vz / tau1 * c->getpi(0, 1)) * dt / 2.);
0679           c->addpiH0(3, 2, -(vz / tau1 * c->getpi(0, 2)) * dt / 2.);
0680           // source from full IS equations (see  draft for the description)
0681           for (int i = 0; i < 4; i++)
0682             for (int j = 0; j <= i; j++) {
0683               c->addpiH0(i, j,
0684                          -4. / 3. * c->getpi(i, j) * du / gamma * dt / 2.);
0685               for (int k = 0; k < 4;
0686                    k++)  // terms to achieve better transverality to u^mu
0687                 for (int l = 0; l < 4; l++)
0688                   c->addpiH0(i, j,
0689                              -c->getpi(i, k) * u[j] * u[l] * dmu[l][k] *
0690                                      gmumu[k] / gamma * dt / 2. -
0691                                  c->getpi(j, k) * u[i] * u[l] * dmu[l][k] *
0692                                      gmumu[k] / gamma * dt / 2.);
0693             }
0694           c->addPiH0(-4. / 3. * c->getPi() * du / gamma * dt / 2.);
0695           // 1) relaxation(piH)+source(piH) terms for full-step
0696           for (int i = 0; i < 4; i++)
0697             for (int j = 0; j <= i; j++) {
0698 #ifdef FORMAL_SOLUTION
0699               c->setpi0(i, j, (c->getpi(i, j) - piNS[i][j]) *
0700                                       exp(-dt / gamma / taupi) +
0701                                   piNS[i][j]);
0702 #else
0703               c->setpi0(i, j, c->getpi(i, j) - (c->getpiH0(i, j) - piNS[i][j]) *
0704                                                    dt / gamma / taupi);
0705 #endif
0706             }
0707 #ifdef FORMAL_SOLUTION
0708           c->setPi0((c->getPi() - PiNS) * exp(-dt / gamma / tauPi) + PiNS);
0709 #else
0710           c->setPi0(c->getPi() - (c->getPiH0() - PiNS) * dt / gamma / tauPi);
0711 #endif
0712           tau1 = tau - dt * 0.5;
0713           c->addpi0(0, 0,
0714                     -2. * vz / tau1 * c->getpiH0(0, 3) * dt);  // *gamma/gamma
0715           c->addpi0(3, 3, -(2. * vz / tau1 * c->getpiH0(0, 3)) * dt);
0716           c->addpi0(3, 0, -(vz / tau1 * c->getpiH0(0, 0) +
0717                             vz / tau1 * c->getpiH0(3, 3)) *
0718                               dt);
0719           c->addpi0(1, 0, -vz / tau1 * c->getpiH0(1, 3) * dt);
0720           c->addpi0(2, 0, -vz / tau1 * c->getpiH0(2, 3) * dt);
0721           c->addpi0(3, 1, -(vz / tau1 * c->getpiH0(0, 1)) * dt);
0722           c->addpi0(3, 2, -(vz / tau1 * c->getpiH0(0, 2)) * dt);
0723           // source from full IS equations (see draft for the description)
0724           for (int i = 0; i < 4; i++)
0725             for (int j = 0; j <= i; j++) {
0726               c->addpi0(i, j, -4. / 3. * c->getpiH0(i, j) * du / gamma * dt);
0727               for (int k = 0; k < 4;
0728                    k++)  // terms to achieve better transverality to u^mu
0729                 for (int l = 0; l < 4; l++)
0730                   c->addpi0(i, j, -c->getpiH0(i, k) * u[j] * u[l] * dmu[l][k] *
0731                                           gmumu[k] / gamma * dt -
0732                                       c->getpiH0(j, k) * u[i] * u[l] *
0733                                           dmu[l][k] * gmumu[k] / gamma * dt);
0734             }
0735           c->addPi0(-4. / 3. * c->getPiH0() * du / gamma * dt);
0736         }  // end non-empty cell
0737       }    // end loop #1
0738 
0739   // 3) -- advection ---
0740   for (int ix = 0; ix < f->getNX(); ix++)
0741     for (int iy = 0; iy < f->getNY(); iy++)
0742       for (int iz = 0; iz < f->getNZ(); iz++) {
0743         Cell *c = f->getCell(ix, iy, iz);
0744         c->getPrimVarHCenter(eos, tau - 0.5 * dt, e, p, nb, nq, ns, vx, vy,
0745                              vz);  // getPrimVar() before
0746         if (e <= 0.) continue;
0747         double xm = -vx * dt / f->getDx();
0748         double ym = -vy * dt / f->getDy();
0749         double zm = -vz * dt / f->getDz() / (tau - 0.5 * dt);
0750         double xmH = -vx * dt / f->getDx() / 2.0;
0751         double ymH = -vy * dt / f->getDy() / 2.0;
0752         double zmH = -vz * dt / f->getDz() / (tau - 0.5 * dt) / 2.0;
0753         double wx[2] = {(1. - fabs(xm)), fabs(xm)};
0754         double wy[2] = {(1. - fabs(ym)), fabs(ym)};
0755         double wz[2] = {(1. - fabs(zm)), fabs(zm)};
0756         double wxH[2] = {(1. - fabs(xmH)), fabs(xmH)};
0757         double wyH[2] = {(1. - fabs(ymH)), fabs(ymH)};
0758         double wzH[2] = {(1. - fabs(zmH)), fabs(zmH)};
0759         for (int i = 0; i < 4; i++)
0760           for (int j = 0; j < 4; j++) {
0761             pi[i][j] = piH[i][j] = 0.0;
0762           }
0763         Pi = PiH = 0.0;
0764         for (int jx = 0; jx < 2; jx++)
0765           for (int jy = 0; jy < 2; jy++)
0766             for (int jz = 0; jz < 2; jz++) {
0767               // pi,Pi-->full step, piH,PiH-->half-step
0768               Cell *c1 = f->getCell(ix + jx * sign(xm), iy + jy * sign(ym),
0769                                     iz + jz * sign(zm));
0770               for (int i = 0; i < 4; i++)
0771                 for (int j = 0; j < 4; j++) {
0772                   pi[i][j] += wx[jx] * wy[jy] * wz[jz] * c1->getpi0(i, j);
0773                   piH[i][j] += wxH[jx] * wyH[jy] * wzH[jz] * c1->getpiH0(i, j);
0774                 }
0775               Pi += wx[jx] * wy[jy] * wz[jz] * c1->getPi0();
0776               PiH += wxH[jx] * wyH[jy] * wzH[jz] * c1->getPiH0();
0777             }
0778 
0779         //--------- debug part: trace check, diag check, transversality check
0780         for (int i = 0; i < 4; i++)
0781           for (int j = 0; j < 4; j++) {
0782             if (pi[i][j] != 0. && fabs(1.0 - pi[j][i] / pi[i][j]) > 1e-10)
0783               cout << "non-diag: " << pi[i][j] << "  " << pi[j][i] << endl;
0784           }
0785         //------ end debug
0786         //======= hydro applicability check (viscous corrections limiter):
0787         // double maxT0 = max(fabs((e+p)*vx*vx/(1.-vx*vx-vy*vy-vz*vz)+p),
0788         //   fabs((e+p)*vy*vy/(1.-vx*vx-vy*vy-vz*vz)+p)) ;
0789         double maxT0 = max((e + p) / (1. - vx * vx - vy * vy - vz * vz) - p,
0790                            (e + p) * (vx * vx + vy * vy + vz * vz) /
0791                                    (1. - vx * vx - vy * vy - vz * vz) +
0792                                p);
0793         // double maxpi = max(fabs(pi[1][1]),fabs(pi[2][2])) ;
0794         double maxpi = 0.;
0795         for (int i = 0; i < 4; i++)
0796           for (int j = 0; j < 4; j++)
0797             if (fabs(pi[i][j]) > maxpi) maxpi = fabs(pi[i][j]);
0798         bool rescaled = false;
0799         if (maxT0 / maxpi < 1.0) {
0800           for (int i = 0; i < 4; i++)
0801             for (int j = 0; j < 4; j++) {
0802               pi[i][j] = 0.1 * pi[i][j] * maxT0 / maxpi;
0803               piH[i][j] = 0.1 * piH[i][j] * maxT0 / maxpi;
0804             }
0805           rescaled = true;
0806         }
0807         if (fabs(Pi) > p) {
0808           if (Pi != 0.) Pi = 0.1 * Pi / fabs(Pi) * p;
0809           if (PiH != 0.) PiH = 0.1 * PiH / fabs(PiH) * p;
0810           rescaled = true;
0811         }
0812         if (rescaled)
0813           c->setViscCorrCutFlag(maxT0 / maxpi);
0814         else
0815           c->setViscCorrCutFlag(1.);
0816         // updating to the new values
0817         for (int i = 0; i < 4; i++)
0818           for (int j = 0; j <= i; j++) {
0819             c->setpi(i, j, pi[i][j]);
0820             c->setpiH(i, j, piH[i][j]);
0821           }
0822         c->setPi(Pi);
0823         c->setPiH(PiH);
0824         // source term  - (tau+dt)*delta_Q_(i+1)/delta_tau
0825         double gamma = 1.0 / sqrt(1.0 - vx * vx - vy * vy - vz * vz);
0826         double u[4];
0827         u[0] = gamma;
0828         u[1] = u[0] * vx;
0829         u[2] = u[0] * vy;
0830         u[3] = u[0] * vz;
0831         double flux[4];
0832         for(int i=0; i<4; i++)
0833          flux[i] = - tau*(c->getpi(0,i) + c->getPi()*u[0]*u[i]);
0834         flux[0] += tau*c->getPi();
0835         c->addFlux(flux[0], flux[1], flux[2], flux[3], 0., 0., 0.);
0836       }  // advection loop (all cells)
0837 }
0838 
0839 
0840 // this procedure explicitly uses T_==0, X_==1, Y_==2, Z_==3
0841 void Hydro::visc_flux(Cell *left, Cell *right, int direction) {
0842   double flux[4];
0843   int ind2 = 0;
0844   double dxa = 0.;
0845   // exit if noth cells are not full with matter
0846   if (left->getM(direction) < 1. && right->getM(direction) < 1.) return;
0847 
0848   if (direction == X_)
0849     dxa = f->getDx();
0850   else if (direction == Y_)
0851     dxa = f->getDy();
0852   else if (direction == Z_)
0853     dxa = f->getDz() * (tau + 0.5 * dt);
0854   double e, p, nb, nq, ns, vxl, vyl, vzl, vxr, vyr, vzr;
0855   // we need to know the velocities at both cell centers at (n+1/2) in order to
0856   // interpolate to
0857   // get the value at the interface
0858   left->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vxl, vyl, vzl);
0859   right->getPrimVarHCenter(eos, tau - dt / 2., e, p, nb, nq, ns, vxr, vyr, vzr);
0860   vxl = 0.5 * (vxl + vxr);
0861   vyl = 0.5 * (vyl + vyr);
0862   vzl = 0.5 * (vzl + vzr);
0863   double v = sqrt(vxl * vxl + vyl * vyl + vzl * vzl);
0864   if (v > 1.) {
0865     vxl = 0.99 * vxl / v;
0866     vyl = 0.99 * vyl / v;
0867     vzl = 0.99 * vzl / v;
0868   }
0869   double gamma = 1. / sqrt(1. - v * v);
0870   double uuu[4] = {gamma, gamma * vxl, gamma * vyl, gamma * vzl};
0871   double gmumu[4] = {1., -1., -1., -1.};
0872   if (direction == X_)
0873     ind2 = 1;
0874   else if (direction == Y_)
0875     ind2 = 2;
0876   else if (direction == Z_)
0877     ind2 = 3;
0878   for (int ind1 = 0; ind1 < 4; ind1++) {
0879     flux[ind1] = 0.5 * (left->getpiH(ind1, ind2) + right->getpiH(ind1, ind2));
0880     if (ind1 == ind2)
0881       flux[ind1] += -0.5 * (left->getPiH() + right->getPiH()) *
0882                     gmumu[ind1];  // gmunu is diagonal
0883     flux[ind1] +=
0884         0.5 * (left->getPiH() + right->getPiH()) * uuu[ind1] * uuu[ind2];
0885   }
0886   for (int i = 0; i < 4; i++) flux[i] = flux[i] * (tau - 0.5*dt) * dt / dxa;
0887   left->addFlux(-flux[T_], -flux[X_], -flux[Y_], -flux[Z_], 0., 0., 0.);
0888   right->addFlux(flux[T_], flux[X_], flux[Y_], flux[Z_], 0., 0., 0.);
0889 }
0890 
0891 void Hydro::performStep(void) {
0892   debugRiemann = false;  // turn off debug output
0893 
0894   f->updateM(tau, dt);
0895 
0896   tau_z = dt / 2. / log(1 + dt / 2. / tau);
0897 
0898   //-----PREDICTOR-ideal
0899   for (int iy = 0; iy < f->getNY(); iy++)
0900     for (int iz = 0; iz < f->getNZ(); iz++)
0901       for (int ix = 0; ix < f->getNX(); ix++) {
0902         Cell *c = f->getCell(ix, iy, iz);
0903         c->saveQprev();
0904         c->clearFlux();
0905       }
0906   // X dir
0907   for (int iy = 0; iy < f->getNY(); iy++)
0908     for (int iz = 0; iz < f->getNZ(); iz++)
0909       for (int ix = 0; ix < f->getNX() - 1; ix++) {
0910         hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix + 1, iy, iz), X_,
0911                   PREDICT);
0912       }
0913   //    cout << "predictor X done\n" ;
0914   // Y dir
0915   for (int iz = 0; iz < f->getNZ(); iz++)
0916     for (int ix = 0; ix < f->getNX(); ix++)
0917       for (int iy = 0; iy < f->getNY() - 1; iy++) {
0918         hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy + 1, iz), Y_,
0919                   PREDICT);
0920       }
0921   //    cout << "predictor Y done\n" ;
0922   // Z dir
0923   for (int ix = 0; ix < f->getNX(); ix++)
0924     for (int iy = 0; iy < f->getNY(); iy++)
0925       for (int iz = 0; iz < f->getNZ() - 1; iz++) {
0926         hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy, iz + 1), Z_,
0927                   PREDICT);
0928       }
0929   //    cout << "predictor Z done\n" ;
0930 
0931   for (int iy = 0; iy < f->getNY(); iy++)
0932     for (int iz = 0; iz < f->getNZ(); iz++)
0933       for (int ix = 0; ix < f->getNX(); ix++) {
0934         Cell *c = f->getCell(ix, iy, iz);
0935         source_step(ix, iy, iz, PREDICT);
0936         c->updateQtoQhByFlux();
0937         c->clearFlux();
0938       }
0939 
0940   //----CORRECTOR-ideal
0941 
0942   tau_z = dt / log(1 + dt / tau);
0943   // X dir
0944   for (int iy = 0; iy < f->getNY(); iy++)
0945     for (int iz = 0; iz < f->getNZ(); iz++)
0946       for (int ix = 0; ix < f->getNX() - 1; ix++) {
0947         hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix + 1, iy, iz), X_,
0948                   CORRECT);
0949       }
0950   //    cout << "corrector X done\n" ;
0951   // Y dir
0952   for (int iz = 0; iz < f->getNZ(); iz++)
0953     for (int ix = 0; ix < f->getNX(); ix++)
0954       for (int iy = 0; iy < f->getNY() - 1; iy++) {
0955         hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy + 1, iz), Y_,
0956                   CORRECT);
0957       }
0958   //    cout << "corrector Y done\n" ;
0959   // Z dir
0960   for (int ix = 0; ix < f->getNX(); ix++)
0961     for (int iy = 0; iy < f->getNY(); iy++)
0962       for (int iz = 0; iz < f->getNZ() - 1; iz++) {
0963         hlle_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy, iz + 1), Z_,
0964                   CORRECT);
0965       }
0966   //    cout << "corrector Z done\n" ;
0967 
0968   for (int iy = 0; iy < f->getNY(); iy++)
0969     for (int iz = 0; iz < f->getNZ(); iz++)
0970       for (int ix = 0; ix < f->getNX(); ix++) {
0971         Cell *c = f->getCell(ix, iy, iz);
0972         source_step(ix, iy, iz, CORRECT);
0973         c->updateByFlux();
0974         c->clearFlux();
0975       }
0976   tau += dt;
0977   f->correctImagCells();
0978 
0979   //===== viscous part ======
0980   if (trcoeff->isViscous()) {
0981     ISformal();  // evolution of viscous quantities according to IS equations
0982 
0983     // X dir
0984     for (int iy = 0; iy < f->getNY(); iy++)
0985       for (int iz = 0; iz < f->getNZ(); iz++)
0986         for (int ix = 0; ix < f->getNX() - 1; ix++) {
0987           visc_flux(f->getCell(ix, iy, iz), f->getCell(ix + 1, iy, iz), X_);
0988         }
0989     //  cout << "visc_flux X done\n" ;
0990     // Y dir
0991     for (int iz = 0; iz < f->getNZ(); iz++)
0992       for (int ix = 0; ix < f->getNX(); ix++)
0993         for (int iy = 0; iy < f->getNY() - 1; iy++) {
0994           visc_flux(f->getCell(ix, iy, iz), f->getCell(ix, iy + 1, iz), Y_);
0995         }
0996     //  cout << "visc_flux Y done\n" ;
0997     // Z dir
0998     for (int ix = 0; ix < f->getNX(); ix++)
0999       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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!