5 double *cmin,
double *cmax,
int beg,
int end)
19 for (i = beg; i <= end; i++) {
27 sroot = sqrt(cs2[i]*(1.0 - vx*vx - vt2*cs2[i])*(1.0 - vel2));
29 cmax[
i] = (vx*(1.0 - cs2[
i]) + sroot)/(1.0 - vel2*cs2[
i]);
30 cmin[
i] = (vx*(1.0 - cs2[
i]) - sroot)/(1.0 - vel2*cs2[
i]);
38 double **LL,
double **RR)
68 real u,
v, w, rp,rm, v2tan;
69 #if CHECK_EIGENVECTORS == YES
81 EXPAND(u = q[VXn]/g; ,
93 v2tan = EXPAND(0.0, + v*v, + w*w);
94 eta = sqrt(1.0 - u*u - cs2*v2tan);
95 g2_1d = 1.0/(1.0 - u*u);
101 rm = u*(1.0 - cs2) - cs*eta/g;
102 rp = u*(1.0 - cs2) + cs*eta/g;
108 for (kk = 2; kk <
NFLX; kk++) lambda[kk] = u;
117 rp = cs*rp/(q[
RHO]*(u*cs - g*
eta));
118 rm = cs*rm/(q[
RHO]*(u*cs + g*
eta));
123 #if RECONSTRUCT_4VEL == NO
132 EXPAND(RR[VXn][0] = -a; ,
140 EXPAND(RR[VXn][1] = a; ,
147 EXPAND(RR[
RHO][2] = 1.0; ,
177 LL[3][
VXn] = u*v*g2_1d;
179 LL[3][PRS] = v*g2_1d/(g2*q[
RHO]*h);
186 LL[4][
VXn] = u*w*g2_1d;
188 LL[4][PRS] = w*g2_1d/(g2*q[
RHO]*h);
192 #elif RECONSTRUCT_4VEL == YES
202 EXPAND(r2 = -(g + g3*u*u)*a + g3*u*rm*v2tan; ,
203 r3 = v*(-g3*u*a + rm*(g + g3*v2tan)); ,
204 r4 = w*(-g3*u*a + rm*(g + g3*v2tan)); )
207 EXPAND(RR[VXn][0] = r2; ,
214 EXPAND(r2 = (g + g3*u*u)*a + g3*u*rp*v2tan; ,
215 r3 = v*(g3*u*a + rp*(g + g3*v2tan)); ,
216 r4 = w*(g3*u*a + rp*(g + g3*v2tan)); )
219 EXPAND(RR[VXn][1] = r2; ,
232 EXPAND(RR[VXn][3] = g3*u*v; ,
233 RR[
VXt][3] = g + g3*v*
v; ,
234 RR[
VXb][3] = g3*v*w;)
243 RR[
VXb][4] = g + g3*w*w;
254 EXPAND(LL[0][VXn] =-0.5*(1.0 - u*u)/a/g; ,
255 LL[0][
VXt] = 0.5*u*v/a/g; ,
256 LL[0][
VXb] = 0.5*u*w/a/g;)
261 EXPAND(LL[1][VXn] = 0.5*(1.0 - u*u)/a/g; ,
262 LL[1][
VXt] = -0.5*u*v/a/g; ,
263 LL[1][
VXb] = -0.5*u*w/a/g;)
275 EXPAND(LL[3][VXn] = -0.5*v*(rp - rm)/a/g/g2_1d - u*v/g; ,
276 LL[3][
VXt] = 0.5*v*v*(rp - rm)*u/a/g + (1.0 - v*v)/g; ,
277 LL[3][
VXb] = 0.5*v*(rp - rm)*u*w/a/g - v*w/g;)
278 LL[3][ENG] = -0.5*v*(rm + rp)/b;
284 LL[4][
VXn] = -0.5*w*(rp - rm)/a/g/g2_1d - u*w/g;
285 LL[4][
VXt] = 0.5*v*(rp - rm)*u*w/a/g - v*w/g;
286 LL[4][
VXb] = 0.5*w*w*(rp - rm)*u/a/g + (1.0 - w*w)/g;
287 LL[4][ENG] = -0.5*w*(rp + rm)/b;
296 #if CHECK_EIGENVECTORS == YES
298 for (ii = 0; ii <
NFLX; ii++){
299 for (jj = 0; jj <
NFLX; jj++){
301 for (kk = 0; kk <
NFLX; kk++){
302 for (ll = 0; ll <
NFLX; ll++){
303 AA[ii][jj] += RR[ii][kk]*(kk==ll)*lambda[kk]*LL[ll][jj];
307 PRIM_RHS (q, q, cs2, h, Aw0);
308 for (ii = 0; ii <
NFLX; ii++){
310 for (jj = 0; jj <
NFLX; jj++){
311 Aw1[ii] += AA[ii][jj]*q[jj];
316 for (ii = 0; ii <
NFLX; ii++){
317 a += fabs(Aw1[ii] - Aw0[ii]);
321 print (
"! Eigenvectors not consistent in EIGENV%12.6e\n",a);
322 for (ii = 0; ii <
NFLX; ii++){
323 for (jj = 0; jj <
NFLX; jj++){
324 print (
"%12.6e ",AA[ii][jj]);
331 for (ii = 0; ii <
NFLX; ii++){
332 print (
"%12.6e , ",Aw0[ii]);
335 for (ii = 0; ii <
NFLX; ii++){
336 print (
"%12.6e , ",Aw1[ii]);
343 for (ii = 0; ii <
NFLX;ii++){
344 for (jj = 0; jj <
NFLX;jj++){
346 for (kk = 0; kk <
NFLX;kk++){
347 a += LL[kk][ii]*RR[jj][kk];
349 if (ii==jj && fabs(a-1.0)>1.e-8) {
350 print (
"! Eigenvectors not orthogonal! %d %d %12.6e \n",ii,jj,a);
354 if(ii!=jj && fabs(a)>1.e-8) {
355 print (
"! Eigenvectors not orthogonal (2) %d %d %12.6e !\n",ii,jj,a);
382 #if RECONSTRUCT_4VEL == NO
384 w[0] = L[0][
VXn]*v[
VXn] + L[0][PRS]*v[PRS];
385 w[1] = L[1][
VXn]*v[
VXn] + L[1][PRS]*v[PRS];
386 EXPAND( w[2] = v[
RHO] + L[2][PRS]*v[PRS]; ,
387 w[3] = L[3][
VXn]*v[
VXn] + v[
VXt] + L[3][PRS]*v[PRS]; ,
388 w[4] = L[4][
VXn]*v[
VXn] + v[
VXb] + L[4][PRS]*v[PRS];)
392 w[0] = EXPAND( L[0][
VXn]*v[
VXn],
394 + L[0][
VXb]*v[
VXb]) + L[0][PRS]*v[PRS];
395 w[1] = EXPAND( L[1][VXn]*v[VXn],
397 + L[1][VXb]*v[VXb]) + L[1][PRS]*v[PRS];;
399 w[2] = v[
RHO] + L[2][PRS]*v[PRS];
402 w[3] = EXPAND( L[3][VXn]*v[VXn],
404 + L[3][VXb]*v[VXb]) + L[3][PRS]*v[PRS];
409 L[4][
VXb]*v[
VXb] + L[4][PRS]*v[PRS];
void PrimEigenvectors(double *q, double cs2, double h, double *lambda, double **LL, double **RR)
double g_maxMach
The maximum Mach number computed during integration.
int g_dir
Specifies the current sweep or direction of integration.
void print(const char *fmt,...)
void PrimToChar(double **L, double *v, double *w)
#define RECONSTRUCT_4VEL
When set to YES, reconstruct 4-velocity rather than 3-velocity (only for RHD and RMHD physics modules...
void MaxSignalSpeed(double **v, double *cs2, double *cmin, double *cmax, int beg, int end)
#define QUIT_PLUTO(e_code)