vhlle/src/rmn.cpp

0001 #include <iostream>
0002 #include <fstream>
0003 #include <iomanip>
0004 #include <cstdlib>
0005 #include <execinfo.h>
0006 #include <signal.h>
0007 #include "rmn.h"
0008
0009 using namespace std;
0010
0011 // flag optionally needed to track
0012 // problems in the transformation procedures
0013 bool debugRiemann;
0014
0015 void handler(int sig) {
0016   void *array[10];
0017   size_t size;
0018
0019   // get void*'s for all entries on the stack
0020   size = backtrace(array, 10);
0021
0022   // print out all the frames to stderr
0023   cout << "Error: signal " << sig << endl;
0024   backtrace_symbols_fd(array, size, 2);
0025   exit(1);
0026 }
0027
0028 void transformPV(EoS *eos, double Q[7], double &e, double &p, double &nb,
0029                  double &nq, double &ns, double &vx, double &vy, double &vz) {
0030   // conserved -> primitive transtormation requires
0031   // a numerical solution to 1D nonlinear algebraic equation:
0032   // v = M / ( Q_t + p(Q_t-M*v, n) )       (A.2)
0033   // M being the modulo of the vector {Q_x, Q_y, Q_z}.
0034   // Bisection/Newton methods are used to solve the equation.
0035   const int MAXIT = 100;        // maximum number of iterations
0036   const double dpe = 1. / 3.;   // dp/de estimate for Newton method
0037   const double corrf = 0.9999;  // corrected value of M
0038   // when it brakes the speed of light limit, M>Q_t
0039   double v, vl = 0., vh = 1., dvold, dv, f, df;
0040   if (debugRiemann) {
0041     cout << "transformPV debug---------------\n";
0042     cout << setw(14) << Q[0] << setw(14) << Q[1] << setw(14) << Q[2] << setw(14)
0043          << Q[3] << endl;
0044     cout << setw(14) << Q[4] << setw(14) << Q[5] << setw(14) << Q[6] << endl;
0045   }
0046   double M = sqrt(Q[X_] * Q[X_] + Q[Y_] * Q[Y_] + Q[Z_] * Q[Z_]);
0047   if (Q[T_] <= 0.) {
0048     e = 0.;
0049     p = 0.;
0050     vx = vy = vz = 0.;
0051     nb = nq = ns = 0.;
0052     return;
0053   }
0054   if (M == 0.) {
0055     e = Q[T_];
0056     vx = 0.;
0057     vy = 0.;
0058     vz = 0.;
0059     nb = Q[NB_];
0060     nq = Q[NQ_];
0061     ns = Q[NS_];
0062     p = eos->p(e, nb, nq, ns);
0063     return;
0064   }
0065   if (M > Q[T_]) {
0066     Q[X_] *= corrf * Q[T_] / M;
0067     Q[Y_] *= corrf * Q[T_] / M;
0068     Q[Z_] *= corrf * Q[T_] / M;
0069     M = Q[T_] * corrf;
0070   }
0071
0072   v = 0.5 * (vl + vh);
0073   e = Q[T_] - M * v;
0074   if (e < 0.) e = 0.;
0075   nb = Q[NB_] * sqrt(1 - v * v);
0076   nq = Q[NQ_] * sqrt(1 - v * v);
0077   ns = Q[NS_] * sqrt(1 - v * v);
0078   p = eos->p(e, nb, nq, ns);
0079   f = (Q[T_] + p) * v - M;
0080   df = (Q[T_] + p) - M * v * dpe;
0081   dvold = vh - vl;
0082   dv = dvold;
0083   for (int i = 0; i < MAXIT; i++) {
0084     if ((f + df * (vh - v)) * (f + df * (vl - v)) > 0. ||
0085         fabs(2. * f) > fabs(dvold * df)) {  // bisection
0086       dvold = dv;
0087       dv = 0.5 * (vh - vl);
0088       v = vl + dv;
0089       //            cout << "BISECTION v = " << setw(12) << v << endl ;
0090     } else {  // Newton
0091       dvold = dv;
0092       dv = f / df;
0093       v -= dv;
0094       //            cout << "NEWTON v = " << setw(12) << v << endl ;
0095     }
0096     if (fabs(dv) < 0.00001) break;
0097
0098     e = Q[T_] - M * v;
0099     if (e < 0.) e = 0.;
0100     nb = Q[NB_] * sqrt(1 - v * v);
0101     nq = Q[NQ_] * sqrt(1 - v * v);
0102     ns = Q[NS_] * sqrt(1 - v * v);
0103     p = eos->p(e, nb, nq, ns);
0104     f = (Q[T_] + p) * v - M;
0105     df = (Q[T_] + p) - M * v * dpe;
0106
0107     if (f > 0.)
0108       vh = v;
0109     else
0110       vl = v;
0111     if (debugRiemann)
0112       cout << "step " << i << "  " << e << "  " << nb << "  " << nq << "  "
0113            << ns << "  " << p << "  " << vx << "  " << vy << "  " << vz << endl;
0114     //      if(i==40) { cout << "error : does not converge\n" ; exit(1) ; } ;
0115   }  // for loop
0116      //----------after
0117      // v = 0.5*(vh+vl) ;
0118   vx = v * Q[X_] / M;
0119   vy = v * Q[Y_] / M;
0120   vz = v * Q[Z_] / M;
0121   e = Q[T_] - M * v;
0122   p = eos->p(e, nb, nq, ns);
0123   nb = Q[NB_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
0124   nq = Q[NQ_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
0125   ns = Q[NS_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
0126   if (e < 0. || sqrt(vx * vx + vy * vy + vz * vz) > 1.) {
0127     cout << Q[T_] << "  " << Q[X_] << "  " << Q[Y_] << "  " << Q[Z_] << "  "
0128          << Q[NB_] << endl;
0129     cout << "transformRF::Error\n";
0130   }
0131   if (!(nb < 0. || nb >= 0.)) {
0132     cout << "transformRF::Error nb=#ind\n";
0133     //      return ;
0134   }
0135 }
0136
0137 void transformPVBulk(EoS *eos, double Pi, double Q[7], double &e, double &p,
0138                      double &nb, double &nq, double &ns, double &vx, double &vy,
0139                      double &vz) {
0140   const int MAXIT = 100;
0141   const double dpe = 1. / 3.;
0142   const double corrf = 0.9999;
0143   double v, vl = 0., vh = 1., dvold, dv, f, df;
0144   if (debugRiemann) {
0145     cout << "transformPVBulk debug---------------\n";
0146     cout << setw(14) << Q[0] << setw(14) << Q[1] << setw(14) << Q[2] << setw(14)
0147          << Q[3] << endl;
0148     cout << setw(14) << Q[4] << setw(14) << Q[5] << setw(14) << Q[6] << endl;
0149   }
0150   double M = sqrt(Q[X_] * Q[X_] + Q[Y_] * Q[Y_] + Q[Z_] * Q[Z_]);
0151   if (Q[T_] <= 0.) {
0152     e = 0.;
0153     p = 0.;
0154     vx = vy = vz = 0.;
0155     nb = nq = ns = 0.;
0156     return;
0157   }
0158   if (M == 0.) {
0159     e = Q[T_];
0160     vx = 0.;
0161     vy = 0.;
0162     vz = 0.;
0163     nb = Q[NB_];
0164     nq = Q[NQ_];
0165     ns = Q[NS_];
0166     p = eos->p(e, nb, nq, ns) + Pi;
0167     return;
0168   }
0169   if (M > Q[T_]) {
0170     Q[X_] *= corrf * Q[T_] / M;
0171     Q[Y_] *= corrf * Q[T_] / M;
0172     Q[Z_] *= corrf * Q[T_] / M;
0173     M = Q[T_] * corrf;
0174   }
0175
0176   v = 0.5 * (vl + vh);
0177   e = Q[T_] - M * v;
0178   if (e < 0.) e = 0.;
0179   nb = Q[NB_] * sqrt(1 - v * v);
0180   nq = Q[NQ_] * sqrt(1 - v * v);
0181   ns = Q[NS_] * sqrt(1 - v * v);
0182   p = eos->p(e, nb, nq, ns) + Pi;
0183   f = (Q[T_] + p) * v - M;
0184   df = (Q[T_] + p) - M * v * dpe;
0185   dvold = vh - vl;
0186   dv = dvold;
0187   for (int i = 0; i < MAXIT; i++) {
0188     if ((f + df * (vh - v)) * (f + df * (vl - v)) > 0. ||
0189         fabs(2. * f) > fabs(dvold * df)) {  // bisection
0190       dvold = dv;
0191       dv = 0.5 * (vh - vl);
0192       v = vl + dv;
0193       //            cout << "BISECTION v = " << setw(12) << v << endl ;
0194     } else {  // Newton
0195       dvold = dv;
0196       dv = f / df;
0197       v -= dv;
0198       //            cout << "NEWTON v = " << setw(12) << v << endl ;
0199     }
0200     if (fabs(dv) < 0.00001) break;
0201
0202     e = Q[T_] - M * v;
0203     if (e < 0.) e = 0.;
0204     nb = Q[NB_] * sqrt(1 - v * v);
0205     nq = Q[NQ_] * sqrt(1 - v * v);
0206     ns = Q[NS_] * sqrt(1 - v * v);
0207     p = eos->p(e, nb, nq, ns) + Pi;
0208     f = (Q[T_] + p) * v - M;
0209     df = (Q[T_] + p) - M * v * dpe;
0210
0211     if (f > 0.)
0212       vh = v;
0213     else
0214       vl = v;
0215     if (debugRiemann)
0216       cout << "step " << i << "  " << e << "  " << nb << "  " << nq << "  "
0217            << ns << "  " << p << "  " << vx << "  " << vy << "  " << vz << endl;
0218     //      if(i==40) { cout << "error : does not converge\n" ; exit(1) ; } ;
0219   }  // for loop
0220      //----------after
0221      // v = 0.5*(vh+vl) ;
0222   vx = v * Q[X_] / M;
0223   vy = v * Q[Y_] / M;
0224   vz = v * Q[Z_] / M;
0225   e = Q[T_] - M * v;
0226   p = eos->p(e, nb, nq, ns);
0227   nb = Q[NB_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
0228   nq = Q[NQ_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
0229   ns = Q[NS_] * sqrt(1 - vx * vx - vy * vy - vz * vz);
0230   if (e < 0. || sqrt(vx * vx + vy * vy + vz * vz) > 1.) {
0231     cout << Q[T_] << "  " << Q[X_] << "  " << Q[Y_] << "  " << Q[Z_] << "  "
0232          << Q[NB_] << endl;
0233     cout << "transformRF::Error\n";
0234   }
0235   if (!(nb < 0. || nb >= 0.)) {
0236     cout << "transformRF::Error nb=#ind\n";
0237     cout << "Q [7]: " << Q[0] << " " << Q[1] << " " << Q[2] << " " << Q[3]
0238          << " " << Q[4] << " " << Q[5] << " " << Q[6] << endl;
0239     cout << "e vx vy vz nb nq ns: " << e << " " << vx << " " << vy << " " << vz
0240          << " " << nb << " " << nq << " " << ns << endl;
0241     handler(333);
0242     //      return ;
0243   }
0244 }
0245
0246 void transformCV(double e, double p, double nb, double nq, double ns, double vx,
0247                  double vy, double vz, double Q[]) {
0248   const double gamma2 = 1. / (1 - vx * vx - vy * vy - vz * vz);
0249   Q[T_] = (e + p * (vx * vx + vy * vy + vz * vz)) * gamma2;
0250   Q[X_] = (e + p) * vx * gamma2;
0251   Q[Y_] = (e + p) * vy * gamma2;
0252   Q[Z_] = (e + p) * vz * gamma2;
0253   Q[NB_] = nb * sqrt(gamma2);
0254   Q[NQ_] = nq * sqrt(gamma2);
0255   Q[NS_] = ns * sqrt(gamma2);
0256 }