Back to home page

sPhenix code displayed by LXR
 
 

    

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

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 <cmath>
0023 #include <algorithm>
0024 #include "inc.h"
0025 #include "rmn.h"
0026 #include "fld.h"
0027 #include "cll.h"
0028 #include "eos.h"
0029 #include "trancoeff.h"
0030 
0031 #define OUTPI
0032 
0033 using namespace std;
0034 
0035 // returns the velocities in cartesian coordinates, fireball rest frame.
0036 // Y=longitudinal rapidity of fluid
0037 void Fluid::getCMFvariables(Cell *c, double tau, double &e, double &nb,
0038                             double &nq, double &ns, double &vx, double &vy,
0039                             double &Y) {
0040   double p, vz;
0041   c->getPrimVar(eos, tau, e, p, nb, nq, ns, vx, vy, vz);
0042   double eta = getZ(c->getZ());
0043   //    Y = eta + TMath::ATanH(vz) ;
0044   Y = eta + 1. / 2. * log((1. + vz) / (1. - vz));
0045   vx = vx * cosh(Y - eta) / cosh(Y);
0046   vy = vy * cosh(Y - eta) / cosh(Y);
0047 }
0048 
0049 Fluid::Fluid(EoS *_eos, EoS *_eosH, TransportCoeff *_trcoeff, int _nx, int _ny,
0050              int _nz, double _minx, double _maxx, double _miny, double _maxy,
0051              double _minz, double _maxz, double _dt, double eCrit) {
0052   eos = _eos;
0053   eosH = _eosH;
0054   trcoeff = _trcoeff;
0055   nx = _nx;
0056   ny = _ny;
0057   nz = _nz;
0058   minx = _minx;
0059   maxx = _maxx;
0060   miny = _miny;
0061   maxy = _maxy;
0062   minz = _minz;
0063   maxz = _maxz;
0064   dx = (maxx - minx) / (nx - 1);
0065   dy = (maxy - miny) / (ny - 1);
0066   dz = (maxz - minz) / (nz - 1);
0067   dt = _dt;
0068 
0069   cell = new Cell[nx * ny * nz];
0070 
0071   cell0 = new Cell;
0072   cell0->setPrimVar(eos, 1.0, 0., 0., 0., 0., 0., 0.,
0073                     0.);  // tau is not important here, since *0
0074   cell0->setAllM(0.);
0075 
0076   for (int ix = 0; ix < nx; ix++)
0077     for (int iy = 0; iy < ny; iy++)
0078       for (int iz = 0; iz < nz; iz++) {
0079         getCell(ix, iy, iz)->setPrev(X_, getCell(ix - 1, iy, iz));
0080         getCell(ix, iy, iz)->setNext(X_, getCell(ix + 1, iy, iz));
0081         getCell(ix, iy, iz)->setPrev(Y_, getCell(ix, iy - 1, iz));
0082         getCell(ix, iy, iz)->setNext(Y_, getCell(ix, iy + 1, iz));
0083         getCell(ix, iy, iz)->setPrev(Z_, getCell(ix, iy, iz - 1));
0084         getCell(ix, iy, iz)->setNext(Z_, getCell(ix, iy, iz + 1));
0085         getCell(ix, iy, iz)->setPos(ix, iy, iz);
0086       }
0087 
0088   output_nt = 0;
0089   output_nx = 0;
0090   output_ny = 0;
0091 }
0092 
0093 Fluid::~Fluid() {
0094   delete[] cell;
0095   delete cell0;
0096 }
0097 
0098 void Fluid::initOutput(char *dir, int maxstep, double tau0, int cmpr2dOut) {
0099   //    directory = dir ;
0100   compress2dOut = cmpr2dOut;
0101   cout << "maxstep=" << maxstep << endl;
0102   char command[255];
0103   sprintf(command, "mkdir -p %s", dir);
0104   int return_mkdir = system(command);
0105   cout << "mkdir returns: " << return_mkdir << endl;
0106   string outx = dir;
0107   outx.append("/outx.dat");
0108   string outxvisc = dir;
0109   outxvisc.append("/outx.visc.dat");
0110   string outyvisc = dir;
0111   outyvisc.append("/outy.visc.dat");
0112   string outdiagvisc = dir;
0113   outdiagvisc.append("/diag.visc.dat");
0114   string outy = dir;
0115   outy.append("/outy.dat");
0116   string outdiag = dir;
0117   outdiag.append("/outdiag.dat");
0118   string outz = dir;
0119   outz.append("/outz.dat");
0120   string outaniz = dir;
0121   outaniz.append("/out.aniz.dat");
0122   string out2d = dir;
0123   out2d.append("/out2D.dat");
0124   string outfreeze = dir;
0125   outfreeze.append("/freezeout.dat");
0126   foutx.open(outx.c_str());
0127   fouty.open(outy.c_str());
0128   foutz.open(outz.c_str());
0129   foutdiag.open(outdiag.c_str());
0130   fout2d.open(out2d.c_str());
0131   foutxvisc.open(outxvisc.c_str());
0132   foutyvisc.open(outyvisc.c_str());
0133   foutdiagvisc.open(outdiagvisc.c_str());
0134   fout_aniz.open(outaniz.c_str());
0135   ffreeze.open(outfreeze.c_str());
0136   //################################################################
0137   // important remark. for correct diagonal output, nx=ny must hold.
0138   //################################################################
0139   foutxvisc << maxstep + 1 << "  " << getNX() << endl;
0140   foutyvisc << maxstep + 1 << "  " << getNY() << endl;
0141   foutdiagvisc << maxstep + 1 << "  " << getNX() << endl;
0142   foutx << "# " << maxstep + 1 << "  " << getNX() << endl;
0143   fouty << "# " << maxstep + 1 << "  " << getNY() << endl;
0144   foutdiag << "# " << maxstep + 1 << "  " << getNX() << endl;
0145   foutz << "# " << maxstep + 1 << "  " << getNZ() << endl;
0146   fout2d << " " << maxstep + 1 << "  " << (getNX() - 5) + 1 << "  "
0147          << (getNY() - 5) + 1 << endl;
0148   fout2d << tau0 << "  " << tau0 + 0.05 * maxstep << "  " << getX(2) << "  "
0149          << getX(getNX() - 3) << "  " << getY(2) << "  " << getY(getNY() - 3)
0150          << endl;
0151   outputGnuplot(tau0);
0152   fout_aniz << "#  tau  <<v_T>>  e_p  e'_p  (to compare with SongHeinz)\n";
0153 }
0154 
0155 void Fluid::correctImagCells(void) {
0156   double Q[7];
0157   // Z
0158   for (int ix = 0; ix < nx; ix++)
0159     for (int iy = 0; iy < ny; iy++) {
0160       // left boundary
0161       getCell(ix, iy, 2)->getQ(Q);
0162       getCell(ix, iy, 1)->setQ(Q);
0163       getCell(ix, iy, 0)->setQ(Q);
0164       // right boundary
0165       getCell(ix, iy, nz - 3)->getQ(Q);
0166       getCell(ix, iy, nz - 2)->setQ(Q);
0167       getCell(ix, iy, nz - 1)->setQ(Q);
0168     }
0169   // Y
0170   for (int ix = 0; ix < nx; ix++)
0171     for (int iz = 0; iz < nz; iz++) {
0172       // left boundary
0173       getCell(ix, 2, iz)->getQ(Q);
0174       getCell(ix, 1, iz)->setQ(Q);
0175       getCell(ix, 0, iz)->setQ(Q);
0176       // right boundary
0177       getCell(ix, ny - 3, iz)->getQ(Q);
0178       getCell(ix, ny - 2, iz)->setQ(Q);
0179       getCell(ix, ny - 1, iz)->setQ(Q);
0180     }
0181   // X
0182   for (int iy = 0; iy < ny; iy++)
0183     for (int iz = 0; iz < nz; iz++) {
0184       // left boundary
0185       getCell(2, iy, iz)->getQ(Q);
0186       getCell(1, iy, iz)->setQ(Q);
0187       getCell(0, iy, iz)->setQ(Q);
0188       // right boundary
0189       getCell(nx - 3, iy, iz)->getQ(Q);
0190       getCell(nx - 2, iy, iz)->setQ(Q);
0191       getCell(nx - 1, iy, iz)->setQ(Q);
0192     }
0193 }
0194 
0195 void Fluid::correctImagCellsFull(void) {
0196   double Q[7], _pi[4][4], _Pi;
0197   // Z
0198   for (int ix = 0; ix < nx; ix++)
0199     for (int iy = 0; iy < ny; iy++) {
0200       // left boundary
0201       getCell(ix, iy, 2)->getQ(Q);
0202       getCell(ix, iy, 1)->setQ(Q);
0203       getCell(ix, iy, 0)->setQ(Q);
0204       for (int i = 0; i < 4; i++)
0205         for (int j = 0; j < 4; j++) _pi[i][j] = getCell(ix, iy, 2)->getpi(i, j);
0206       _Pi = getCell(ix, iy, 2)->getPi();
0207 
0208       for (int i = 0; i < 4; i++)
0209         for (int j = 0; j <= i; j++) {
0210           getCell(ix, iy, 0)->setpi(i, j, _pi[i][j]);
0211           getCell(ix, iy, 1)->setpi(i, j, _pi[i][j]);
0212         }
0213       getCell(ix, iy, 0)->setPi(_Pi);
0214       getCell(ix, iy, 1)->setPi(_Pi);
0215       // right boundary
0216       getCell(ix, iy, nz - 3)->getQ(Q);
0217       getCell(ix, iy, nz - 2)->setQ(Q);
0218       getCell(ix, iy, nz - 1)->setQ(Q);
0219       for (int i = 0; i < 4; i++)
0220         for (int j = 0; j < 4; j++)
0221           _pi[i][j] = getCell(ix, iy, nz - 3)->getpi(i, j);
0222       _Pi = getCell(ix, iy, nz - 3)->getPi();
0223 
0224       for (int i = 0; i < 4; i++)
0225         for (int j = 0; j <= i; j++) {
0226           getCell(ix, iy, nz - 2)->setpi(i, j, _pi[i][j]);
0227           getCell(ix, iy, nz - 1)->setpi(i, j, _pi[i][j]);
0228         }
0229       getCell(ix, iy, nz - 2)->setPi(_Pi);
0230       getCell(ix, iy, nz - 1)->setPi(_Pi);
0231     }
0232   // Y
0233   for (int ix = 0; ix < nx; ix++)
0234     for (int iz = 0; iz < nz; iz++) {
0235       // left boundary
0236       getCell(ix, 2, iz)->getQ(Q);
0237       getCell(ix, 1, iz)->setQ(Q);
0238       getCell(ix, 0, iz)->setQ(Q);
0239       for (int i = 0; i < 4; i++)
0240         for (int j = 0; j < 4; j++) _pi[i][j] = getCell(ix, 2, iz)->getpi(i, j);
0241       _Pi = getCell(ix, 2, iz)->getPi();
0242 
0243       for (int i = 0; i < 4; i++)
0244         for (int j = 0; j <= i; j++) {
0245           getCell(ix, 0, iz)->setpi(i, j, _pi[i][j]);
0246           getCell(ix, 1, iz)->setpi(i, j, _pi[i][j]);
0247         }
0248       getCell(ix, 0, iz)->setPi(_Pi);
0249       getCell(ix, 1, iz)->setPi(_Pi);
0250       // right boundary
0251       getCell(ix, ny - 3, iz)->getQ(Q);
0252       getCell(ix, ny - 2, iz)->setQ(Q);
0253       getCell(ix, ny - 1, iz)->setQ(Q);
0254       for (int i = 0; i < 4; i++)
0255         for (int j = 0; j < 4; j++)
0256           _pi[i][j] = getCell(ix, ny - 3, iz)->getpi(i, j);
0257       _Pi = getCell(ix, ny - 3, iz)->getPi();
0258 
0259       for (int i = 0; i < 4; i++)
0260         for (int j = 0; j <= i; j++) {
0261           getCell(ix, ny - 2, iz)->setpi(i, j, _pi[i][j]);
0262           getCell(ix, ny - 1, iz)->setpi(i, j, _pi[i][j]);
0263         }
0264       getCell(ix, ny - 2, iz)->setPi(_Pi);
0265       getCell(ix, ny - 1, iz)->setPi(_Pi);
0266     }
0267   // X
0268   for (int iy = 0; iy < ny; iy++)
0269     for (int iz = 0; iz < nz; iz++) {
0270       // left boundary
0271       getCell(2, iy, iz)->getQ(Q);
0272       getCell(1, iy, iz)->setQ(Q);
0273       getCell(0, iy, iz)->setQ(Q);
0274       for (int i = 0; i < 4; i++)
0275         for (int j = 0; j < 4; j++) _pi[i][j] = getCell(2, iy, iz)->getpi(i, j);
0276       _Pi = getCell(2, iy, iz)->getPi();
0277 
0278       for (int i = 0; i < 4; i++)
0279         for (int j = 0; j <= i; j++) {
0280           getCell(0, iy, iz)->setpi(i, j, _pi[i][j]);
0281           getCell(1, iy, iz)->setpi(i, j, _pi[i][j]);
0282         }
0283       getCell(0, iy, iz)->setPi(_Pi);
0284       getCell(1, iy, iz)->setPi(_Pi);
0285       // right boundary
0286       getCell(nx - 3, iy, iz)->getQ(Q);
0287       getCell(nx - 2, iy, iz)->setQ(Q);
0288       getCell(nx - 1, iy, iz)->setQ(Q);
0289       for (int i = 0; i < 4; i++)
0290         for (int j = 0; j < 4; j++)
0291           _pi[i][j] = getCell(nx - 3, iy, iz)->getpi(i, j);
0292       _Pi = getCell(nx - 3, iy, iz)->getPi();
0293 
0294       for (int i = 0; i < 4; i++)
0295         for (int j = 0; j <= i; j++) {
0296           getCell(nx - 2, iy, iz)->setpi(i, j, _pi[i][j]);
0297           getCell(nx - 1, iy, iz)->setpi(i, j, _pi[i][j]);
0298         }
0299       getCell(nx - 2, iy, iz)->setPi(_Pi);
0300       getCell(nx - 1, iy, iz)->setPi(_Pi);
0301     }
0302 }
0303 
0304 void Fluid::updateM(double tau, double dt) {
0305   for (int ix = 0; ix < getNX(); ix++)
0306     for (int iy = 0; iy < getNY(); iy++)
0307       for (int iz = 0; iz < getNZ(); iz++) {
0308         Cell *c = getCell(ix, iy, iz);
0309         c->setDM(X_, 0.);
0310         c->setDM(Y_, 0.);
0311         c->setDM(Z_, 0.);
0312         if (getCell(ix, iy, iz)->getMaxM() < 1.) {
0313           if (getCell(ix + 1, iy, iz)->getM(X_) >= 1. ||
0314               getCell(ix - 1, iy, iz)->getM(X_) >= 1.)
0315             c->setDM(X_, dt / dx);
0316           if (getCell(ix, iy + 1, iz)->getM(Y_) >= 1. ||
0317               getCell(ix, iy - 1, iz)->getM(Y_) >= 1.)
0318             c->setDM(Y_, dt / dy);
0319           if (getCell(ix, iy, iz + 1)->getM(Z_) >= 1. ||
0320               getCell(ix, iy, iz - 1)->getM(Z_) >= 1.)
0321             c->setDM(Z_, dt / dz / tau);
0322 
0323           if (c->getDM(X_) == 0. && c->getDM(Y_) == 0.) {
0324             if (getCell(ix + 1, iy + 1, iz)->getMaxM() >= 1. ||
0325                 getCell(ix + 1, iy - 1, iz)->getMaxM() >= 1. ||
0326                 getCell(ix - 1, iy + 1, iz)->getMaxM() >= 1. ||
0327                 getCell(ix - 1, iy - 1, iz)->getMaxM() >= 1.) {
0328               c->setDM(X_, 0.707 * dt / dx);
0329               c->setDM(Y_, 0.707 * dt / dy);
0330             }
0331           }
0332         }  // if
0333       }
0334 
0335   for (int ix = 0; ix < getNX(); ix++)
0336     for (int iy = 0; iy < getNY(); iy++)
0337       for (int iz = 0; iz < getNZ(); iz++) {
0338         Cell *c = getCell(ix, iy, iz);
0339         c->addM(X_, c->getDM(X_));
0340         c->addM(Y_, c->getDM(Y_));
0341         c->addM(Z_, c->getDM(Z_));
0342       }
0343 }
0344 
0345 
0346 void Fluid::outputGnuplot(double tau) {
0347   double e, p, nb, nq, ns, t, mub, muq, mus, vx, vy, vz;
0348 
0349   // X direction
0350   for (int ix = 0; ix < nx; ix++) {
0351     double x = getX(ix);
0352     Cell *c = getCell(ix, ny / 2, nz / 2);
0353     getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
0354     eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
0355     foutx << setw(14) << tau << setw(14) << x << setw(14) << vx << setw(14)
0356           << vy << setw(14) << e << setw(14) << nb << setw(14) << t << setw(14)
0357           << mub;
0358     foutx << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
0359           << setw(14) << c->getpi(0, 2);
0360     foutx << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
0361           << setw(14) << c->getpi(1, 2);
0362     foutx << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
0363           << setw(14) << c->getpi(2, 3);
0364     foutx << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
0365           << c->getViscCorrCutFlag() << endl;
0366   }
0367   foutx << endl;
0368 
0369   // Y direction
0370   for (int iy = 0; iy < ny; iy++) {
0371     double y = getY(iy);
0372     Cell *c = getCell(nx / 2, iy, nz / 2);
0373     getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
0374     eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
0375     fouty << setw(14) << tau << setw(14) << y << setw(14) << vy << setw(14)
0376           << vx << setw(14) << e << setw(14) << nb << setw(14) << t << setw(14)
0377           << mub;
0378     fouty << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
0379           << setw(14) << c->getpi(0, 2);
0380     fouty << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
0381           << setw(14) << c->getpi(1, 2);
0382     fouty << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
0383           << setw(14) << c->getpi(2, 3);
0384     fouty << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
0385           << c->getViscCorrCutFlag() << endl;
0386   }
0387   fouty << endl;
0388 
0389   // diagonal
0390   for (int ix = 0; ix < nx; ix++) {
0391     double x = getY(ix);
0392     Cell *c = getCell(ix, ix, nz / 2);
0393     getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
0394     eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
0395     foutdiag << setw(14) << tau << setw(14) << sqrt(2.) * x << setw(14) << vx
0396              << setw(14) << vy << setw(14) << e << setw(14) << nb << setw(14)
0397              << t << setw(14) << mub << endl;
0398     foutdiag << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
0399              << setw(14) << c->getpi(0, 2);
0400     foutdiag << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
0401              << setw(14) << c->getpi(1, 2);
0402     foutdiag << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
0403              << setw(14) << c->getpi(2, 3);
0404     foutdiag << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
0405              << c->getViscCorrCutFlag() << endl;
0406   }
0407   foutdiag << endl;
0408 
0409   // Z direction
0410   for (int iz = 0; iz < nz; iz++) {
0411     double z = getZ(iz);
0412     Cell *c = getCell(nx / 2, ny / 2, iz);
0413     getCMFvariables(getCell(nx / 2, ny / 2, iz), tau, e, nb, nq, ns, vx, vy,
0414                     vz);
0415     eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
0416     foutz << setw(14) << tau << setw(14) << z << setw(14) << vz << setw(14)
0417           << vx << setw(14) << e << setw(14) << nb << setw(14) << t << setw(14)
0418           << mub;
0419     foutz << setw(14) << c->getpi(0, 0) << setw(14) << c->getpi(0, 1)
0420           << setw(14) << c->getpi(0, 2);
0421     foutz << setw(14) << c->getpi(0, 3) << setw(14) << c->getpi(1, 1)
0422           << setw(14) << c->getpi(1, 2);
0423     foutz << setw(14) << c->getpi(1, 3) << setw(14) << c->getpi(2, 2)
0424           << setw(14) << c->getpi(2, 3);
0425     foutz << setw(14) << c->getpi(3, 3) << setw(14) << c->getPi() << setw(14)
0426           << c->getViscCorrCutFlag() << endl;
0427   }
0428   foutz << endl;
0429 
0430   // averaged quantities for y=0
0431   double T0xx = 0.0, T0yy = 0.0, Txx = 0.0, Tyy = 0.0, vtNum = 0.0, vtDen = 0.0;
0432   for (int ix = 0; ix < nx; ix++)
0433     for (int iy = 0; iy < ny; iy++) {
0434       Cell *c = getCell(ix, iy, nz / 2);
0435       getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
0436       eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
0437       double dTxx =
0438           (e + p) * vx * vx / (1.0 - vx * vx - vy * vy - pow(tanh(vz), 2)) + p;
0439       double dTyy =
0440           (e + p) * vy * vy / (1.0 - vx * vx - vy * vy - pow(tanh(vz), 2)) + p;
0441       T0xx += dTxx;
0442       T0yy += dTyy;
0443       Txx += dTxx + c->getpi(1, 1);
0444       Tyy += dTyy + c->getpi(2, 2);
0445       vtNum += sqrt(vx * vx + vy * vy) * e /
0446                sqrt(1.0 - vx * vx - vy * vy - pow(tanh(vz), 2));
0447       vtDen += e / sqrt(1.0 - vx * vx - vy * vy - pow(tanh(vz), 2));
0448     }
0449   fout_aniz << setw(14) << tau << setw(14) << vtNum / vtDen << setw(14)
0450             << (T0xx - T0yy) / (T0xx + T0yy) << setw(14)
0451             << (Txx - Tyy) / (Txx + Tyy) << endl;
0452 }
0453 
0454 // unput: geom. rapidity + velocities in Bjorken frame, --> output: velocities
0455 // in lab.frame
0456 void transformToLab(double eta, double &vx, double &vy, double &vz) {
0457   const double Y = eta + 1. / 2. * log((1. + vz) / (1. - vz));
0458   vx = vx * cosh(Y - eta) / cosh(Y);
0459   vy = vy * cosh(Y - eta) / cosh(Y);
0460   vz = tanh(Y);
0461 }
0462 
0463 void Fluid::calcTotals(double tau) {
0464   double e, p, nb, nq, ns, t, mub, muq, mus, vx, vy, vz, Q[7];
0465   double E = 0., Efull = 0., Px = 0., Nb1 = 0., Nb2 = 0., S = 0.;
0466   double eta = 0;
0467 
0468   fout2d << endl;
0469   for (int ix = 2; ix < nx - 2; ix++)
0470     for (int iy = 2; iy < ny - 2; iy++)
0471       for (int iz = 2; iz < nz - 2; iz++) {
0472         Cell *c = getCell(ix, iy, iz);
0473         getCMFvariables(c, tau, e, nb, nq, ns, vx, vy, vz);
0474         c->getQ(Q);
0475         eos->eos(e, nb, nq, ns, t, mub, muq, mus, p);
0476         double s = eos->s(e, nb, nq, ns);
0477         eta = getZ(iz);
0478         E += tau * (e + p) / (1. - vx * vx - vy * vy - tanh(vz) * tanh(vz)) *
0479                  (cosh(eta) - tanh(vz) * sinh(eta)) -
0480              tau * p * cosh(eta);
0481         Nb1 += Q[NB_];
0482         Nb2 += tau * nb * (cosh(eta) - tanh(vz) * sinh(eta)) /
0483                sqrt(1. - vx * vx - vy * vy - tanh(vz) * tanh(vz));
0484         // -- noneq. corrections to entropy flux
0485         const double gmumu[4] = {1., -1., -1., -1.};
0486         double deltas = 0.;
0487         for (int i = 0; i < 4; i++)
0488           for (int j = 0; j < 4; j++)
0489             deltas += pow(c->getpi(i, j), 2) * gmumu[i] * gmumu[j];
0490         if (t > 0.05) s += 1.5 * deltas / ((e + p) * t);
0491         S += tau * s * (cosh(eta) - tanh(vz) * sinh(eta)) /
0492              sqrt(1. - vx * vx - vy * vy - tanh(vz) * tanh(vz));
0493         Efull +=
0494             tau * (e + p) / (1. - vx * vx - vy * vy - tanh(vz) * tanh(vz)) *
0495                 (cosh(eta) - tanh(vz) * sinh(eta)) -
0496             tau * p * cosh(eta);
0497         if (trcoeff->isViscous())
0498           Efull += tau * c->getpi(0, 0) * cosh(eta) +
0499                    tau * c->getpi(0, 3) * sinh(eta);
0500       }
0501   E = E * dx * dy * dz;
0502   Efull = Efull * dx * dy * dz;
0503   Nb1 *= dx * dy * dz;
0504   Nb2 *= dx * dy * dz;
0505   S *= dx * dy * dz;
0506   cout << setw(16) << "calcTotals: E = " << setw(14) << E
0507        << "  Efull = " << setw(14) << Efull << endl;
0508   cout << setw(16) << "Px = " << setw(14) << Px << "      S = " << setw(14) << S
0509        << endl;
0510 }
[ 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!