PLUTO
eigenv.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Compute the eigenvalues for the relativisitc MHD equations.
5 
6  \author A. Mignone (mignone@ph.unito.it)
7  \date Sep 10, 2012
8 */
9 /* ///////////////////////////////////////////////////////////////////// */
10 #include "pluto.h"
11 
12 /* ********************************************************************* */
13 int MaxSignalSpeed (double **v, double *a2, double *h,
14  double *cmin, double *cmax, int beg, int end)
15 /*!
16  * Return the rightmost (cmax) and leftmost (cmin)
17  * wave propagation speed in the Riemann fan
18  *
19  *********************************************************************** */
20 {
21  int i, err;
22  double lambda[NFLX];
23 
24  for (i = beg; i <= end; i++) {
25  err = Magnetosonic (v[i], a2[i], h[i], lambda);
26  if (err != 0) return i;
27  cmax[i] = lambda[KFASTP];
28  cmin[i] = lambda[KFASTM];
29  }
30  return 0;
31 }
32 
33 /* ********************************************************************* */
34 int Magnetosonic (double *vp, double cs2, double h, double *lambda)
35 /*!
36  * Find the four magneto-sonic speeds (fast and slow) for the
37  * relativistic MHD equations (RMHD).
38  * We follow the notations introduced in Eq. (16) in
39  *
40  * Del Zanna, Bucciantini and Londrillo, A&A, 400, 397--413, 2003
41  *
42  * and also Eq. (57-58) of Mignone & Bodo, MNRAS, 2006 for the
43  * degenerate cases.
44  *
45  * The quartic equation is solved analytically and the solver
46  * (function 'QuarticSolve') assumes all roots are real
47  * (which should be always the case here, since eigenvalues
48  * are recovered from primitive variables). The coefficients
49  * of the quartic were found through the following MAPLE
50  * script:
51  * \verbatim
52  * ------------------------------------------------
53 restart;
54 
55 u[0] := gamma;
56 u[x] := gamma*v[x];
57 b[0] := gamma*vB;
58 b[x] := B[x]/gamma + b[0]*v[x];
59 wt := w + b^2;
60 
61 #############################################################
62 "fdZ will be equation (16) of Del Zanna 2003 times w_{tot}";
63 
64 e2 := c[s]^2 + b^2*(1 - c[s]^2)/wt;
65 fdZ := (1-e2)*(u[0]*lambda - u[x])^4 + (1-lambda^2)*
66  (c[s]^2*(b[0]*lambda - b[x])^2/wt - e2*(u[0]*lambda - u[x])^2);
67 
68 ########################################################
69 "print the coefficients of the quartic polynomial fdZ";
70 
71 coeff(fMB,lambda,4);
72 coeff(fMB,lambda,3);
73 coeff(fMB,lambda,2);
74 coeff(fMB,lambda,1);
75 coeff(fMB,lambda,0);
76 
77 fdZ := fdZ*wt;
78 
79 ######################################################
80 "fMB will be equation (56) of Mignone & Bodo (2006)";
81 
82 a := gamma*(lambda-v[x]);
83 BB := b[x] - lambda*b[0];
84 fMB := w*(1-c[s]^2)*a^4 - (1-lambda^2)*((b^2+w*c[s]^2)*a^2-c[s]^2*BB^2);
85 
86 ######################################
87 "check that fdZ = fMB";
88 
89 simplify(fdZ - fMB);
90 
91 ########################################################
92 "print the coefficients of the quartic polynomial fMB";
93 
94 coeff(fMB,lambda,4);
95 coeff(fMB,lambda,3);
96 coeff(fMB,lambda,2);
97 coeff(fMB,lambda,1);
98 coeff(fMB,lambda,0);
99 
100 
101 ###############################################
102 "Degenerate case 2: Bx = 0, Eq (58) in MB06";
103 
104 B[x] := 0;
105 
106 a2 := w*(c[s]^2 + gamma^2*(1-c[s]^2)) + Q;
107 a1 := -2*w*gamma^2*v[x]*(1-c[s]^2);
108 a0 := w*(-c[s]^2 + gamma^2*v[x]^2*(1-c[s]^2))-Q;
109 
110 "delc is a good way to evaluate the determinant so that it is > 0";
111 del := a1^2 - 4*a2*a0;
112 delc := 4*(w*c[s]^2 + Q)*(w*(1-c[s]^2)*gamma^2*(1-v[x]^2) + (w*c[s]^2+Q));
113 simplify(del-delc);
114 
115 ############################################################################
116 "check the correctness of the coefficients of Eq. (58) in MB06";
117 
118 Q := b^2 - c[s]^2*vB^2;
119 divide(fMB, (lambda-v[x])^2,'fMB2');
120 simplify(a2 - coeff(fMB2,lambda,2)/gamma^2);
121 simplify(a1 - coeff(fMB2,lambda,1)/gamma^2);
122 simplify(a0 - coeff(fMB2,lambda,0)/gamma^2);
123 
124 \endverbatim
125  *
126  *********************************************************************** */
127 {
128  int iflag;
129  double vB, vB2, Bm2, vm2, w, w_1;
130  double eps2, one_m_eps2, a2_w;
131  double vx, vx2, u0, u02;
132  double a4, a3, a2, a1, a0;
133  double scrh1, scrh2, scrh;
134  double b2, del, z[4];
135 
136 #if RMHD_FAST_EIGENVALUES
137  Eigenvalues (vp, cs2, h, lambda);
138  return 0;
139 #endif
140 
141  scrh = fabs(vp[VXn])/sqrt(cs2);
142  g_maxMach = MAX(scrh, g_maxMach);
143 
144  vB = EXPAND(vp[VX1]*vp[BX1], + vp[VX2]*vp[BX2], + vp[VX3]*vp[BX3]);
145  u02 = EXPAND(vp[VX1]*vp[VX1], + vp[VX2]*vp[VX2], + vp[VX3]*vp[VX3]);
146  Bm2 = EXPAND(vp[BX1]*vp[BX1], + vp[BX2]*vp[BX2], + vp[BX3]*vp[BX3]);
147  vm2 = u02;
148 
149  if (u02 >= 1.0){
150  print ("! MAGNETOSONIC: |v|= %f > 1\n",u02);
151  return 1;
152  }
153 
154  if (u02 < 1.e-14) {
155 
156  /* -----------------------------------------------------
157  if total velocity is = 0, the eigenvalue equation
158  reduces to a biquadratic:
159 
160  x^4 + a1*x^2 + a0 = 0
161 
162  with
163  a0 = cs^2*bx*bx/W > 0,
164  a1 = - a0 - eps^2 < 0.
165  ----------------------------------------------------- */
166 
167  w_1 = 1.0/(vp[RHO]*h + Bm2);
168  eps2 = cs2 + Bm2*w_1*(1.0 - cs2);
169  a0 = cs2*vp[BXn]*vp[BXn]*w_1;
170  a1 = - a0 - eps2;
171  del = a1*a1 - 4.0*a0;
172  del = MAX(del, 0.0);
173  del = sqrt(del);
174  lambda[KFASTP] = sqrt(0.5*(-a1 + del));
175  lambda[KFASTM] = -lambda[KFASTP];
176  #if COMPONENTS > 1
177  lambda[KSLOWP] = sqrt(0.5*(-a1 - del));
178  lambda[KSLOWM] = -lambda[KSLOWP];
179  #endif
180  return 0;
181  }
182 
183  vB2 = vB*vB;
184  u02 = 1.0/(1.0 - u02);
185  b2 = Bm2/u02 + vB2;
186  u0 = sqrt(u02);
187  w_1 = 1.0/(vp[RHO]*h + b2);
188  vx = vp[VXn];
189  vx2 = vx*vx;
190 
191  if (fabs(vp[BXn]) < 1.e-14){
192 
193  w = vp[RHO]*h;
194  scrh1 = b2 + cs2*(w - vB2); /* wc_s^2 + Q */
195  scrh2 = w*(1.0 - cs2)*u02;
196  del = 4.0*scrh1*(scrh2*(1.0 - vx2) + scrh1);
197 
198  a2 = scrh1 + scrh2;
199  a1 = -2.0*vx*scrh2;
200 
201  lambda[KFASTP] = 0.5*(-a1 + sqrt(del))/a2;
202  lambda[KFASTM] = 0.5*(-a1 - sqrt(del))/a2;
203  #if COMPONENTS > 1
204  lambda[KSLOWP] = lambda[KSLOWM] = vp[VXn];
205  #endif
206 
207  return 0;
208 
209  }else{
210 
211  scrh1 = vp[BXn]/u02 + vB*vx; /* -- this is bx/u0 -- */
212  scrh2 = scrh1*scrh1;
213 
214  a2_w = cs2*w_1;
215  eps2 = (cs2*vp[RHO]*h + b2)*w_1;
216  one_m_eps2 = u02*vp[RHO]*h*(1.0 - cs2)*w_1;
217 
218  /* ---------------------------------------
219  Define coefficients for the quartic
220  --------------------------------------- */
221 
222  scrh = 2.0*(a2_w*vB*scrh1 - eps2*vx);
223  a4 = one_m_eps2 - a2_w*vB2 + eps2;
224  a3 = - 4.0*vx*one_m_eps2 + scrh;
225  a2 = 6.0*vx2*one_m_eps2 + a2_w*(vB2 - scrh2) + eps2*(vx2 - 1.0);
226  a1 = - 4.0*vx*vx2*one_m_eps2 - scrh;
227  a0 = vx2*vx2*one_m_eps2 + a2_w*scrh2 - eps2*vx2;
228 
229  if (a4 < 1.e-12){
230  print ("! Magnetosonic: cannot divide by a4\n");
231  return 1;
232  }
233 
234  scrh = 1.0/a4;
235 
236  a3 *= scrh;
237  a2 *= scrh;
238  a1 *= scrh;
239  a0 *= scrh;
240  iflag = QuarticSolve(a3, a2, a1, a0, z);
241 
242  if (iflag){
243  print ("! Magnetosonic: Cannot find max speed [dir = %d]\n", g_dir);
244  print ("! rho = %12.6e\n",vp[RHO]);
245  print ("! prs = %12.6e\n",vp[PRS]);
246  EXPAND(print ("! vx = %12.6e\n",vp[VX1]); ,
247  print ("! vy = %12.6e\n",vp[VX2]); ,
248  print ("! vz = %12.6e\n",vp[VX3]);)
249  EXPAND(print ("! Bx = %12.6e\n",vp[BX1]); ,
250  print ("! By = %12.6e\n",vp[BX2]); ,
251  print ("! Bz = %12.6e\n",vp[BX3]);)
252  print ("! f(x) = %12.6e + x*(%12.6e + x*(%12.6e ",
253  a0*a4, a1*a4, a2*a4);
254  print (" + x*(%12.6e + x*%12.6e)))\n", a3*a4, a4);
255  return 1;
256  }
257 /*
258 if (z[3] < z[2] || z[3] < z[1] || z[3] < z[0] ||
259  z[2] < z[1] || z[2] < z[0] || z[1] < z[0]){
260  printf ("!Sorting not correct\n");
261  printf ("z = %f %f %f %f\n",z[0],z[1],z[2],z[3]);
262  QUIT_PLUTO(1);
263 }
264 */
265 
266  lambda[KFASTM] = z[0];
267  lambda[KFASTP] = z[3];
268  #if COMPONENTS > 1
269  lambda[KSLOWP] = z[2];
270  lambda[KSLOWM] = z[1];
271  #endif
272  }
273  return 0;
274 }
275 
276 /* ********************************************************************* */
277 int Eigenvalues(double *v, double cs2, double h, double *lambda)
278 /*!
279  * Compute an approximate expression for the fast magnetosonic speed
280  * using the upper-bound estimated outlined by
281  * Leismann et al. (A&A 2005, 436, 503), Eq. [27].
282  *
283  *********************************************************************** */
284 {
285  double vel2, lor2, Bmag2, b2, ca2, om2, vB2;
286  double vB, scrh, vl, vx, vx2;
287 
288  vel2 = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
289  vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
290  Bmag2 = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]);
291 
292  if (vel2 >= 1.0){
293  print ("! Eigenavalues(): |v|^2 = %f > 1\n",vel2);
294  return 1;
295  }
296 
297  vx = v[VXn];
298  vx2 = vx*vx;
299  vB2 = vB*vB;
300  b2 = Bmag2*(1.0 - vel2) + vB2;
301  ca2 = b2/(v[RHO]*h + b2);
302 
303  om2 = cs2 + ca2 - cs2*ca2;
304  vl = vx*(1.0 - om2)/(1.0 - vel2*om2);
305  scrh = om2*(1.0 - vel2)*((1.0 - vel2*om2) - vx2*(1.0 - om2));
306  scrh = sqrt(scrh)/(1.0 - vel2*om2);
307  lambda[KFASTM] = vl - scrh;
308  lambda[KFASTP] = vl + scrh;
309 
310  if (fabs(lambda[KFASTM])>1.0 || fabs(lambda[KFASTP]) > 1.0){
311  print ("! Eigenvalues(): vm, vp = %8.3e, %8.3e\n",lambda[KFASTM],lambda[KFASTP]);
312  QUIT_PLUTO(1);
313  }
314  if (lambda[KFASTM] != lambda[KFASTM] || lambda[KFASTP] != lambda[KFASTP]){
315  print ("! Eigenvalues(): nan, vm, vp = %8.3e, %8.3e\n",lambda[KFASTM],lambda[KFASTP]);
316  QUIT_PLUTO(1);
317  }
318 
319  scrh = fabs(vx)/sqrt(cs2);
320  g_maxMach = MAX(scrh, g_maxMach);
321 }
322 
323 
324 
325 
326 
327 #if 1==0 /* this part of the code is beta-testing
328  and is excluded from normal usage */
329 #define DEBUG_MODE NO
330 /* **************************************************************** */
331 void PRIM_EIGENVECTORS(real *q, real cs2, real h, real *lambda,
332  real **LL, real **RR)
333 /*
334  *
335  * PURPOSE
336  *
337  * Provide right and left eigenvectors for the primitive
338  * variables (\rho, v, p, B).
339  * We use the procedure outlined in
340  *
341  * "Total Variation Diminishing Scheme for Relativistic MHD"
342  * Balsara, Apj (2001), 132:83.
343  *
344  * The notations are identical, except for the covariant
345  * magnetic field for which we use b^\mu instead of h^\mu
346  * We use square brackets to describe four-vector, always
347  * intended with an upper subscript,
348  *
349  * v[mu] = v^\mu
350  *
351  * !Beware of the distinction between lower and upper index!
352  *
353  * u^\mu = \gamma*(1,v) u_\mu = \gamma(-1,v)
354  * phi^\mu = (\lambda,1,0,0) phi_\mu = (-\lambda,1,0,0)
355  *
356  *
357  *
358 ##############################################
359 #
360 # Maple script to check the behavior of
361 # the eigenvalue equations
362 # (Useful to check degeneracies)
363 #
364 ##############################################
365 
366 restart;
367 
368 v[x] := 0; v[y] := 0; v[z] := 0;
369 B[x] := BB; B[y] := 0; B[z] := 0;
370 
371 vB := v[x]*B[x] + v[y]*B[y] + v[z]*B[z];
372 v2 := v[x]*v[x] + v[y]*v[y] + v[z]*v[z];
373 
374 u[0] := 1/sqrt(1-v2);
375 u[1] := u[0]*v[x];
376 u[2] := u[0]*v[y];
377 u[3] := u[0]*v[z];
378 
379 b[0] := u[0]*vB;
380 b[1] := B[x]/u[0] + u[x]*vB;
381 b[2] := B[y]/u[0] + u[y]*vB;
382 b[3] := B[z]/u[0] + u[z]*vB;
383 
384 b2 := (B[x]^2 + B[y]^2 + B[z]^2)/u[0]^2 + vB^2;
385 
386 #lambda := b[1]/(u[0]*sqrt(w));
387 e2 := cs^2*(1-b2/w) + b2/w;
388 f := (1-e2)*(u[0]*lambda - u[1])^4 +
389  (1-lambda^2)*(cs^2/w*(b[0]*lambda - b[1])^2 -
390  e2*(u[0]*lambda - u[1])^2);
391 simplify(f);
392 #################################################
393  *
394  * LAST MODIFIED
395  *
396  * 11 Aug 2009 by A. Mignone (mignone@ph.unito.it)
397  *
398  ****************************************************************** */
399 {
400  int nv, mu, k, indx[NFLX];
401  int degen1 = 0, degen2[NFLX];
402  double vB, Bmag2, eta, b2, E;
403  double eta_1, e_P, rho_P, rho_S;
404  double a, A, B, betay, betaz;
405  double u[4], b[4], l[4], f[4], d[4];
406  double la[4], lB[4], fa[4], fB[4];
407  double phi[4];
408  double s, Ru[10], col[NFLX];
409  static double **tmp;
410 
411 #if DEBUG_MODE == YES
412 q[RHO] = 1.7;
413 q[VXn] = 0.;
414 q[VXt] = 0.0;
415 q[VXb] = 0.0;
416 q[BXn] = 0.65723;
417 q[BXt] = 0.0;
418 q[BXb] = 0.0;
419 q[PRS] = 2.;
420 cs2 = g_gamma*q[PRS]/(q[PRS]+g_gamma/(g_gamma-1.0)*q[RHO]);
421 h = 1.0+g_gamma/(g_gamma-1.0)*q[PRS]/q[RHO];
422 
423 printf ("AA = %12.6e\n",q[BXn]/sqrt(q[BXn]*q[BXn] + q[RHO]*h));
424 #endif
425 
426  if (tmp == NULL) tmp = ARRAY_2D(NFLX, NFLX, double);
427 
428  u[0] = EXPAND(q[VX1]*q[VX1], + q[VX2]*q[VX2], + q[VX3]*q[VX3]);
429  u[0] = 1.0/sqrt(1.0 - u[0]);
430  vB = EXPAND(q[VX1]*q[BX1], + q[VX2]*q[BX2], + q[VX3]*q[BX3]);
431  Bmag2 = EXPAND(q[BX1]*q[BX1], + q[BX2]*q[BX2], + q[BX3]*q[BX3]);
432 
433  EXPAND(u[1] = u[0]*q[VXn]; ,
434  u[2] = u[0]*q[VXt]; ,
435  u[3] = u[0]*q[VXb];)
436 
437  b[0] = u[0]*vB;
438  EXPAND(b[1] = q[BXn]/u[0] + b[0]*q[VXn]; ,
439  b[2] = q[BXt]/u[0] + b[0]*q[VXt]; ,
440  b[3] = q[BXb]/u[0] + b[0]*q[VXb];)
441  b2 = Bmag2/(u[0]*u[0]) + vB*vB;
442 
443 /* -- thermodynamics relations -- */
444 
445  eta = q[RHO] + g_gamma/(g_gamma - 1.0)*q[PRS];
446  E = eta + b2; /* -- this is the total enthalpy -- */
447  eta_1 = 1.0/eta;
448 
449  a = q[RHO];
450  s = q[PRS]/pow(a,g_gamma);
451  a = 1.0/g_gamma;
452  e_P = a*q[RHO]/q[PRS] + 1.0/(g_gamma - 1.0);
453  rho_P = a*q[RHO]/q[PRS];
454  rho_S = -a*q[RHO]/s;
455 
456 /* ---------------------------------------
457  Find eigenvalues and store their
458  content into lambda[k]
459  --------------------------------------- */
460 
461  MAGNETOSONIC (q, cs2, h, lambda);
462  lambda[KALFVP] = (b[1] + u[1]*sqrt(E))/(b[0] + u[0]*sqrt(E));
463  lambda[KALFVM] = (b[1] - u[1]*sqrt(E))/(b[0] - u[0]*sqrt(E));
464  lambda[KENTRP] = q[VXn];
465  #ifdef KDIVB
466  #if DIVB_CONTROL != EIGHT_WAVES
467  lambda[KDIVB] = 0.0;
468  #endif
469  #endif
470 
471 /* ---------------------------------------
472  spot degeneracies
473  --------------------------------------- */
474 
475  for (k = 0; k < NFLX; k++) degen2[k] = 0;
476 
477  B = sqrt(Bmag2);
478  if (fabs(q[BXn]) < 1.e-9*B){
479  degen1 = 1;
480  }
481  if (fabs(lambda[KALFVM] - lambda[KFASTM]) < 1.e-6)
482  degen2[KFASTM] = 1;
483 
484  if (fabs(lambda[KALFVP] - lambda[KFASTP]) < 1.e-6)
485  degen2[KFASTP] = 1;
486 
487  if (fabs(lambda[KALFVM] - lambda[KSLOWM]) < 1.e-6)
488  degen2[KSLOWM] = 1;
489 
490  if (fabs(lambda[KALFVP] - lambda[KSLOWP]) < 1.e-6)
491  degen2[KSLOWP] = 1;
492 
493 
494 
495 #if DEBUG_MODE == YES
496 PRINT_STATE(q,lambda,cs2,h);
497 #endif
498 
499 /* ----------------------------------------
500  ~ Magnetosonic waves ~
501  ------------------
502 
503  Here "a", "B" depends on the wave
504  speed. For efficiency, we split
505  l^\mu = l[mu] and f^\mu = f[mu] into
506 
507  l[mu] = phi[mu] + la[mu]*a + B*lB[mu]
508  f[mu] = fa[mu]*a + fB[mu]*B
509  ---------------------------------------- */
510 
511 
512  for (mu = 0; mu <= COMPONENTS; mu++) {
513  la[mu] = (eta - e_P*b2)*u[mu]*eta_1;
514  lB[mu] = b[mu]*eta_1;
515  fa[mu] = b[mu]*e_P*eta_1;
516  fB[mu] = -u[mu]*eta_1;
517  }
518 
519  /* ----------------------------
520  ~ Fast waves ~
521  ---------------------------- */
522 
523  for (k = 0; k < NFLX; k++){
524  if (k != KFASTM && k != KFASTP) continue;
525 
526  a = u[0]*(q[VXn] - lambda[k]);
527  B = -b[0]*lambda[k] + b[1];
528  A = E*a*a - B*B;
529 
530 #if DEBUG_MODE == YES
531 printf ("FAST, a = %12.6e, A = %12.6e, B = %12.6e\n",a,A,B);
532 #endif
533 
534  phi[0] = lambda[k];
535  EXPAND(phi[1] = 1.0;, phi[2] = 0.0;, phi[3] = 0.0;)
536 
537  if (degen2[k]) {
538  for (mu = 0; mu <= COMPONENTS; mu++) {
539  d[mu] = B*b[mu] - b2*(phi[mu] + a*u[mu]);
540  }
541  }else {
542  for (mu = 0; mu <= COMPONENTS; mu++) {
543  l[mu] = phi[mu] + la[mu]*a + lB[mu]*B;
544  f[mu] = fa[mu]*a + fB[mu]*B;
545  d[mu] = a*(B*f[mu] - a*l[mu])
546  + (B*B - e_P*b2*a*a)*(phi[mu] + 2.0*a*u[mu])*eta_1;
547  }
548  }
549 
550  for (mu = 0; mu <= COMPONENTS; mu++) {
551  Ru[ mu] = a*d[mu];
552  Ru[4+mu] = B*d[mu] + a*A*f[mu];
553  }
554  Ru[8] = a*a*A;
555  Ru[9] = 0.0;
556 
557  /* -- make physical vector -- */
558 
559  RR[RHO][k] = rho_P*Ru[8] + rho_S*Ru[9];
560  EXPAND(RR[VXn][k] = (-q[VXn]*Ru[0] + Ru[1])/u[0]; ,
561  RR[VXt][k] = (-q[VXt]*Ru[0] + Ru[2])/u[0]; ,
562  RR[VXb][k] = (-q[VXb]*Ru[0] + Ru[3])/u[0];)
563  RR[PRS][k] = Ru[8];
564 
565  EXPAND( ,
566  RR[BXt][k] = b[2]*Ru[0] - b[0]*Ru[2] - u[2]*Ru[4] + u[0]*Ru[6]; ,
567  RR[BXb][k] = b[3]*Ru[0] - b[0]*Ru[3] - u[3]*Ru[4] + u[0]*Ru[7];)
568  }
569 
570  /* ----------------------------
571  ~ Slow waves ~
572  ---------------------------- */
573 
574 /* -----------------------------------------
575  When Bx -> 0 redefine b^\mu
576  ----------------------------------------- */
577 
578  if (degen1){
579  B = sqrt(Bmag2);
580  if (Bmag2 < 1.e-9*q[PRS]){
581  betay = betaz = 1.0/sqrt(2.0);
582  }else{
583  betay = q[BXt]/B;
584  betaz = q[BXb]/B;
585  }
586 
587  A = q[VXt]*betay + q[VXb]*betaz;
588  b[0] = u[0]*A;
589  b[1] = u[1]*A;
590  b[2] = betay/u[0] + u[2]*A;
591  b[3] = betaz/u[0] + u[3]*A;
592 
593  b2 = (betay*betay + betaz*betaz)/(u[0]*u[0]) + A*A;
594 
595 #if DEBUG_MODE == YES
596 printf ("!!!! Degen 1: Bx -> 0, b^2 = %12.6e\n", b2);
597 #endif
598  }
599 
600  for (k = 0; k < NFLX; k++){
601  if (k != KSLOWP && k != KSLOWM) continue;
602 
603  if (degen1) { /* -- when Bx->0 also a, A, B -> 0 -- */
604 
605  if (k == KSLOWM) {
606  for (mu = 0; mu <= COMPONENTS; mu++) {
607  Ru[mu] = -b[mu];
608  Ru[mu+4] = b[mu];
609  }
610  }else{
611  for (mu = 0; mu <= COMPONENTS; mu++) {
612  Ru[mu] = b[mu];
613  Ru[mu+4] = b[mu];
614  }
615  }
616 
617  Ru[8] = - b2*sqrt(Bmag2);
618  Ru[9] = 0.0;
619  }else{
620 
621  a = u[0]*(q[VXn] - lambda[k]);
622  B = -b[0]*lambda[k] + b[1];
623  A = E*a*a - B*B;
624 
625 #if DEBUG_MODE == YES
626 printf ("SLOW, a = %12.6e, A = %12.6e, B = %12.6e\n",a,A,B);
627 #endif
628  phi[0] = lambda[k];
629  EXPAND(phi[1] = 1.0;, phi[2] = 0.0;, phi[3] = 0.0;)
630 
631  if (degen2[k]) {
632  A = 0.0;
633  for (mu = 0; mu <= COMPONENTS; mu++) {
634  d[mu] = B*b[mu] - b2*(phi[mu] + a*u[mu]);
635  }
636  }else {
637  for (mu = 0; mu <= COMPONENTS; mu++) {
638  l[mu] = phi[mu] + la[mu]*a + lB[mu]*B;
639  f[mu] = fa[mu]*a + fB[mu]*B;
640  d[mu] = a*(B*f[mu] - a*l[mu])
641  + (B*B - e_P*b2*a*a)*(phi[mu] + 2.0*a*u[mu])*eta_1;
642  }
643  }
644 
645  for (mu = 0; mu <= COMPONENTS; mu++) {
646  Ru[ mu] = a*d[mu];
647  Ru[4+mu] = B*d[mu] + a*A*f[mu];
648  }
649  Ru[8] = a*a*A;
650  Ru[9] = 0.0;
651  }
652 
653  /* -- make physical vector -- */
654 
655  RR[RHO][k] = rho_P*Ru[8] + rho_S*Ru[9];
656  EXPAND(RR[VXn][k] = (-q[VXn]*Ru[0] + Ru[1])/u[0]; ,
657  RR[VXt][k] = (-q[VXt]*Ru[0] + Ru[2])/u[0]; ,
658  RR[VXb][k] = (-q[VXb]*Ru[0] + Ru[3])/u[0];)
659  RR[PRS][k] = Ru[8];
660 
661  EXPAND( ,
662  RR[BXt][k] = b[2]*Ru[0] - b[0]*Ru[2] - u[2]*Ru[4] + u[0]*Ru[6]; ,
663  RR[BXb][k] = b[3]*Ru[0] - b[0]*Ru[3] - u[3]*Ru[4] + u[0]*Ru[7];)
664  }
665 
666 /* -----------------------------------------------------
667  ~ Alfven waves ~
668  ------------
669 
670  They are defined through the inner tensor product
671 
672  d^a = eps^{abcd}\phi_b u_c b_d
673 
674  Since \phi_a = (-lambda,1), the only non-vanishing
675  components, for a=0,1,2,3 are:
676 
677  d^0 = eps^{0123} \phi_1 (u_2b_3 - u_3b_2) =
678  = (u^2b^3 - u^3b^2)
679 
680  d^1 = eps^{1023} \phi_0 (u_2b_3 - u_3b_2) =
681  = lambda*(u^2b^3 - u^3b^2)
682 
683  d^2 = eps^{2013} \phi_0 (u_1b_3 - u_3b_1)
684  + eps^{2103} \phi_1 (u_0b_3 - u_3b_0) =
685  = -lambda*(u_1b_3 - u_3b_1) + B3
686 
687  d^3 = eps^{3012} \phi_0 (u_1b_2 - u_2b_1)
688  + eps^{3102} \phi_1 (u_0b_2 - u_2b_0) =
689  = lambda*(u_1b_2 - u_2b_1) - B2
690 
691  Where
692 
693  \eps^{0123} = \eps^{2013} = \eps{3102} = 1
694  \eps^{1023} = \eps^{2103} = \eps{3012} = -1
695 
696  and B^k = b^ku^0 - u^kb^0 is the lab frame field.
697  ----------------------------------------------------- */
698 
699  for (k = 0; k < NFLX; k++){
700  if (k != KALFVP && k != KALFVM) continue;
701 
702  a = u[0]*(q[VXn] - lambda[k]);
703  B = -b[0]*lambda[k] + b[1];
704  phi[0] = lambda[k];
705  EXPAND(phi[1] = 1.0;, phi[2] = 0.0;, phi[3] = 0.0;)
706 
707  Ru[0] = (u[2]*b[3] - u[3]*b[2]);
708  Ru[1] = lambda[k]*(u[2]*b[3] - u[3]*b[2]);
709  Ru[2] = -lambda[k]*(u[1]*b[3] - u[3]*b[1])
710  + (u[0]*b[3] - u[3]*b[0]);
711  Ru[3] = lambda[k]*(u[1]*b[2] - u[2]*b[1])
712  - (u[0]*b[2] - u[2]*b[0]);
713 
714  s = B/a;
715  if (degen1 && k == KALFVP) s = 1.0;
716  if (degen1 && k == KALFVM) s = -1.0;
717 
718  Ru[4] = Ru[0]*s;
719  Ru[5] = Ru[1]*s;
720  Ru[6] = Ru[2]*s;
721  Ru[7] = Ru[3]*s;
722  Ru[8] = Ru[9] = 0.0;
723 
724  /* -- make physical vector -- */
725 
726  RR[RHO][k] = 0.0;
727  EXPAND(RR[VXn][k] = (-q[VXn]*Ru[0] + Ru[1])/u[0]; ,
728  RR[VXt][k] = (-q[VXt]*Ru[0] + Ru[2])/u[0]; ,
729  RR[VXb][k] = (-q[VXb]*Ru[0] + Ru[3])/u[0];)
730  RR[PRS][k] = Ru[8];
731 
732  EXPAND( ,
733  RR[BXt][k] = b[2]*Ru[0] - b[0]*Ru[2] - u[2]*Ru[4] + u[0]*Ru[6]; ,
734  RR[BXb][k] = b[3]*Ru[0] - b[0]*Ru[3] - u[3]*Ru[4] + u[0]*Ru[7];)
735  }
736 
737  /* -------------------------
738  ~ Entropy wave ~
739  ------------
740  ------------------------- */
741 
742  k = KENTRP;
743  RR[RHO][k] = 1.0;
744 
745  /* -------------------------
746  ~ Divergence wave ~
747  ---------------
748  ------------------------- */
749 
750  #ifdef KDIVB
751  #if DIVB_CONTROL != EIGHT_WAVES
752  k = KDIVB;
753  RR[BXn][k] = 1.0;
754  #endif
755  #endif
756 
757 /* ----------------------------------------------------
758  Now find the left eigenvectors by performing
759  LU decomposition
760  ---------------------------------------------------- */
761 
762  for (nv = NFLX; nv--; ){
763  for (k = NFLX; k--; ){
764  tmp[nv][k] = RR[nv][k];
765  }}
766 
767 #if DEBUG_MODE == YES
768 ShowMatrix(RR);
769 #endif
770 
771  if (!LUDecompose(tmp, NFLX, indx, &a)) {
772 /* printf ("! EIGENV: Singular eigenvectors\n"); */
773  for (nv = NFLX; nv--; ){
774  for (k = NFLX; k--; ){
775  RR[nv][k] = LL[k][nv] = (k == nv ? 1.0:0.0);
776  }}
777  #ifdef KDIVB
778  #if DIVB_CONTROL != EIGHT_WAVES
779  LL[KDIVB][BXn] = RR[BXn][KDIVB] = 0.0;
780  #endif
781  #endif
782  return;
783  }
784 
785  for (mu = 0; mu < NFLX; mu++){
786  for (k = 0; k < NFLX; k++) col[k] = 0.0;
787  col[mu] = 1.0;
788  LUBackSubst(tmp, NFLX, indx, col);
789  for (k = 0; k < NFLX; k++) LL[k][mu] = col[k];
790  }
791 
792  #ifdef KDIVB
793  #if DIVB_CONTROL != EIGHT_WAVES
794  LL[KDIVB][BXn] = RR[BXn][KDIVB] = 0.0;
795  #endif
796  #endif
797 
798 #if DEBUG_MODE == YES
799 printf ("-------------------------------\n");
800 printf (" INVERSE MATRIX\n");
801 printf ("-------------------------------\n");
802 ShowMatrix(LL);
803 exit(1);
804 #endif
805 }
806 #undef DEBUG_MODE
807 
808 void PrimToChar (double **Lp, double *v, double *w)
809 {
810  int k, nv;
811 
812  for (k = NFLX; k--; ){
813  w[k] = 0.0;
814  for (nv = NFLX; nv--; ){
815  w[k] += Lp[k][nv]*v[nv];
816  }
817  }
818 }
819 
820 
821 int PRINT_STATE(double *q, double *lambda, double cs2, double hh)
822 {
823 
824 printf ("q[RHO] = %12.6e;\n", q[RHO]);
825 printf ("q[VXn] = %12.6e;\n", q[VXn]);
826 printf ("q[VXt] = %12.6e;\n", q[VXt]);
827 printf ("q[VXb] = %12.6e;\n", q[VXb]);
828 
829 printf ("q[BXn] = %12.6e;\n", q[BXn]);
830 printf ("q[BXt] = %12.6e;\n", q[BXt]);
831 printf ("q[BXb] = %12.6e;\n", q[BXb]);
832 
833 printf ("q[PRS] = %12.6e;\n", q[PRS]);
834 printf ("cs2 = %12.6e; hh = %12.6e;\n",cs2,hh);
835 
836 
837 printf ("fm, %d %12.6e\n", KFASTM, lambda[KFASTM]);
838 printf ("am, %d %12.6e\n", KALFVM, lambda[KALFVM]);
839 printf ("sm, %d %12.6e\n", KSLOWM, lambda[KSLOWM]);
840 printf ("en, %d %12.6e\n", KENTRP, lambda[KENTRP]);
841 printf ("sp, %d %12.6e\n", KSLOWP, lambda[KSLOWP]);
842 printf ("ap, %d %12.6e\n", KALFVP, lambda[KALFVP]);
843 printf ("fp, %d %12.6e\n", KFASTP, lambda[KFASTP]);
844 
845 }
846 
847 #endif /* 1==0 */
848 
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
g_gamma
double g_gamma
Definition: globals.h:112
QuarticSolve
int QuarticSolve(double, double, double, double, double *)
Definition: quartic.c:23
VX2
#define VX2
Definition: mod_defs.h:29
Eigenvalues
void Eigenvalues(double **v, double *csound2, double **lambda, int beg, int end)
Definition: eigenv.c:67
real
double real
Definition: pluto.h:488
PRIM_EIGENVECTORS
void PRIM_EIGENVECTORS(double *, double, double, double *, double **, double **)
ShowMatrix
void ShowMatrix(double **, int n, double)
Definition: tools.c:182
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
eta
static double *** eta[3]
Definition: res_functions.c:94
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
VXb
int VXb
Definition: globals.h:73
BXn
int BXn
Definition: globals.h:75
g_maxMach
double g_maxMach
The maximum Mach number computed during integration.
Definition: globals.h:119
LUBackSubst
void LUBackSubst(double **a, int n, int *indx, double b[])
Definition: math_lu_decomp.c:84
NFLX
#define NFLX
Definition: mod_defs.h:32
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
VXt
int VXt
Definition: globals.h:73
KFASTP
Definition: mod_defs.h:58
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
KFASTM
Definition: mod_defs.h:58
PrimToChar
void PrimToChar(double **L, double *v, double *w)
Definition: eigenv.c:593
VXn
int VXn
Definition: globals.h:73
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
COMPONENTS
#define COMPONENTS
Definition: definitions_01.h:3
pluto.h
PLUTO main header file.
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
MaxSignalSpeed
void MaxSignalSpeed(double **v, double *cs2, double *cmin, double *cmax, int beg, int end)
Definition: eigenv.c:34
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
BX2
#define BX2
Definition: mod_defs.h:26
LUDecompose
int LUDecompose(double **a, int n, int *indx, double *d)
Definition: math_lu_decomp.c:12
BXt
int BXt
Definition: globals.h:75
KENTRP
Definition: mod_defs.h:59
q
static Runtime q
Definition: runtime_setup.c:490
BXb
int BXb
Definition: globals.h:75
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
Magnetosonic
int Magnetosonic(double *vp, double cs2, double h, double *lambda)
Definition: eigenv.c:34