vhlle/src/cll.h

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 #pragma once
0020 #include <iosfwd>
0021 #include <algorithm>
0022 #include "inc.h"
0023 class EoS;
0024
0025 // returns an index of pi^{mu nu} mu,nu component in a plain 1D array
0026 int index44(const int &i, const int &j);
0027
0028 // this class stores the information about an individual hydro cell
0029 class Cell {
0030  private:
0031   // Q usually denotes the conserved quantities, T^{0i}
0032   // here, Q, Qh, Qprev etc. ~tau*T^{0i}, like Hirano'01
0033   double Q[7];      // final values at a given timestep
0034   double Qh[7];     // half-step updated values
0035   double Qprev[7];  // values at the end of previous timestep
0036   double pi[10], piH[10];  // pi^{mu nu}, WITHOUT tau factor, final (pi) and
0037                            // half-step updated (piH)
0038   double Pi,
0039       PiH;  // Pi, WITHOUT tau factor, final (Pi) and half-step updated (PiH)
0040   double pi0[10], piH0[10];  // // pi^{mu nu}, WITHOUT tau factor, auxiliary
0041   double Pi0, PiH0;  // viscous, WITHOUT tau factor, auxiliary
0042   double flux[7];  // cumulative fluxes
0043   Cell *next[3];   // pointer to the next cell in a given direction
0044   Cell *prev[3];   // pointer to the previous cell in a given direction
0045   double m[3];     // extend of matter propagation inside cell [0...1]
0046   double dm[3];    // auxiliary
0047   int ix, iy, iz;  // cell coordinate on the grid
0048   // viscCorrCut: flag if the viscous corrections are cut for this cell:
0049   // 1.0 = uncut, < 1 :  cut by this factor
0050   double viscCorrCut;
0051
0052  public:
0053   Cell();
0054   ~Cell() {};
0055   inline void setPos(int iix, int iiy, int iiz) {
0056     ix = iix;
0057     iy = iiy;
0058     iz = iiz;
0059   }
0060   inline int getX(void) { return ix; }
0061   inline int getY(void) { return iy; }
0062   inline int getZ(void) { return iz; }
0063
0064   inline void setQ(double *_Q) {
0065     for (int i = 0; i < 7; i++) Q[i] = _Q[i];
0066     if (Q[T_] < 0.) {
0067       for (int i = 0; i < 7; i++) Q[i] = 0.;
0068     }
0069   }
0070   inline void setQh(double *_Qh) {
0071     for (int i = 0; i < 7; i++) Qh[i] = _Qh[i];
0072     if (Qh[T_] < 0.) {
0073       for (int i = 0; i < 7; i++) Qh[i] = 0.;
0074     }
0075   }
0076
0077   // getter and setter methods for the class members
0078   inline double getpi(const int &i, const int &j) { return pi[index44(i, j)]; }
0079   inline double getpiH(const int &i, const int &j) {
0080     return piH[index44(i, j)];
0081   }
0082   inline double getpi0(const int &i, const int &j) {
0083     return pi0[index44(i, j)];
0084   }
0085   inline double getpiH0(const int &i, const int &j) {
0086     return piH0[index44(i, j)];
0087   }
0088   inline double getPi(void) { return Pi; }
0089   inline double getPiH(void) { return PiH; }
0090   inline double getPi0(void) { return Pi0; }
0091   inline double getPiH0(void) { return PiH0; }
0092
0093   inline void setpi(const int &i, const int &j, const double &val) {
0094     pi[index44(i, j)] = val;
0095   }
0096   inline void setpiH(const int &i, const int &j, const double &val) {
0097     piH[index44(i, j)] = val;
0098   }
0099   inline void setpi0(const int &i, const int &j, const double &val) {
0100     pi0[index44(i, j)] = val;
0101   }
0102   inline void setpiH0(const int &i, const int &j, const double &val) {
0103     piH0[index44(i, j)] = val;
0104   }
0105   inline void addpi0(const int &i, const int &j, const double &val) {
0106     pi0[index44(i, j)] += val;
0107   }
0108   inline void addpiH0(const int &i, const int &j, const double &val) {
0109     piH0[index44(i, j)] += val;
0110   }
0111   inline void setPi(const double &val) { Pi = val; }
0112   inline void setPiH(const double &val) { PiH = val; }
0113   inline void setPi0(const double &val) { Pi0 = val; }
0114   inline void setPiH0(const double &val) { PiH0 = val; }
0115   inline void addPi0(const double &val) { Pi0 += val; }
0116   inline void addPiH0(const double &val) { PiH0 += val; }
0117
0118   inline void getQ(double *_Q) {
0119     for (int i = 0; i < 7; i++) _Q[i] = Q[i];
0120   }
0121   inline void getQh(double *_Qh) {
0122     for (int i = 0; i < 7; i++) _Qh[i] = Qh[i];
0123   }
0124   inline void saveQprev(void) {
0125     for (int i = 0; i < 7; i++) Qprev[i] = Q[i];
0126   }
0127   inline void setNext(int i, Cell *c) { next[i - 1] = c; }
0128   inline void setPrev(int i, Cell *c) { prev[i - 1] = c; }
0129   inline Cell *getNext(int i) { return next[i - 1]; }
0130   inline Cell *getPrev(int i) { return prev[i - 1]; }
0131
0132   inline void setAllM(double value) { m[0] = m[1] = m[2] = value; }
0133   inline void addM(int dir, double inc) {
0134     m[dir - 1] += inc;
0135     if (m[dir - 1] > 0.9)
0136       for (int i = 0; i < 3; i++) m[i] = 1.;
0137   }
0138   inline double getM(int dir) { return m[dir - 1]; }
0139   inline double getMaxM(void) { return std::max(m[0], std::max(m[1], m[2])); }
0140   inline void setDM(int dir, double value) { dm[dir - 1] = value; }
0141   inline double getDM(int dir) { return dm[dir - 1]; }
0142
0143   inline void setpi0(double values[4][4]) {
0144     for (int i = 0; i < 4; i++)
0145       for (int j = 0; j < 4; j++) pi0[index44(i, j)] = values[i][j];
0146   }
0147   inline void setpiH0(double values[4][4]) {
0148     for (int i = 0; i < 4; i++)
0149       for (int j = 0; j < 4; j++) piH0[index44(i, j)] = values[i][j];
0150   }
0151
0152   // get the energy density, pressure, charge densities and flow velocity
0153   // components (e,p,n,v) from conserved quantities Q in the centre of the cell
0154   void getPrimVar(EoS *eos, double tau, double &_e, double &_p, double &_nb,
0155                   double &_nq, double &_ns, double &_vx, double &_vy,
0156                   double &_vz);
0157   // (e,p,n,v) at cell's left boundary in a given direction dir
0158   void getPrimVarLeft(EoS *eos, double tau, double &_e, double &_p, double &_nb,
0159                       double &_nq, double &_ns, double &_vx, double &_vy,
0160                       double &_vz, int dir);
0161   // (e,p,n,v) at cell's right boundary in a given direction dir
0162   void getPrimVarRight(EoS *eos, double tau, double &_e, double &_p,
0163                        double &_nb, double &_nq, double &_ns, double &_vx,
0164                        double &_vy, double &_vz, int dir);
0165
0166   // (e,p,n,v) from half-step updated Qh at cell's left boundary in a given
0167   // direction
0168   void getPrimVarHLeft(EoS *eos, double tau, double &_e, double &_p,
0169                        double &_nb, double &_nq, double &_ns, double &_vx,
0170                        double &_vy, double &_vz, int dir);
0171   // (e,p,n,v) from half-step updated Qh at cell's right boundary in a given
0172   // direction
0173   void getPrimVarHRight(EoS *eos, double tau, double &_e, double &_p,
0174                         double &_nb, double &_nq, double &_ns, double &_vx,
0175                         double &_vy, double &_vz, int dir);
0176   // (e,p,n,v) from half-step updated Qh at cell's centre
0177   void getPrimVarHCenter(EoS *eos, double tau, double &_e, double &_p,
0178                          double &_nb, double &_nq, double &_ns, double &_vx,
0179                          double &_vy, double &_vz);
0180   // (e,p,n,v) at the previous timestep and cell's centre
0181   void getPrimVarPrev(EoS *eos, double tau, double &_e, double &_p, double &_nb,
0182                       double &_nq, double &_ns, double &_vx, double &_vy,
0183                       double &_vz);
0184   // calculate and set Q from (e,n,v)
0185   void setPrimVar(EoS *eos, double tau, double _e, double _nb, double _nq,
0186                   double _ns, double _vx, double _vy, double _vz);
0187
0188   // update the cumulative fluxes through the cell
0189   inline void addFlux(double Ft, double Fx, double Fy, double Fz, double Fnb,
0190                       double Fnq, double Fns) {
0191     flux[T_] += Ft;
0192     flux[X_] += Fx;
0193     flux[Y_] += Fy;
0194     flux[Z_] += Fz;
0195     flux[NB_] += Fnb;
0196     flux[NQ_] += Fnq;
0197     flux[NS_] += Fns;
0198   }
0199   inline void clearFlux(void) {
0200     for (int i = 0; i < 7; i++) flux[i] = 0.;
0201   }
0202   void updateByFlux();       // Q = Q + flux
0203   void updateQtoQhByFlux();  // Qh = Q + flux
0204   inline void setViscCorrCutFlag(double value) { viscCorrCut = value; }
0205   inline double getViscCorrCutFlag(void) { return viscCorrCut; }
0206   void Dump(double tau);  // dump the contents of the cell into dump.dat
0207 };