vhlle/src/cll.cpp

0001 #include <iostream>
0002 #include <fstream>
0003 #include <cmath>
0004 #include <iomanip>
0005 #include "rmn.h"
0006 #include "cll.h"
0007
0008 using namespace std;
0009
0010 // slope limiter; chooses minimal abs of the neighbouring slopes
0011 double minmod(double a, double b) {
0012   if (a * b <= 0.) return 0.;
0013   //    else return (a*a*b+a*b*b)/(a*a+b*b) ;
0014   if (fabs(a) > fabs(b))
0015     return b;
0016   else
0017     return a;
0018 }
0019
0020 // index44: returns an index of pi^{mu nu} mu,nu component in a plain 1D array
0021 int index44(const int &i, const int &j) {
0022   if (i > 3 || j > 3 || i < 0 || j < 0) {
0023     std::cout << "index44: i j " << i << " " << j << endl;
0024     exit(1);
0025   }
0026   if (j < i)
0027     return (i * (i + 1)) / 2 + j;
0028   else
0029     return (j * (j + 1)) / 2 + i;
0030 }
0031
0032 Cell::Cell() {
0033   for (int i = 0; i < 7; i++) {
0034     Q[i] = 0.;
0035     Qh[i] = 0.;
0036     Qprev[i] = 0.;
0037     flux[i] = 0.;
0038   }
0039   viscCorrCut = 1.;
0040   for (int i = 0; i < 10; i++) {
0041     pi[i] = 0.0;
0042     piH[i] = 0.0;
0043     pi0[i] = 0.0;
0044     piH0[i] = 0.0;
0045   }
0046   Pi = 0.0;
0047   PiH = 0.0;
0048   Pi0 = 0.0;
0049   PiH0 = 0.0;
0050   setAllM(0.);
0051 }
0052
0053 void Cell::updateByFlux() {
0054   for (int i = 0; i < 7; i++) Q[i] += flux[i];
0055 }
0056
0057 void Cell::updateQtoQhByFlux() {
0058   for (int i = 0; i < 7; i++) Qh[i] = Q[i] + flux[i];
0059 }
0060
0061 void Cell::getPrimVar(EoS *eos, double tau, double &_e, double &_p, double &_nb,
0062                       double &_nq, double &_ns, double &_vx, double &_vy,
0063                       double &_vz) {
0064   double _Q[7];
0065   for (int i = 0; i < 7; i++) _Q[i] = Q[i] / tau;
0066   transformPV(eos, _Q, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0067   //-------------------- debug ---------------
0068   if (_nb != _nb) {
0069     cout << "---error in getPrimVar:\n";
0070     Dump(tau);
0071   }
0072   //------------------------------------------
0073 }
0074
0075 void Cell::getPrimVarLeft(EoS *eos, double tau, double &_e, double &_p,
0076                           double &_nb, double &_nq, double &_ns, double &_vx,
0077                           double &_vy, double &_vz, int dir) {
0078   double Qr[7], Ql[7], dQ[7];
0079
0080   next[dir - 1]->getQ(Qr);
0081   prev[dir - 1]->getQ(Ql);
0082
0083   for (int i = 0; i < 7; i++)
0084     dQ[i] = minmod((Qr[i] - Q[i]) / 2., (Q[i] - Ql[i]) / 2.);
0085
0086   for (int i = 0; i < 7; i++) Ql[i] = (Q[i] - dQ[i]) / tau;
0087   //    if(Ql[T_]<sqrt(Ql[X_]*Ql[X_]+Ql[Y_]*Ql[Y_]+Ql[Z_]*Ql[Z_]))
0088   transformPV(eos, Ql, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0089   //    else    transformPV(eos, Q, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz) ;
0090   //-------------------- debug ---------------
0091   if (_nb != _nb) {
0092     cout << "---error in getPrimVarLeft:\n";
0093     Dump(tau);
0094     next[dir - 1]->Dump(tau);
0095     prev[dir - 1]->Dump(tau);
0096   }
0097   //------------------------------------------
0098 }
0099
0100 void Cell::getPrimVarRight(EoS *eos, double tau, double &_e, double &_p,
0101                            double &_nb, double &_nq, double &_ns, double &_vx,
0102                            double &_vy, double &_vz, int dir) {
0103   double Qr[7], Ql[7], dQ[7];
0104
0105   next[dir - 1]->getQ(Qr);
0106   prev[dir - 1]->getQ(Ql);
0107
0108   for (int i = 0; i < 7; i++)
0109     dQ[i] = minmod((Qr[i] - Q[i]) / 2., (Q[i] - Ql[i]) / 2.);
0110
0111   for (int i = 0; i < 7; i++) Qr[i] = (Q[i] + dQ[i]) / tau;
0112   //    if(Qr[T_]<sqrt(Qr[X_]*Qr[X_]+Qr[Y_]*Qr[Y_]+Qr[Z_]*Qr[Z_]))
0113   transformPV(eos, Qr, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0114   //    else    transformPV(eos, Q, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz) ;
0115   //-------------------- debug ---------------
0116   if (_nb != _nb) {
0117     cout << "---error in getPrimVarRight:\n";
0118     Dump(tau);
0119     next[dir - 1]->Dump(tau);
0120     prev[dir - 1]->Dump(tau);
0121   }
0122   //------------------------------------------
0123 }
0124
0125 void Cell::getPrimVarHLeft(EoS *eos, double tau, double &_e, double &_p,
0126                            double &_nb, double &_nq, double &_ns, double &_vx,
0127                            double &_vy, double &_vz, int dir) {
0128   double Qr[7], Ql[7], dQ[7];
0129
0130   next[dir - 1]->getQh(Qr);
0131   prev[dir - 1]->getQh(Ql);
0132
0133   for (int i = 0; i < 7; i++)
0134     dQ[i] = minmod((Qr[i] - Qh[i]) / 2., (Qh[i] - Ql[i]) / 2.);
0135
0136   for (int i = 0; i < 7; i++) Ql[i] = (Qh[i] - dQ[i]) / tau;
0137   //    if(Ql[T_]<sqrt(Ql[X_]*Ql[X_]+Ql[Y_]*Ql[Y_]+Ql[Z_]*Ql[Z_]))
0138   transformPV(eos, Ql, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0139   //    else    transformPV(eos, Qh, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz) ;
0140   //-------------------- debug ---------------
0141   if (_nb != _nb) {
0142     cout << "---error in getPrimVarHLeft:\n";
0143     Dump(tau);
0144     next[dir - 1]->Dump(tau);
0145     prev[dir - 1]->Dump(tau);
0146   }
0147   //------------------------------------------
0148 }
0149
0150 void Cell::getPrimVarHRight(EoS *eos, double tau, double &_e, double &_p,
0151                             double &_nb, double &_nq, double &_ns, double &_vx,
0152                             double &_vy, double &_vz, int dir) {
0153   double Qr[7], Ql[7], dQ[7];
0154
0155   next[dir - 1]->getQh(Qr);
0156   prev[dir - 1]->getQh(Ql);
0157
0158   for (int i = 0; i < 7; i++)
0159     dQ[i] = minmod((Qr[i] - Qh[i]) / 2., (Qh[i] - Ql[i]) / 2.);
0160
0161   for (int i = 0; i < 7; i++) Qr[i] = (Qh[i] + dQ[i]) / tau;
0162   //    if(Qr[T_]<sqrt(Qr[X_]*Qr[X_]+Qr[Y_]*Qr[Y_]+Qr[Z_]*Qr[Z_]))
0163   transformPV(eos, Qr, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0164   //    else    transformPV(eos, Qh, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz) ;
0165   //-------------------- debug ---------------
0166   if (_nb != _nb) {
0167     cout << "---error in getPrimVarHRight:\n";
0168     Dump(tau);
0169     next[dir - 1]->Dump(tau);
0170     prev[dir - 1]->Dump(tau);
0171   }
0172   //------------------------------------------
0173 }
0174
0175 void Cell::getPrimVarHCenter(EoS *eos, double tau, double &_e, double &_p,
0176                              double &_nb, double &_nq, double &_ns, double &_vx,
0177                              double &_vy, double &_vz) {
0178   double _Q[7];
0179   for (int i = 0; i < 7; i++) _Q[i] = Qh[i] / tau;
0180   transformPV(eos, _Q, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0181 }
0182
0183 void Cell::getPrimVarPrev(EoS *eos, double tau, double &_e, double &_p,
0184                           double &_nb, double &_nq, double &_ns, double &_vx,
0185                           double &_vy, double &_vz) {
0186   double _Q[7];
0187   for (int i = 0; i < 7; i++) _Q[i] = Qprev[i] / tau;
0188   transformPV(eos, _Q, _e, _p, _nb, _nq, _ns, _vx, _vy, _vz);
0189 }
0190
0191
0192 void Cell::setPrimVar(EoS *eos, double tau, double _e, double _nb, double _nq,
0193                       double _ns, double _vx, double _vy, double _vz) {
0194   const double gamma2 = 1. / (1 - _vx * _vx - _vy * _vy - _vz * _vz);
0195   const double p = eos->p(_e, _nb, _nq, _ns);
0196   Q[T_] = tau * (_e + p * (_vx * _vx + _vy * _vy + _vz * _vz)) * gamma2;
0197   Q[X_] = tau * (_e + p) * _vx * gamma2;
0198   Q[Y_] = tau * (_e + p) * _vy * gamma2;
0199   Q[Z_] = tau * (_e + p) * _vz * gamma2;
0200   Q[NB_] = tau * _nb * sqrt(gamma2);
0201   Q[NQ_] = tau * _nq * sqrt(gamma2);
0202   Q[NS_] = tau * _ns * sqrt(gamma2);
0203   if (Q[NB_] != Q[NB_]) {
0204     cout << "init error!\n";
0205     eos->p(_e, _nb, _nq, _ns);
0206     cout << "e = " << _e << " p = " << p << " vx = " << _vx << " vy = " << _vy
0207          << " vz = " << _vz << " gamma2 = " << gamma2 << endl;
0208     //      exit(1) ;
0209     return;
0210   }
0211 }
0212
0213 void Cell::Dump(double tau) {
0214   cout << "---------cell values dump-------\n";
0215   cout << setw(5) << ix << setw(5) << iy << setw(5) << iz << endl;
0216   cout << setw(14) << Q[0] / tau << setw(14) << Q[1] / tau << setw(14)
0217        << Q[2] / tau << setw(14) << Q[3] / tau << endl;
0218   cout << setw(14) << Q[4] / tau << setw(14) << Q[5] / tau << setw(14)
0219        << Q[6] / tau << endl;
0220   cout << setw(14) << Qh[0] / tau << setw(14) << Qh[1] / tau << setw(14)
0221        << Qh[2] / tau << setw(14) << Qh[3] / tau << endl;
0222   cout << setw(14) << Qh[4] / tau << setw(14) << Qh[5] / tau << setw(14)
0223        << Qh[6] / tau << endl;
0224
0225   cout << "--------------------------------\n";
0226 }