PLUTO
eigenv.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Wave-speeds and characteristic decomposition for the MHD equations.
5 
6  This file contains various functions containing Jacobian-related information
7  such as characteristic signal speeds, eigenvalues and eigenvector
8  decomposition for the MHD module.
9 
10  The function MaxSignalSpeed() computes the maximum and minimum
11  characteristic signal velocity for the MHD equations.
12 
13  The function Eigenvalues() computes the 7 characteristic waves.
14 
15  The function PrimEigenvectors() calculates left and right eigenvectors
16  and the corresponding eigenvalues for the \e primitive form the the
17  MHD equations with an adiabatic or isothermal EoS.
18 
19  The function ConsEigenvectors() provides the characteristic decomposition
20  of the convervative MHD equations.
21 
22  The function PrimToChar() compute the matrix-vector multiplcation
23  between the L matrix (containing primitive left eigenvectors)
24  and the vector v. The result containing the characteristic
25  variables is stored in the vector w = L.v
26 
27  \authors A. Mignone (mignone@ph.unito.it)\n
28  P. Tzeferacos (petros.tzeferacos@ph.unito.it)
29  \date April 02, 2015
30 */
31 /* ///////////////////////////////////////////////////////////////////// */
32 #include "pluto.h"
33 
34 /* ********************************************************************* */
35 void MaxSignalSpeed (double **v, double *cs2, double *cmin, double *cmax,
36  double **bgf, int beg, int end)
37 /*!
38  * Compute the maximum and minimum characteristic velocities for the
39  * MHD equation, cmin= v - cf, cmax = v + cf
40  *
41  * \param [in] v 1-D array of primitive variables
42  * \param [in] cs2 1-D array containing the square of the sound speed
43  * \param [out] cmin 1-D array containing the leftmost characteristic
44  * \param [out] cmin 1-D array containing the rightmost characteristic
45  * \param [in] bgf 1-D array containing the background magnetic field
46  * \param [in] beg starting index of computation
47  * \param [in] end final index of computation
48  *
49  *********************************************************************** */
50 {
51  int i;
52  double gpr, Bmag2, Btmag2;
53  double cf;
54  double b1, b2, b3;
55 
56  b1 = b2 = b3 = 0.0;
57  for (i = beg; i <= end; i++) {
58 
59  /* ----------------------------------------------------------
60  the following are equivalent, but round-off may prevent
61  a couple of tests from reproducing the same log files
62  ---------------------------------------------------------- */
63 
64 #if EOS == IDEAL
65  gpr = g_gamma*v[i][PRS];
66 #else
67  gpr = cs2[i]*v[i][RHO];
68 #endif
69 
70  /* -- get total field -- */
71 
72  EXPAND (b1 = v[i][BXn]; ,
73  b2 = v[i][BXt]; ,
74  b3 = v[i][BXb];)
75 
76 #if BACKGROUND_FIELD == YES
77  EXPAND (b1 += bgf[i][BXn]; ,
78  b2 += bgf[i][BXt]; ,
79  b3 += bgf[i][BXb];)
80 #endif
81  Btmag2 = b2*b2 + b3*b3;
82  Bmag2 = b1*b1 + Btmag2;
83 
84  cf = gpr - Bmag2;
85  cf = gpr + Bmag2 + sqrt(cf*cf + 4.0*gpr*Btmag2);
86  cf = sqrt(0.5*cf/v[i][RHO]);
87  cmin[i] = v[i][VXn] - cf;
88  cmax[i] = v[i][VXn] + cf;
89 /*
90  g_maxMach = MAX(fabs (w[ii][VXn] / sqrt(a2)), g_maxMach);
91 */
92  }
93 }
94 
95 /* ********************************************************************* */
96 void Eigenvalues(double **v, double *csound2, double **lambda,
97  int beg, int end)
98 /*!
99  * Compute eigenvalues for the MHD equations
100  *
101  * \param [in] v 1-D array of primitive variables
102  * \param [out] csound2 1-D array containing the square of sound speed
103  * \param [out] lambda 1-D array [i][nv] containing the eigenvalues
104  * \param [in] beg starting index of computation
105  * \param [in] end final index of computation
106  *
107  *********************************************************************** */
108 {
109  int i, k;
110  double *q;
111  double scrh0, scrh1, scrh2, scrh3, scrh4;
112  double u, a, a2, ca2, cf2, cs2;
113  double cs, ca, cf, b2, A2, At2;
114  double tau;
115  double sqrt_rho;
116 
117  for (i = beg; i <= end; i++){
118  q = v[i];
119 
120  u = q[VXn];
121  tau = 1.0/q[RHO];
122  sqrt_rho = sqrt(q[RHO]);
123 
124  a2 = csound2[i];
125 
126  scrh2 = q[BXn]*q[BXn]; /* > 0 */
127  scrh3 = EXPAND(0.0, + q[BXt]*q[BXt], + q[BXb]*q[BXb]); /* > 0 */
128 
129  b2 = scrh2 + scrh3; /* >0 */
130  ca2 = scrh2*tau; /* >0 if tau >0 */
131  A2 = b2*tau; /* >0 '' '' */
132  At2 = scrh3*tau;
133 
134  scrh1 = a2 - A2;
135  scrh0 = sqrt(scrh1*scrh1 + 4.0*a2*At2); /* >0 */
136 
137 /* Now get fast and slow speeds */
138 
139  cf2 = 0.5*(a2 + A2 + scrh0);
140  cs2 = a2*ca2/cf2;
141 
142  cf = sqrt(cf2);
143  cs = sqrt(cs2);
144  ca = sqrt(ca2);
145  a = sqrt(a2);
146 
147  lambda[i][KFASTM] = u - cf;
148  lambda[i][KFASTP] = u + cf;
149 
150 #if EOS == IDEAL
151  lambda[i][KENTRP] = u;
152 #endif
153 
154 #ifndef GLM_MHD
155  #if DIVB_CONTROL == EIGHT_WAVES
156  lambda[i][KDIVB] = u;
157  #else
158  lambda[i][KDIVB] = 0.0;
159  #endif
160 #endif
161 
162 #if COMPONENTS > 1
163  lambda[i][KSLOWM] = u - cs;
164  lambda[i][KSLOWP] = u + cs;
165 #endif
166 
167 #if COMPONENTS == 3
168  lambda[i][KALFVM] = u - ca;
169  lambda[i][KALFVP] = u + ca;
170 #endif
171 
172 #ifdef GLM_MHD
173  lambda[i][KPSI_GLMM] = -glm_ch;
174  lambda[i][KPSI_GLMP] = glm_ch;
175 #endif
176  }
177 }
178 /* ********************************************************************* */
179 void PrimEigenvectors(double *q, double a2, double h, double *lambda,
180  double **LL, double **RR)
181 /*!
182  * Provide left and right eigenvectors and corresponding
183  * eigenvalues for the PRIMITIVE form of the MHD equations
184  * (adiabatic & isothermal cases).
185  *
186  * \param [in] q vector of primitive variables
187  * \param [in] a2 sound speed squared
188  * \param [in] h enthalpy
189  * \param [out] lambda eigenvalues
190  * \param [out] LL left primitive eigenvectors
191  * \param [out] RR right primitive eigenvectors
192  *
193  * \note It is highly recommended that LL and RR be initialized to
194  * zero \em before since only non-zero entries are treated here.
195  *
196  * Wave names and their order are defined as enumeration constants in
197  * mod_defs.h.
198  * Notice that, the characteristic decomposition may differ
199  * depending on the way div.B is treated.
200  *
201  * Advection modes associated with passive scalars are simple cases
202  * for which lambda = u (entropy mode) and l = r = (0, ... , 1, 0, ...).
203  * For this reason they are NOT defined here and must be treated
204  * seperately.
205  *
206  * \b References:
207  * - "Notes on the Eigensystem of Magnetohydrodynamics", \n
208  * P.L. Roe, D.S. Balsara,
209  * SIAM Journal on Applied Mathematics, 56, 57 (1996)
210  *
211  * - "A solution adaptive upwind scheme for ideal MHD", \n
212  * K. G. Powell, P. L. Roe, and T. J. Linde,
213  * Journal of Computational Physics, 154, 284-309, (1999).
214  *
215  * - "A second-order unsplit Godunov scheme for cell-centered MHD:
216  * the CTU-GLM scheme"\n
217  * Mignone \& Tzeferacos, JCP (2010) 229, 2117
218  *
219  * - "ATHENA: A new code for astrophysical MHD",
220  * J. Stone, T. Gardiner,
221  * ApJS, 178, 137 (2008)
222  *
223  *
224  * The derivation of the isothermal eigenvectors follows the
225  * consideration given in roe.c
226  *
227  *********************************************************************** */
228 #define sqrt_1_2 (0.70710678118654752440)
229 {
230  int i, j, k;
231  double scrh0, scrh1, scrh2, scrh3, scrh4;
232  double u, a, ca2, cf2, cs2;
233  double cs, ca, cf, b2, A2, At2;
234  double tau, S;
235  double alpha_f, alpha_s, beta_y, beta_z;
236  double sqrt_rho;
237 
238 #if BACKGROUND_FIELD == YES
239  print1 ("! Background field does not support characteristic limiting\n");
240  QUIT_PLUTO(1);
241 #endif
242 
243  u = q[VXn];
244  tau = 1.0/q[RHO];
245  sqrt_rho = sqrt(q[RHO]);
246 
247  scrh2 = q[BXn]*q[BXn]; /* Bx^2 */
248  scrh3 = EXPAND(0.0, + q[BXt]*q[BXt], + q[BXb]*q[BXb]); /* Bt^2 */
249 
250  b2 = scrh2 + scrh3; /* B^2 = Bx^2 + Bt^2 */
251  ca2 = scrh2*tau; /* Bx^2/rho */
252  A2 = b2*tau; /* B^2/rho */
253  At2 = scrh3*tau; /* Bt^2/rho */
254 
255  scrh1 = a2 - A2;
256  scrh0 = sqrt(scrh1*scrh1 + 4.0*a2*At2); /* sqrt( (g*p/rho-B^2/rho)^2
257  + 4*g*p/rho*Bt^2/rho) */
258 /* -- Now get fast and slow speeds -- */
259 
260  cf2 = 0.5*(a2 + A2 + scrh0);
261  cs2 = a2*ca2/cf2;
262 
263  cf = sqrt(cf2);
264  cs = sqrt(cs2);
265  ca = sqrt(ca2);
266  a = sqrt(a2);
267 
268  if (cf == cs) {
269  alpha_f = 1.0;
270  alpha_s = 0.0;
271  }else{
272  scrh0 = 1.0/scrh0;
273  alpha_f = (a2 - cs2)*scrh0;
274  alpha_s = (cf2 - a2)*scrh0;
275 
276  alpha_f = MAX(0.0, alpha_f);
277  alpha_s = MAX(0.0, alpha_s);
278 
279  alpha_f = sqrt(alpha_f);
280  alpha_s = sqrt(alpha_s);
281  }
282 
283  scrh0 = sqrt(scrh3);
284  if (scrh0 > 1.e-9) {
285  SELECT ( ,
286  beta_y = DSIGN(q[BXt]); ,
287  beta_y = q[BXt] / scrh0;
288  beta_z = q[BXb] / scrh0;)
289  } else {
290  SELECT ( ,
291  beta_y = 1.0; ,
292  beta_z = beta_y = sqrt_1_2;)
293  }
294 
295  S = (q[BXn] >= 0.0 ? 1.0 : -1.0);
296 
297 /* ------------------------------------------------------------
298  define primitive right and left eigenvectors (RR and LL),
299  for all of the 8 waves;
300  left eigenvectors for fast & slow waves can be defined
301  in terms of right eigenvectors (see page 296)
302  ------------------------------------------------------------ */
303 
304  /* -------------------------
305  FAST WAVE, (u - c_f)
306  ------------------------- */
307 
308  k = KFASTM;
309  lambda[k] = u - cf;
310  scrh0 = alpha_s*cs*S;
311  scrh1 = alpha_s*sqrt_rho*a;
312  scrh2 = 0.5 / a2;
313  scrh3 = scrh2*tau;
314 
315  RR[RHO][k] = q[RHO]*alpha_f;
316  EXPAND(RR[VXn][k] = -cf*alpha_f; ,
317  RR[VXt][k] = scrh0*beta_y; ,
318  RR[VXb][k] = scrh0*beta_z;)
319  EXPAND( ; ,
320  RR[BXt][k] = scrh1*beta_y; ,
321  RR[BXb][k] = scrh1*beta_z;)
322 
323 #if HAVE_ENERGY
324 // scrh4 = alpha_f*g_gamma*q[PRS];
325  scrh4 = alpha_f*a2*q[RHO];
326  RR[PRS][k] = scrh4;
327 #endif
328 
329 #if EOS == ISOTHERMAL
330  LL[k][RHO] = 0.5*alpha_f/q[RHO];
331 #endif
332  EXPAND(LL[k][VXn] = RR[VXn][k]*scrh2; ,
333  LL[k][VXt] = RR[VXt][k]*scrh2; ,
334  LL[k][VXb] = RR[VXb][k]*scrh2;)
335  EXPAND( ; ,
336  LL[k][BXt] = RR[BXt][k]*scrh3; ,
337  LL[k][BXb] = RR[BXb][k]*scrh3;)
338 #if HAVE_ENERGY
339  LL[k][PRS] = alpha_f*scrh3;
340 #endif
341 
342  /* -------------------------
343  FAST WAVE, (u + c_f)
344  ------------------------- */
345 
346  k = KFASTP;
347  lambda[k] = u + cf;
348  RR[RHO][k] = RR[RHO][KFASTM];
349  EXPAND(RR[VXn][k] = -RR[VXn][KFASTM]; ,
350  RR[VXt][k] = -RR[VXt][KFASTM]; ,
351  RR[VXb][k] = -RR[VXb][KFASTM];)
352  EXPAND( ; ,
353  RR[BXt][k] = RR[BXt][KFASTM]; ,
354  RR[BXb][k] = RR[BXb][KFASTM];)
355 #if HAVE_ENERGY
356  RR[PRS][k] = RR[PRS][KFASTM];
357 #endif
358 
359 #if EOS == ISOTHERMAL
360  LL[k][RHO] = LL[KFASTM][RHO];
361 #endif
362  EXPAND(LL[k][VXn] = -LL[KFASTM][VXn]; ,
363  LL[k][VXt] = -LL[KFASTM][VXt]; ,
364  LL[k][VXb] = -LL[KFASTM][VXb];)
365  EXPAND( ; ,
366  LL[k][BXt] = LL[KFASTM][BXt]; ,
367  LL[k][BXb] = LL[KFASTM][BXb];)
368 #if HAVE_ENERGY
369  LL[k][PRS] = LL[KFASTM][PRS];
370 #endif
371 
372  /* -------------------------
373  entropy wave, (u) only
374  in ideal MHD
375  ------------------------- */
376 
377 #if HAVE_ENERGY
378  k = KENTRP;
379  lambda[k] = u;
380  RR[RHO][k] = 1.0;
381  LL[k][RHO] = 1.0;
382  LL[k][PRS] = - 1.0/a2;
383 #endif
384 
385  /* -------------------------
386  magnetic flux, (u)
387  ------------------------- */
388 
389 #ifndef GLM_MHD
390  k = KDIVB;
391  #if DIVB_CONTROL == EIGHT_WAVES
392  lambda[k] = u;
393  RR[BXn][k] = 1.0;
394  LL[k][BXn] = 1.0;
395  #else
396  lambda[k] = 0.0;
397  #endif
398 #endif
399 
400 #if COMPONENTS > 1
401 
402  /* -------------------------
403  SLOW WAVE, (u - c_s)
404  ------------------------- */
405 
406  k = KSLOWM;
407  lambda[k] = u - cs;
408  scrh0 = alpha_f*cf*S;
409  scrh1 = alpha_f*sqrt_rho*a;
410 
411  RR[RHO][k] = q[RHO]*alpha_s;
412  EXPAND(RR[VXn][k] = -cs*alpha_s; ,
413  RR[VXt][k] = -scrh0*beta_y; ,
414  RR[VXb][k] = -scrh0*beta_z;)
415  EXPAND( ; ,
416  RR[BXt][k] = -scrh1*beta_y; ,
417  RR[BXb][k] = -scrh1*beta_z;)
418  #if HAVE_ENERGY
419 // scrh4 = alpha_s*g_gamma*q[PRS];
420  scrh4 = alpha_s*a2*q[RHO];
421  RR[PRS][k] = scrh4;
422  #endif
423 
424  #if EOS == ISOTHERMAL
425  LL[k][RHO] = 0.5*alpha_s/q[RHO];
426  #endif
427  EXPAND(LL[k][VXn] = RR[VXn][k]*scrh2; ,
428  LL[k][VXt] = RR[VXt][k]*scrh2; ,
429  LL[k][VXb] = RR[VXb][k]*scrh2;)
430  EXPAND( ; ,
431  LL[k][BXt] = RR[BXt][k]*scrh3; ,
432  LL[k][BXb] = RR[BXb][k]*scrh3;)
433 
434  #if HAVE_ENERGY
435  LL[k][PRS] = alpha_s*scrh3;
436  #endif
437 
438  /* -------------------------
439  SLOW WAVE, (u + c_s)
440  ------------------------- */
441 
442  k = KSLOWP;
443  lambda[k] = u + cs;
444 
445  RR[RHO][k] = RR[RHO][KSLOWM];
446  EXPAND(RR[VXn][k] = -RR[VXn][KSLOWM]; ,
447  RR[VXt][k] = -RR[VXt][KSLOWM]; ,
448  RR[VXb][k] = -RR[VXb][KSLOWM];)
449 
450  EXPAND( ; ,
451  RR[BXt][k] = RR[BXt][KSLOWM]; ,
452  RR[BXb][k] = RR[BXb][KSLOWM];)
453  #if HAVE_ENERGY
454  RR[PRS][k] = scrh4;
455  #endif
456 
457  #if EOS == ISOTHERMAL
458  LL[k][RHO] = LL[KSLOWM][RHO];
459  #endif
460  EXPAND(LL[k][VXn] = -LL[KSLOWM][VXn]; ,
461  LL[k][VXt] = -LL[KSLOWM][VXt]; ,
462  LL[k][VXb] = -LL[KSLOWM][VXb];)
463  EXPAND( ; ,
464  LL[k][BXt] = LL[KSLOWM][BXt]; ,
465  LL[k][BXb] = LL[KSLOWM][BXb];)
466 
467  #if HAVE_ENERGY
468  LL[k][PRS] = LL[KSLOWM][PRS];
469  #endif
470 #endif
471 
472 #if COMPONENTS == 3
473 
474  /* -------------------------
475  Alfven WAVE, (u - c_a)
476  ------------------------- */
477 
478  k = KALFVM;
479  lambda[k] = u - ca;
480  scrh2 = beta_y*sqrt_1_2;
481  scrh3 = beta_z*sqrt_1_2;
482 
483  RR[VXt][k] = -scrh3;
484  RR[VXb][k] = scrh2;
485  RR[BXt][k] = -scrh3*sqrt_rho*S;
486  RR[BXb][k] = scrh2*sqrt_rho*S;
487 
488  LL[k][VXt] = RR[VXt][k];
489  LL[k][VXb] = RR[VXb][k];
490  LL[k][BXt] = RR[BXt][k]*tau;
491  LL[k][BXb] = RR[BXb][k]*tau;
492 
493  /* -------------------------
494  Alfven WAVE, (u + c_a)
495  (-R, -L of the eigenv
496  defined by Gardiner Stone)
497  ------------------------- */
498 
499  k = KALFVP;
500  lambda[k] = u + ca;
501  RR[VXt][k] = RR[VXt][KALFVM];
502  RR[VXb][k] = RR[VXb][KALFVM];
503  RR[BXt][k] = - RR[BXt][KALFVM];
504  RR[BXb][k] = - RR[BXb][KALFVM];
505 
506  LL[k][VXt] = LL[KALFVM][VXt];
507  LL[k][VXb] = LL[KALFVM][VXb];
508  LL[k][BXt] = - LL[KALFVM][BXt];
509  LL[k][BXb] = - LL[KALFVM][BXb];
510 
511 /* -------------------------------------------------------------------
512  HotFix: when B=0 in a 2D plane and only Bz != 0, enforce zero
513  jumps in BXt. This is useful with CharTracing since 2 equal and
514  opposite jumps may accumulate due to round-off.
515  It also perfectly consistent, since no (Bx,By) can be generated.
516  ------------------------------------------------------------------- */
517 
518  #if DIMENSIONS <= 2
519  if (q[BXn] == 0 && q[BXt] == 0.0) for (k = 0; k < NFLX; k++) RR[BXt][k] = 0.0;
520  #endif
521 #endif
522 
523 
524 #ifdef GLM_MHD
525 
526 /* -------------------------
527  GLM wave, -glm_ch
528  ------------------------- */
529 
530  k = KPSI_GLMM;
531  lambda[k] = -glm_ch;
532  RR[BXn][k] = 1.0;
533  RR[PSI_GLM][k] = -glm_ch;
534 
535  LL[k][BXn] = 0.5;
536  LL[k][PSI_GLM] = -0.5/glm_ch;
537 
538 /* -------------------------
539  GLM wave, +glm_ch
540  ------------------------- */
541 
542  k = KPSI_GLMP;
543  lambda[k] = glm_ch;
544  RR[BXn][k] = 1.0;
545  RR[PSI_GLM][k] = glm_ch;
546 
547  LL[k][BXn] = 0.5;
548  LL[k][PSI_GLM] = 0.5/glm_ch;
549 
550 #endif
551 
552 /* ------------------------------------------------------------
553  Verify eigenvectors consistency by
554 
555  1) checking that A = L.Lambda.R, where A is
556  the Jacobian dF/dU
557  2) verify orthonormality, L.R = R.L = I
558  ------------------------------------------------------------ */
559 
560 #if CHECK_EIGENVECTORS == YES
561 {
562  static double **A, **ALR;
563  double dA;
564 
565  if (A == NULL){
566  A = ARRAY_2D(NFLX, NFLX, double);
567  ALR = ARRAY_2D(NFLX, NFLX, double);
568  #if COMPONENTS != 3
569  print ("! PrimEigenvectors: eigenvector check requires 3 components\n");
570  #endif
571  }
572  #if COMPONENTS != 3
573  return;
574  #endif
575 
576  /* --------------------------------------
577  Construct the Jacobian analytically
578  -------------------------------------- */
579 
580  for (i = 0; i < NFLX; i++){
581  for (j = 0; j < NFLX; j++){
582  A[i][j] = ALR[i][j] = 0.0;
583  }}
584 
585  #if HAVE_ENERGY
586 
587  A[RHO][RHO] = q[VXn]; A[RHO][VXn] = q[RHO];
588  A[VXn][VXn] = q[VXn]; A[VXn][BXt] = q[BXt]*tau; A[VXn][BXb] = q[BXb]*tau; A[VXn][PRS] = tau;
589  A[VXt][VXt] = q[VXn]; A[VXt][BXt] = -q[BXn]*tau;
590  A[VXb][VXb] = q[VXn]; A[VXb][BXb] = -q[BXn]*tau;
591  A[BXt][VXn] = q[BXt]; A[BXt][VXt] = -q[BXn]; A[BXt][BXn] = -q[VXt]; A[BXt][BXt] = q[VXn];
592  A[BXb][VXn] = q[BXb]; A[BXb][VXb] = -q[BXn]; A[BXb][BXn] = -q[VXb]; A[BXb][BXb] = q[VXn];
593  A[PRS][VXn] = a2*q[RHO]; A[PRS][PRS] = q[VXn];
594 
595  #elif EOS == ISOTHERMAL
596 
597  A[RHO][RHO] = q[VXn] ; A[RHO][VXn] = q[RHO];
598  A[VXn][VXn] = q[VXn] ; A[VXn][BXt] = q[BXt]*tau; A[VXn][BXb] = q[BXb]*tau;
599  A[VXn][RHO] = tau*a2;
600  A[VXt][VXt] = q[VXn] ; A[VXt][BXt] = -q[BXn]*tau;
601  A[VXb][VXb] = q[VXn] ; A[VXb][BXb] = -q[BXn]*tau;
602  A[BXt][VXn] = q[BXt] ; A[BXt][VXt] = -q[BXn]; A[BXt][BXn] = -q[VXt]; A[BXt][BXt] = q[VXn];
603  A[BXb][VXn] = q[BXb] ; A[BXb][VXb] = -q[BXn]; A[BXb][BXn] = -q[VXb]; A[BXb][BXb] = q[VXn];
604 
605  #endif
606 
607  #ifdef GLM_MHD
608  A[BXn][PSI_GLM] = 1.0;
609  A[PSI_GLM][BXn] = glm_ch*glm_ch;
610  #endif
611 
612  for (i = 0; i < NFLX; i++){
613  for (j = 0; j < NFLX; j++){
614  ALR[i][j] = 0.0;
615  for (k = 0; k < NFLX; k++){
616  ALR[i][j] += RR[i][k]*lambda[k]*LL[k][j];
617  }
618  }}
619 
620  for (i = 0; i < NFLX; i++){
621  for (j = 0; j < NFLX; j++){
622 
623  /* --------------------------------------------------------
624  NOTE: if the standard 7-wave formulation is adopted,
625  the column and the row corresponding to B(normal)
626  do not exist.
627  However, PLUTO uses a primitive matrix with 8
628  entries and the B(normal) column contains two
629  entries (-vy, -vz) which cannot be recovered using
630  a 7x7 system.
631  -------------------------------------------------------- */
632 
633  if (j == BXn) continue;
634  dA = ALR[i][j] - A[i][j];
635  if (fabs(dA) > 1.e-8){
636  print ("! PrimEigenvectors: eigenvectors not consistent\n");
637  print ("! g_dir = %d\n",g_dir);
638  print ("! A[%d][%d] = %16.9e, R.Lambda.L[%d][%d] = %16.9e\n",
639  i,j, A[i][j], i,j,ALR[i][j]);
640  print ("\n\n A = \n"); ShowMatrix(A, NFLX, 1.e-8);
641  print ("\n\n R.Lambda.L = \n"); ShowMatrix(ALR, NFLX, 1.e-8);
642  QUIT_PLUTO(1);
643  }
644  }}
645 
646 /* -- check orthornomality -- */
647 
648  for (i = 0; i < NFLX; i++){
649  for (j = 0; j < NFLX; j++){
650  #if (DIVB_CONTROL == NO) || (DIVB_CONTROL == CONSTRAINED_TRANSPORT)
651  if (i == KDIVB || j == KDIVB) continue;
652  #endif
653  a = 0.0;
654  for (k = 0; k < NFLX; k++) a += LL[i][k]*RR[k][j];
655  if ( (i == j && fabs(a-1.0) > 1.e-8) ||
656  (i != j && fabs(a)>1.e-8) ) {
657  print ("! PrimEigenvectors: Eigenvectors not orthogonal\n");
658  print ("! i,j = %d, %d %12.6e \n",i,j,a);
659  print ("! g_dir: %d\n",g_dir);
660  QUIT_PLUTO(1);
661  }
662  }}
663 }
664 #endif
665 }
666 #undef sqrt_1_2
667 
668 /* ********************************************************************* */
669 void ConsEigenvectors (double *u, double *v, double a2,
670  double **Lc, double **Rc, double *lambda)
671 /*!
672  * Provide conservative eigenvectors for MHD equations.
673  *
674  * \param [in] u array of conservative variables
675  * \param [in] v array of primitive variables
676  * \param [in] a2 square of sound speed
677  * \param [out] Lc left conservative eigenvectors
678  * \param [out] Rc right conservative eigenvectors
679  * \param [out] lambda eigenvalues
680  *
681  * \b Reference
682  * - "High-order conservative finite difference GLM-MHD schemes for
683  * cell-centered MHD"\n
684  * Mignone, Tzeferacos \& Bodo, JCP (2010) 229, 5896
685  * - "A High-order WENO Finite Difference Scheme for the Equations
686  * of Ideal MHD"\n
687  * Jiang,Wu, JCP 150,561 (1999)
688  *
689  * With the following corrections:
690  *
691  * The components (By, Bz) of L_{1,7} page 571 should be
692  * +\sqrt{rho}a\alpha_s instead of -.
693  * Also, since Bx is not considered as a parameter, one must also add a
694  * component in the fast and slow waves.
695  * This can be seen by forming the conservative eigenvectors from the
696  * primitive ones, see the paper from Powell, Roe Linde.
697  *
698  * The Lc_{1,3,5,7}[Bx] components, according to Serna 09
699  * are -(1-g_gamma)*a_{f,s}*Bx and not (1-g_gamma)*a_{f,s}*Bx. Both are
700  * orthonormal though. She is WRONG!
701  * -Petros-
702  *********************************************************************** */
703 {
704  int nv, i,j,k;
705  double rho;
706  double vx, vy, vz;
707  double Bx, By, Bz;
708  double bx, by, bz;
709  double beta_y, beta_z;
710  double one_sqrho;
711  double cf2, cs2, ca2, v2;
712  double cf, cs, ca, a;
713  double alpha_f, alpha_s;
714  double g1, g2, tau, Gf, Ga, Gs;
715  double scrh0, scrh1, S, bt2, b2, Btmag;
716  double one_gmm;
717  double LBX = 0.0; /* setting LBX to 1.0 will include Powell's */
718  /* eight wave. Set it to 0 to recover the */
719  /* standard 7 wave formulation */
720 
721  #if BACKGROUND_FIELD == YES
722  print ("! ConsEigenvectors: Background field does not support\n");
723  print (" characteristic limiting.\n");
724  QUIT_PLUTO(1);
725  #endif
726 
727  #if EOS == PVTE_LAW
728  print1( "! ConsEigenvectors: cannot be used presently with PVTE_LAW EoS\n");
729  QUIT_PLUTO(1);
730  #endif
731 
732 /* --------------------------------------------------------------
733  If eigenvector check is required, we make sure that
734  U = U(V) is exactly converted from V.
735  This is not the case, with the present version of the code,
736  with FD scheme where V = 0.5*(VL + VR) and U = 0.5*(UL + UR)
737  hence U != U(V).
738  --------------------------------------------------------------- */
739 
740  #if CHECK_EIGENVECTORS == YES
741  u[RHO] = v[RHO];
742  u[MX1] = v[RHO]*v[VX1];
743  u[MX2] = v[RHO]*v[VX2];
744  u[MX3] = v[RHO]*v[VX3];
745  u[BX1] = v[BX1];
746  u[BX2] = v[BX2];
747  u[BX3] = v[BX3];
748  #if EOS == IDEAL
749  u[ENG] = v[PRS]/(g_gamma-1.0)
750  + 0.5*v[RHO]*(v[VX1]*v[VX1] + v[VX2]*v[VX2] + v[VX3]*v[VX3])
751  + 0.5*(v[BX1]*v[BX1] + v[BX2]*v[BX2] + v[BX3]*v[BX3]);
752  #endif
753  #endif
754 
755 /* ----------------------------------------------------------
756  compute speed and eigenvalues
757  ---------------------------------------------------------- */
758 
759  #if EOS == BAROTROPIC
760  print ("! ConsEigenvectors: not defined for barotropic MHD\n");
761  QUIT_PLUTO(1);
762  #endif
763 
764  rho = v[RHO];
765 
766  EXPAND (vx = v[VXn]; ,
767  vy = v[VXt]; ,
768  vz = v[VXb];)
769 
770  EXPAND (Bx = v[BXn]; ,
771  By = v[BXt]; ,
772  Bz = v[BXb];)
773 
774  one_sqrho = 1.0/sqrt(rho);
775 
776  S = (Bx >= 0.0 ? 1.0 : -1.0);
777 
778  EXPAND(bx = Bx*one_sqrho; ,
779  by = By*one_sqrho; ,
780  bz = Bz*one_sqrho;)
781 
782  bt2 = EXPAND(0.0 , + by*by, + bz*bz);
783  b2 = bx*bx + bt2;
784  Btmag = sqrt(bt2*rho);
785  v2 = EXPAND(vx*vx , + vy*vy, + vz*vz);
786 
787 /* ------------------------------------------------------------
788  Compute fast and slow magnetosonic speeds.
789 
790  The following expression appearing in the definitions
791  of the fast magnetosonic speed
792 
793  (a^2 - b^2)^2 + 4*a^2*bt^2 = (a^2 + b^2)^2 - 4*a^2*bx^2
794 
795  is always positive and avoids round-off errors.
796  ------------------------------------------------------------ */
797 
798  scrh0 = a2 - b2;
799  ca2 = bx*bx;
800  scrh0 = scrh0*scrh0 + 4.0*bt2*a2;
801  scrh0 = sqrt(scrh0);
802 
803  cf2 = 0.5*(a2 + b2 + scrh0);
804  cs2 = a2*ca2/cf2; /* -- same as 0.5*(a2 + b2 - scrh0) -- */
805 
806  cf = sqrt(cf2);
807  cs = sqrt(cs2);
808  ca = sqrt(ca2);
809  a = sqrt(a2);
810 
811  if (Btmag > 1.e-9) {
812  SELECT( ,
813  beta_y = DSIGN(By); ,
814  beta_y = By/Btmag;
815  beta_z = Bz/Btmag;)
816  } else {
817  SELECT( ,
818  beta_y = 1.0; ,
819  beta_z = beta_y = sqrt(0.5);)
820  }
821 
822  if (cf == cs) {
823  alpha_f = 1.0;
824  alpha_s = 0.0;
825  }else if (a <= cs) {
826  alpha_f = 0.0;
827  alpha_s = 1.0;
828  }else if (cf <= a){
829  alpha_f = 1.0;
830  alpha_s = 0.0;
831  }else{
832  scrh0 = 1.0/(cf2 - cs2);
833  alpha_f = (a2 - cs2)*scrh0;
834  alpha_s = (cf2 - a2)*scrh0;
835  alpha_f = MAX(0.0, alpha_f);
836  alpha_s = MAX(0.0, alpha_s);
837  alpha_f = sqrt(alpha_f);
838  alpha_s = sqrt(alpha_s);
839  }
840 
841 /* --------------------------------------------------------
842  Compute non-zero entries of conservative
843  eigenvectors (Rc, Lc).
844  -------------------------------------------------------- */
845 
846  #if EOS == IDEAL
847  one_gmm = 1.0 - g_gamma;
848  g1 = 0.5*(g_gamma - 1.0);
849  g2 = (g_gamma - 2.0)/(g_gamma - 1.0);
850  tau = (g_gamma - 1.0)/a2;
851  #elif EOS == ISOTHERMAL
852  one_gmm = g1 = tau = 0.0;
853  #endif
854 
855  scrh0 = EXPAND(0.0, + beta_y*vy, + beta_z*vz);
856  Gf = alpha_f*cf*vx - alpha_s*cs*S*scrh0;
857  Ga = EXPAND(0.0, , + S*(beta_z*vy - beta_y*vz));
858  Gs = alpha_s*cs*vx + alpha_f*cf*S*scrh0;
859 
860 /* -----------------------
861  FAST WAVE (u - c_f)
862  ----------------------- */
863 
864  k = KFASTM;
865  lambda[k] = vx - cf;
866 
867  scrh0 = alpha_s*cs*S;
868  scrh1 = alpha_s*a*one_sqrho;
869  Rc[RHO][k] = alpha_f;
870  EXPAND( Rc[MXn][k] = alpha_f*lambda[k]; ,
871  Rc[MXt][k] = alpha_f*vy + scrh0*beta_y; ,
872  Rc[MXb][k] = alpha_f*vz + scrh0*beta_z; )
873  EXPAND( ; ,
874  Rc[BXt][k] = scrh1*beta_y; ,
875  Rc[BXb][k] = scrh1*beta_z; )
876  #if EOS == IDEAL
877  Rc[ENG][k] = alpha_f*(0.5*v2 + cf2 - g2*a2) - Gf;
878  #endif
879 
880  Lc[k][RHO] = (g1*alpha_f*v2 + Gf)*0.5/a2;
881  #if EOS == ISOTHERMAL
882  Lc[k][RHO] += alpha_f*0.5;
883  #endif
884  EXPAND( Lc[k][MXn] = (one_gmm*alpha_f*vx - alpha_f*cf) *0.5/a2; ,
885  Lc[k][MXt] = (one_gmm*alpha_f*vy + scrh0*beta_y)*0.5/a2; ,
886  Lc[k][MXb] = (one_gmm*alpha_f*vz + scrh0*beta_z)*0.5/a2;)
887  EXPAND( Lc[k][BXn] = LBX*one_gmm*alpha_f*Bx*0.5/a2; ,
888  Lc[k][BXt] = (one_gmm*alpha_f*By + scrh1*rho*beta_y)*0.5/a2; ,
889  Lc[k][BXb] = (one_gmm*alpha_f*Bz + scrh1*rho*beta_z)*0.5/a2; )
890  #if EOS == IDEAL
891  Lc[k][ENG] = alpha_f*(g_gamma - 1.0)*0.5/a2;
892  #endif
893 
894  /* -----------------------
895  FAST WAVE (u + c_f)
896  ----------------------- */
897 
898  k = KFASTP;
899  lambda[k] = vx + cf;
900 
901  Rc[RHO][k] = alpha_f;
902  EXPAND( Rc[MXn][k] = alpha_f*lambda[k]; ,
903  Rc[MXt][k] = alpha_f*vy - scrh0*beta_y; ,
904  Rc[MXb][k] = alpha_f*vz - scrh0*beta_z; )
905  EXPAND( ; ,
906  Rc[BXt][k] = scrh1*beta_y; ,
907  Rc[BXb][k] = scrh1*beta_z; )
908  #if EOS == IDEAL
909  Rc[ENG][k] = alpha_f*(0.5*v2 + cf2 - g2*a2) + Gf;
910  #endif
911 
912  Lc[k][RHO] = (g1*alpha_f*v2 - Gf)*0.5/a2;
913  #if EOS == ISOTHERMAL
914  Lc[k][RHO] += alpha_f*0.5;
915  #endif
916  EXPAND( Lc[k][MXn] = (one_gmm*alpha_f*vx + alpha_f*cf) *0.5/a2; ,
917  Lc[k][MXt] = (one_gmm*alpha_f*vy - scrh0*beta_y)*0.5/a2; ,
918  Lc[k][MXb] = (one_gmm*alpha_f*vz - scrh0*beta_z)*0.5/a2;)
919  EXPAND( Lc[k][BXn] = LBX*one_gmm*alpha_f*Bx*0.5/a2; ,
920  Lc[k][BXt] = (one_gmm*alpha_f*By + sqrt(rho)*a*alpha_s*beta_y)*0.5/a2; ,
921  Lc[k][BXb] = (one_gmm*alpha_f*Bz + sqrt(rho)*a*alpha_s*beta_z)*0.5/a2; )
922  #if EOS == IDEAL
923  Lc[k][ENG] = alpha_f*(g_gamma - 1.0)*0.5/a2;
924  #endif
925 
926  /* -----------------------
927  Entropy wave (u)
928  ----------------------- */
929 
930  #if EOS == IDEAL
931  k = KENTRP;
932  lambda[k] = vx;
933 
934  Rc[RHO][k] = 1.0;
935  EXPAND( Rc[MXn][k] = vx; ,
936  Rc[MXt][k] = vy; ,
937  Rc[MXb][k] = vz; )
938  Rc[ENG][k] = 0.5*v2;
939 
940  Lc[k][RHO] = 1.0 - 0.5*tau*v2;
941  EXPAND(Lc[k][VXn] = tau*vx; ,
942  Lc[k][VXt] = tau*vy; ,
943  Lc[k][VXb] = tau*vz;)
944  EXPAND(Lc[k][BXn] = LBX*tau*Bx; ,
945  Lc[k][BXt] = tau*By; ,
946  Lc[k][BXb] = tau*Bz;)
947  Lc[k][ENG] = -tau;
948  #endif
949 
950  /* --------------------------------
951  div.B wave (u)
952  -------------------------------- */
953 
954  #ifndef GLM_MHD
955  k = KDIVB;
956  #if DIVB_CONTROL == EIGHT_WAVES
957  lambda[k] = vx;
958  Rc[BXn][k] = 1.0;
959  Lc[k][BXn] = 1.0;
960  #else
961  lambda[k] = 0.0;
962  Rc[BXn][k] = Lc[k][BXn] = 0.0;
963  #endif
964  #endif
965 
966  #if COMPONENTS > 1
967 
968  /* -----------------------
969  SLOW WAVE (u - c_s)
970  ----------------------- */
971 
972  k = KSLOWM;
973  lambda[k] = vx - cs;
974  scrh0 = alpha_f*cf*S;
975 
976  Rc[RHO][k] = alpha_s;
977  EXPAND( Rc[MXn][k] = alpha_s*lambda[k]; ,
978  Rc[MXt][k] = alpha_s*vy - scrh0*beta_y; ,
979  Rc[MXb][k] = alpha_s*vz - scrh0*beta_z; )
980  EXPAND( ; ,
981  Rc[BXt][k] = - alpha_f*a*beta_y*one_sqrho; ,
982  Rc[BXb][k] = - alpha_f*a*beta_z*one_sqrho; )
983 
984  #if EOS == IDEAL
985  Rc[ENG][k] = alpha_s*(0.5*v2 + cs2 - g2*a2) - Gs;
986  #endif
987 
988  Lc[k][RHO] = (g1*alpha_s*v2 + Gs)*0.5/a2;
989 
990  #if EOS == ISOTHERMAL
991  Lc[k][RHO] += alpha_s*0.5;
992  #endif
993 
994  EXPAND( Lc[k][MXn] = (one_gmm*alpha_s*vx - alpha_s*cs) *0.5/a2; ,
995  Lc[k][MXt] = (one_gmm*alpha_s*vy - scrh0*beta_y)*0.5/a2; ,
996  Lc[k][MXb] = (one_gmm*alpha_s*vz - scrh0*beta_z)*0.5/a2;)
997  EXPAND( Lc[k][BXn] = LBX*one_gmm*alpha_s*Bx*0.5/a2; ,
998  Lc[k][BXt] = (one_gmm*alpha_s*By - sqrt(rho)*a*alpha_f*beta_y)*0.5/a2; ,
999  Lc[k][BXb] = (one_gmm*alpha_s*Bz - sqrt(rho)*a*alpha_f*beta_z)*0.5/a2; )
1000  #if EOS == IDEAL
1001  Lc[k][ENG] = alpha_s*(g_gamma - 1.0)*0.5/a2;
1002  #endif
1003 
1004  /* -----------------------
1005  SLOW WAVE (u + c_s)
1006  ----------------------- */
1007 
1008  k = KSLOWP;
1009  lambda[k] = vx + cs;
1010 
1011  Rc[RHO][k] = alpha_s;
1012  EXPAND( Rc[MXn][k] = alpha_s*lambda[k]; ,
1013  Rc[MXt][k] = alpha_s*vy + scrh0*beta_y; ,
1014  Rc[MXb][k] = alpha_s*vz + scrh0*beta_z; )
1015  EXPAND( ; ,
1016  Rc[BXt][k] = - alpha_f*a*beta_y*one_sqrho; ,
1017  Rc[BXb][k] = - alpha_f*a*beta_z*one_sqrho; )
1018 
1019  #if EOS == IDEAL
1020  Rc[ENG][k] = alpha_s*(0.5*v2 + cs2 - g2*a2) + Gs;
1021  #endif
1022 
1023  Lc[k][RHO] = (g1*alpha_s*v2 - Gs)*0.5/a2;
1024 
1025  #if EOS == ISOTHERMAL
1026  Lc[k][RHO] += alpha_s*0.5;
1027  #endif
1028 
1029  EXPAND( Lc[k][MXn] = (one_gmm*alpha_s*vx + alpha_s*cs) *0.5/a2; ,
1030  Lc[k][MXt] = (one_gmm*alpha_s*vy + scrh0*beta_y)*0.5/a2; ,
1031  Lc[k][MXb] = (one_gmm*alpha_s*vz + scrh0*beta_z)*0.5/a2;)
1032  EXPAND( Lc[k][BXn] = LBX*one_gmm*alpha_s*Bx*0.5/a2; ,
1033  Lc[k][BXt] = (one_gmm*alpha_s*By - sqrt(rho)*a*alpha_f*beta_y)*0.5/a2; ,
1034  Lc[k][BXb] = (one_gmm*alpha_s*Bz - sqrt(rho)*a*alpha_f*beta_z)*0.5/a2; )
1035  #if EOS == IDEAL
1036  Lc[k][ENG] = alpha_s*(g_gamma - 1.0)*0.5/a2;
1037  #endif
1038  #endif
1039 
1040  #if COMPONENTS == 3
1041 
1042  /* ------------------------
1043  Alfven WAVE (u - c_a)
1044  ------------------------ */
1045 
1046  k = KALFVM;
1047  lambda[k] = vx - ca;
1048 
1049  Rc[MXt][k] = - beta_z*S;
1050  Rc[MXb][k] = + beta_y*S;
1051  Rc[BXt][k] = - beta_z*one_sqrho;
1052  Rc[BXb][k] = beta_y*one_sqrho;
1053  #if EOS == IDEAL
1054  Rc[ENG][k] = - Ga;
1055  #endif
1056 
1057  Lc[k][RHO] = 0.5*Ga;
1058  Lc[k][MXt] = - 0.5*beta_z*S;
1059  Lc[k][MXb] = 0.5*beta_y*S;
1060  Lc[k][BXt] = - 0.5*sqrt(rho)*beta_z;
1061  Lc[k][BXb] = 0.5*sqrt(rho)*beta_y;
1062 
1063  /* -----------------------
1064  Alfven WAVE (u + c_a)
1065  ----------------------- */
1066 
1067  k = KALFVP;
1068  lambda[k] = vx + ca;
1069 
1070  Rc[MXt][k] = Rc[MXt][KALFVM];
1071  Rc[MXb][k] = Rc[MXb][KALFVM];
1072  Rc[BXt][k] = - Rc[BXt][KALFVM];
1073  Rc[BXb][k] = - Rc[BXb][KALFVM];
1074  #if EOS == IDEAL
1075  Rc[ENG][k] = Rc[ENG][KALFVM];
1076  #endif
1077 
1078  Lc[k][RHO] = Lc[KALFVM][RHO];
1079  Lc[k][MXt] = Lc[KALFVM][MXt];
1080  Lc[k][MXb] = Lc[KALFVM][MXb];
1081  Lc[k][BXt] = - Lc[KALFVM][BXt];
1082  Lc[k][BXb] = - Lc[KALFVM][BXb];
1083  #if EOS == IDEAL
1084  Lc[k][ENG] = Lc[KALFVM][ENG];
1085  #endif
1086  #endif
1087 
1088  #ifdef GLM_MHD
1089 
1090  /* -------------------------
1091  GLM wave, -glm_ch
1092  ------------------------- */
1093 
1094  k = KPSI_GLMM;
1095  lambda[k] = -glm_ch;
1096  Rc[BXn][k] = 1.0;
1097  Rc[PSI_GLM][k] = -glm_ch;
1098 
1099  Lc[k][BXn] = 0.5;
1100  Lc[k][PSI_GLM] = -0.5/glm_ch;
1101 
1102  /* -------------------------
1103  GLM wave, +glm_ch
1104  ------------------------- */
1105 
1106  k = KPSI_GLMP;
1107  lambda[k] = glm_ch;
1108  Rc[BXn][k] = 1.0;
1109  Rc[PSI_GLM][k] = glm_ch;
1110 
1111  Lc[k][BXn] = 0.5;
1112  Lc[k][PSI_GLM] = 0.5/glm_ch;
1113 
1114  #endif
1115 
1116 /* --------------------------------------------------------------
1117  Verify eigenvectors consistency by
1118 
1119  1) checking that A = L.Lambda.R, where A is
1120  the Jacobian dF/dU
1121  2) verify orthonormality, L.R = R.L = I
1122 
1123  IMPORTANT: condition 1) can be satisfied only if U and
1124  V are consistently defined, that is, U = U(V).
1125  This is why we perform an extra conversion at the
1126  beginning of this function only when this check is enabled.
1127  -------------------------------------------------------------- */
1128 
1129 #if CHECK_EIGENVECTORS == YES
1130 {
1131  static double **A, **ALR;
1132  double dA, vel2, Bmag2, vB;
1133 
1134  if (A == NULL){
1135  A = ARRAY_2D(NFLX, NFLX, double);
1136  ALR = ARRAY_2D(NFLX, NFLX, double);
1137  #if COMPONENTS != 3
1138  print ("! ConsEigenvectors: eigenvector check requires 3 components\n");
1139  #endif
1140  }
1141  #if COMPONENTS != 3
1142  return;
1143  #endif
1144 
1145  /* --------------------------------------
1146  Construct the Jacobian analytically
1147  -------------------------------------- */
1148 
1149  for (i = 0; i < NFLX; i++){
1150  for (j = 0; j < NFLX; j++){
1151  A[i][j] = ALR[i][j] = 0.0;
1152  }}
1153 
1154  vel2 = v[VXn]*v[VXn] + v[VXt]*v[VXt] + v[VXb]*v[VXb];
1155  Bmag2 = v[BXn]*v[BXn] + v[BXt]*v[BXt] + v[BXb]*v[BXb];
1156  vB = v[VXn]*v[BXn] + v[VXt]*v[BXt] + v[VXb]*v[BXb];
1157  #if EOS == IDEAL
1158  A[RHO][MXn] = 1.0;
1159 
1160  A[MXn][RHO] = -v[VXn]*v[VXn] + (g_gamma-1.0)*vel2*0.5;
1161  A[MXn][MXn] = (3.0 - g_gamma)*v[VXn];
1162  A[MXn][MXt] = (1.0 - g_gamma)*v[VXt];
1163  A[MXn][MXb] = (1.0 - g_gamma)*v[VXb];
1164  A[MXn][BXt] = (2.0 - g_gamma)*v[BXt];
1165  A[MXn][BXb] = (2.0 - g_gamma)*v[BXb];
1166  A[MXn][ENG] = g_gamma - 1.0;
1167 
1168  A[MXt][RHO] = - v[VXn]*v[VXt];
1169  A[MXt][MXn] = v[VXt];
1170  A[MXt][MXt] = v[VXn];
1171  A[MXt][BXt] = - v[BXn];
1172 
1173  A[MXb][RHO] = - v[VXn]*v[VXb];
1174  A[MXb][MXn] = v[VXb];
1175  A[MXb][MXb] = v[VXn];
1176  A[MXb][BXb] = - v[BXn];
1177 
1178  #ifdef GLM_MHD
1179  A[BXn][PSI_GLM] = 1.0;
1180  #endif
1181 
1182  A[BXt][RHO] = (v[BXn]*v[VXt] - v[BXt]*v[VXn])/v[RHO];
1183  A[BXt][MXn] = v[BXt]/v[RHO];
1184  A[BXt][MXt] = -v[BXn]/v[RHO];
1185  A[BXt][BXt] = v[VXn];
1186 
1187  A[BXb][RHO] = (v[BXn]*v[VXb] - v[BXb]*v[VXn])/v[RHO];
1188  A[BXb][MXn] = v[BXb]/v[RHO];
1189  A[BXb][MXb] = -v[BXn]/v[RHO];
1190  A[BXb][BXb] = v[VXn];
1191 
1192  A[ENG][RHO] = 0.5*(g_gamma - 1.0)*vel2*v[VXn]
1193  - (u[ENG] + v[PRS] + 0.5*Bmag2)/v[RHO]*v[VXn] + vB*v[BXn]/v[RHO];
1194 
1195  A[ENG][MXn] = (1.0 - g_gamma)*v[VXn]*v[VXn]
1196  + (u[ENG] + v[PRS] + 0.5*Bmag2)/v[RHO] - v[BXn]*v[BXn]/v[RHO];
1197 
1198  A[ENG][MXt] = (1.0 - g_gamma)*v[VXn]*v[VXt] - v[BXt]*v[BXn]/v[RHO];
1199  A[ENG][MXb] = (1.0 - g_gamma)*v[VXn]*v[VXb] - v[BXb]*v[BXn]/v[RHO];
1200 
1201  A[ENG][BXt] = (2.0 - g_gamma)*v[BXt]*v[VXn] - v[VXt]*v[BXn];
1202  A[ENG][BXb] = (2.0 - g_gamma)*v[BXb]*v[VXn] - v[VXb]*v[BXn];
1203  A[ENG][ENG] = g_gamma*v[VXn];
1204  #ifdef PSI_GLM
1205  A[PSI_GLM][BXn] = glm_ch;
1206  #endif
1207 
1208  #elif EOS == ISOTHERMAL
1209  A[RHO][MXn] = 1.0;
1210 
1211  A[MXn][RHO] = -v[VXn]*v[VXn] + a2;
1212  A[MXn][MXn] = 2.0*v[VXn];
1213  A[MXn][BXt] = v[BXt];
1214  A[MXn][BXb] = v[BXb];
1215 
1216  A[MXt][RHO] = - v[VXn]*v[VXt];
1217  A[MXt][MXn] = v[VXt];
1218  A[MXt][MXt] = v[VXn];
1219  A[MXt][BXt] = - v[BXn];
1220 
1221  A[MXb][RHO] = - v[VXn]*v[VXb];
1222  A[MXb][MXn] = v[VXb];
1223  A[MXb][MXb] = v[VXn];
1224  A[MXb][BXb] = - v[BXn];
1225 
1226  #ifdef GLM_MHD
1227  A[BXn][PSI_GLM] = 1.0;
1228  #endif
1229 
1230  A[BXt][RHO] = (v[BXn]*v[VXt] - v[BXt]*v[VXn])/v[RHO];
1231  A[BXt][MXn] = v[BXt]/v[RHO];
1232  A[BXt][MXt] = -v[BXn]/v[RHO];
1233  A[BXt][BXt] = v[VXn];
1234 
1235  A[BXb][RHO] = (v[BXn]*v[VXb] - v[BXb]*v[VXn])/v[RHO];
1236  A[BXb][MXn] = v[BXb]/v[RHO];
1237  A[BXb][MXb] = -v[BXn]/v[RHO];
1238  A[BXb][BXb] = v[VXn];
1239 
1240  #ifdef PSI_GLM
1241  A[PSI_GLM][BXn] = glm_ch;
1242  #endif
1243 
1244  #endif
1245 
1246  for (i = 0; i < NFLX; i++){
1247  for (j = 0; j < NFLX; j++){
1248  ALR[i][j] = 0.0;
1249  for (k = 0; k < NFLX; k++){
1250  ALR[i][j] += Rc[i][k]*lambda[k]*Lc[k][j];
1251  }
1252  }}
1253 
1254  for (i = 0; i < NFLX; i++){
1255  for (j = 0; j < NFLX; j++){
1256  if (j == BXn) continue;
1257  if (fabs(ALR[i][j] - A[i][j]) > 1.e-6){
1258  print ("! ConsEigenvectors: eigenvectors not consistent\n");
1259  print ("! g_dir = %d\n",g_dir);
1260  print ("! A[%d][%d] = %16.9e, R.Lambda.L[%d][%d] = %16.9e\n",
1261  i,j, A[i][j], i,j,ALR[i][j]);
1262  print ("\n\n A = \n"); ShowMatrix(A, NFLX, 1.e-8);
1263  print ("\n\n R.Lambda.L = \n"); ShowMatrix(ALR, NFLX, 1.e-8);
1264  QUIT_PLUTO(1);
1265  }
1266  }}
1267 
1268 /* -- check orthornomality -- */
1269 
1270  for (i = 0; i < NFLX; i++){
1271  for (j = 0; j < NFLX; j++){
1272  #if (DIVB_CONTROL == NO) || (DIVB_CONTROL == CONSTRAINED_TRANSPORT)
1273  if (i == KDIVB || j == KDIVB) continue;
1274  #endif
1275  a = 0.0;
1276  for (k = 0; k < NFLX; k++) a += Lc[i][k]*Rc[k][j];
1277  if ( (i == j && fabs(a-1.0) > 1.e-8) ||
1278  (i != j && fabs(a)>1.e-8) ) {
1279  print ("! ConsEigenvectors: Eigenvectors not orthogonal\n");
1280  print ("! i,j = %d, %d %12.6e \n",i,j,a);
1281  print ("! g_dir: %d\n",g_dir);
1282  QUIT_PLUTO(1);
1283  }
1284  }}
1285 }
1286 #endif
1287 }
1288 
1289 /* ********************************************************************* */
1290 void PrimToChar (double **Lp, double *v, double *w)
1291 /*!
1292  * Compute the matrix-vector multiplcation between the
1293  * the L matrix (containing primitive left eigenvectors)
1294  * and the vector v. The result containing the characteristic
1295  * variables is stored in the vector w = L.v
1296  *
1297  * For efficiency purpose, multiplication is done
1298  * explicitly, so that only nonzero entries
1299  * of the left primitive eigenvectors are considered.
1300  *
1301  * \param [in] Lp Left eigenvectors
1302  * \param [in] v (difference of) primitive variables
1303  * \param [out] w (difference of) characteristic variables
1304  *
1305  *********************************************************************** */
1306 {
1307  int k;
1308  double wv, wB, *L;
1309 
1310 /* -- fast waves -- */
1311 
1312  L = Lp[KFASTM];
1313  wv = EXPAND(L[VXn]*v[VXn], + L[VXt]*v[VXt], + L[VXb]*v[VXb]);
1314  #if HAVE_ENERGY
1315  wB = EXPAND(L[PRS]*v[PRS], + L[BXt]*v[BXt], + L[BXb]*v[BXb]);
1316  #elif EOS == ISOTHERMAL
1317  wB = EXPAND(L[RHO]*v[RHO], + L[BXt]*v[BXt], + L[BXb]*v[BXb]);
1318  #endif
1319  w[KFASTM] = wv + wB;
1320  w[KFASTP] = -wv + wB;
1321 
1322 /* -- entropy -- */
1323 
1324  #if HAVE_ENERGY
1325  L = Lp[KENTRP];
1326  w[KENTRP] = L[RHO]*v[RHO] + L[PRS]*v[PRS];
1327  #endif
1328 
1329  #ifndef GLM_MHD
1330  L = Lp[KDIVB];
1331  #if DIVB_CONTROL == EIGHT_WAVES
1332  w[KDIVB] = v[BXn];
1333  #else
1334  w[KDIVB] = 0.0;
1335  #endif
1336  #endif
1337 
1338  #if COMPONENTS > 1
1339  L = Lp[KSLOWM];
1340  wv = EXPAND(L[VXn]*v[VXn], + L[VXt]*v[VXt], + L[VXb]*v[VXb]);
1341 
1342  #if HAVE_ENERGY
1343  wB = EXPAND(L[PRS]*v[PRS], + L[BXt]*v[BXt], + L[BXb]*v[BXb]);
1344  #elif EOS == ISOTHERMAL
1345  wB = EXPAND(L[RHO]*v[RHO], + L[BXt]*v[BXt], + L[BXb]*v[BXb]);
1346  #endif
1347 
1348  w[KSLOWM] = wv + wB;
1349  w[KSLOWP] = -wv + wB;
1350 
1351  #if COMPONENTS == 3
1352  L = Lp[KALFVM];
1353  wv = L[VXt]*v[VXt] + L[VXb]*v[VXb];
1354  wB = L[BXt]*v[BXt] + L[BXb]*v[BXb];
1355  w[KALFVM] = wv + wB;
1356  w[KALFVP] = wv - wB;
1357  #endif
1358  #endif
1359 
1360  #ifdef GLM_MHD
1361  L = Lp[KPSI_GLMP];
1362  w[KPSI_GLMM] = L[BXn]*v[BXn] - L[PSI_GLM]*v[PSI_GLM];
1363  w[KPSI_GLMP] = L[BXn]*v[BXn] + L[PSI_GLM]*v[PSI_GLM];
1364  #endif
1365 
1366 /* -----------------------------------------------
1367  For passive scalars, the characteristic
1368  variable is equal to the primitive one,
1369  since l = r = (0,..., 1 , 0 ,....)
1370  ----------------------------------------------- */
1371 
1372 #if NSCL > 0
1373  NSCL_LOOP(k) w[k] = v[k];
1374 #endif
1375 
1376 /* -------------------------------------------------
1377  verify that the previous simplified expressions
1378  are indeed w = Lp.v
1379  ------------------------------------------------- */
1380 
1381 #if CHECK_EIGENVECTORS == YES
1382 {
1383  int i;
1384  double w2[NVAR];
1385 
1386  for (k = 0; k < NVAR; k++){
1387  w2[k] = 0.0;
1388  for (i = 0; i < NVAR; i++) w2[k] += Lp[k][i]*v[i];
1389  if (fabs(w[k] - w2[k]) > 1.e-8){
1390  printf ("! PrimToChar: projection not correct, k = %d\n",k);
1391  QUIT_PLUTO(1);
1392  }
1393  }
1394 }
1395 #endif
1396 }
IDEAL
#define IDEAL
Definition: pluto.h:45
MX3
#define MX3
Definition: mod_defs.h:22
EOS
#define EOS
Definition: pluto.h:341
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
g_gamma
double g_gamma
Definition: globals.h:112
MX1
#define MX1
Definition: mod_defs.h:20
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
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
PrimEigenvectors
void PrimEigenvectors(double *q, double cs2, double h, double *lambda, double **LL, double **RR)
Definition: eigenv.c:91
ShowMatrix
void ShowMatrix(double **, int n, double)
Definition: tools.c:182
RHO
#define RHO
Definition: mod_defs.h:19
YES
#define YES
Definition: pluto.h:25
ca
static double ca
Definition: init.c:75
sqrt_1_2
#define sqrt_1_2
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
func_param::v
double v[NVAR]
Definition: eos.h:106
DSIGN
#define DSIGN(x)
Definition: macros.h:110
VX1
#define VX1
Definition: mod_defs.h:28
VXb
int VXb
Definition: globals.h:73
KPSI_GLMP
Definition: mod_defs.h:75
BXn
int BXn
Definition: globals.h:75
MXt
int MXt
Definition: globals.h:74
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
Bx
static double Bx
Definition: hlld.c:62
j
int j
Definition: analysis.c:2
MX2
#define MX2
Definition: mod_defs.h:21
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
MXn
int MXn
Definition: globals.h:74
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
BX3
#define BX3
Definition: mod_defs.h:27
pluto.h
PLUTO main header file.
i
int i
Definition: analysis.c:2
ConsEigenvectors
void ConsEigenvectors(double *u, double *v, double a2, double **LL, double **RR, double *lambda)
Definition: eigenv.c:279
BACKGROUND_FIELD
void BACKGROUND_FIELD(real x1, real x2, real x3, real *B0)
Definition: init.c:95
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
glm_ch
double glm_ch
The propagation speed of divergence error.
Definition: glm.c:21
BXt
int BXt
Definition: globals.h:75
KENTRP
Definition: mod_defs.h:59
NVAR
#define NVAR
Definition: pluto.h:609
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
MXb
int MXb
Definition: globals.h:74
KPSI_GLMM
Definition: mod_defs.h:75
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395