PLUTO
eigenv.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Wave-speeds and characteristic decomposition for the HD 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 HD module.
9 
10  The function MaxSignalSpeed() computes the maximum and minimum
11  characteristic signal velocity for the HD equations.
12 
13  The function Eigenvalues() computes the 3 characteristic wave speed
14 
15  The function PrimEigenvectors() calculates left and right eigenvectors
16  and the corresponding eigenvalues for the \e primitive form the the
17  HD equations with an adiabatic, general convex or isothermal EoS.
18 
19  The function ConsEigenvectors() provides the characteristic decomposition
20  of the convervative HD 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  \author A. Mignone (mignone@ph.unito.it)
28  \date April 02, 2015
29 */
30 /* ///////////////////////////////////////////////////////////////////// */
31 #include "pluto.h"
32 
33 /* ********************************************************************* */
34 void MaxSignalSpeed (double **v, double *cs2, double *cmin, double *cmax,
35  int beg, int end)
36 /*!
37  * Compute the maximum and minimum characteristic velocities for the
38  * HD equation from the vector of primitive variables v.
39  *
40  * \param [in] v 1-D array of primitive variables
41  * \param [in] cs2 1-D array containing the square of the sound speed
42  * \param [out] cmin 1-D array containing the leftmost characteristic
43  * \param [out] cmin 1-D array containing the rightmost characteristic
44  * \param [in] beg starting index of computation
45  * \param [in] end final index of computation
46  *
47  *********************************************************************** */
48 {
49  int i;
50  double a;
51 
52  for (i = beg; i <= end; i++) {
53  #if HAVE_ENERGY
54 /* a = sqrt(g_gamma*v[i][PRS]/v[i][RHO]); */
55  a = sqrt(cs2[i]);
56  #elif EOS == ISOTHERMAL
57 /* a = g_isoSoundSpeed; */
58  a = sqrt(cs2[i]);
59  #endif
60 
61  cmin[i] = v[i][VXn] - a;
62  cmax[i] = v[i][VXn] + a;
63  }
64 }
65 
66 /* ********************************************************************* */
67 void Eigenvalues(double **v, double *csound2, double **lambda, int beg, int end)
68 /*!
69  * Compute eigenvalues for the HD equations
70  *
71  * \param [in] v 1-D array of primitive variables
72  * \param [out] csound2 1-D array containing the square of sound speed
73  * \param [out] lambda 1-D array [i][nv] containing the eigenvalues
74  * \param [in] beg starting index of computation
75  * \param [in] end final index of computation
76  *
77  *********************************************************************** */
78 {
79  int i, k;
80  double cs;
81 
82  for (i = beg; i <= end; i++){
83  cs = sqrt(csound2[i]);
84  lambda[i][0] = v[i][VXn] - cs;
85  lambda[i][1] = v[i][VXn] + cs;
86  for (k = 2; k < NFLX; k++) lambda[i][k] = v[i][VXn];
87  }
88 }
89 
90 /* ********************************************************************* */
91 void PrimEigenvectors (double *q, double cs2, double h, double *lambda,
92  double **LL, double **RR)
93 /*!
94  * Provide left and right eigenvectors and corresponding
95  * eigenvalues for the primitive form of the HD equations
96  * (adiabatic, pvte & isothermal cases).
97  *
98  * \param [in] q vector of primitive variables
99  * \param [in] cs2 sound speed squared
100  * \param [in] h enthalpy
101  * \param [out] lambda eigenvalues
102  * \param [out] LL matrix containing left primitive eigenvectors (rows)
103  * \param [out] RR matrix containing right primitive eigenvectors (columns)
104  *
105  * \note It is highly recommended that LL and RR be initialized to
106  * zero *BEFORE* since only non-zero entries are treated here.
107  *
108  * Wave names and their order are defined as enumeration constants in
109  * mod_defs.h.
110  *
111  * Advection modes associated with passive scalars are simple cases
112  * for which lambda = u (entropy mode) and l = r = (0, ... , 1, 0, ...).
113  * For this reason they are NOT defined here and must be treated
114  * seperately.
115  *
116  * \b References:
117  * - "Riemann Solvers and Numerical Methods for Fluid Dynamics", \n
118  * Toro, 1997 Springer-Verlag, Eq. [3.18]
119  *
120  *********************************************************************** */
121 {
122  int i, j, k;
123  double rhocs, rho_cs, cs;
124 #if CHECK_EIGENVECTORS == YES
125  double Aw1[NFLX], Aw0[NFLX], AA[NFLX][NFLX], a;
126 #endif
127 
128  cs = sqrt(cs2);
129 
130  rhocs = q[RHO]*cs;
131  rho_cs = q[RHO]/cs;
132 
133 /* ------------------------------------------------------
134  1. Compute RIGHT eigenvectors
135  ------------------------------------------------------ */
136 
137  lambda[0] = q[VXn] - cs; /* lambda = u - c */
138  RR[RHO][0] = 0.5*rho_cs;
139  RR[VXn][0] = -0.5;
140 #if HAVE_ENERGY
141  RR[PRS][0] = 0.5*rhocs;
142 #endif
143 
144  lambda[1] = q[VXn] + cs; /* lambda = u + c */
145  RR[RHO][1] = 0.5*rho_cs;
146  RR[VXn][1] = 0.5;
147 #if HAVE_ENERGY
148  RR[PRS][1] = 0.5*rhocs;
149 #endif
150 
151 
152  for (k = 2; k < NFLX; k++) lambda[k] = q[VXn]; /* lambda = u */
153 
154 #if HAVE_ENERGY
155  EXPAND(RR[RHO][2] = 1.0; ,
156  RR[VXt][3] = 1.0; ,
157  RR[VXb][4] = 1.0;)
158 #elif EOS == ISOTHERMAL
159  EXPAND( ,
160  RR[VXt][2] = 1.0; ,
161  RR[VXb][3] = 1.0;)
162 #endif
163 
164 /* ------------------------------------------------------
165  2. Compute LEFT eigenvectors
166  ------------------------------------------------------ */
167 
168  LL[0][VXn] = -1.0;
169 #if HAVE_ENERGY
170  LL[0][PRS] = 1.0/rhocs;
171 #elif EOS == ISOTHERMAL
172  LL[0][RHO] = 1.0/rhocs;
173 #endif
174 
175  LL[1][VXn] = 1.0;
176 #if HAVE_ENERGY
177  LL[1][PRS] = 1.0/rhocs;
178 #elif EOS == ISOTHERMAL
179  LL[1][RHO] = 1.0/rhocs;
180 #endif
181 
182 #if HAVE_ENERGY
183  EXPAND(LL[2][RHO] = 1.0; ,
184  LL[3][VXt] = 1.0; ,
185  LL[4][VXb] = 1.0;)
186 
187  LL[2][PRS] = -1.0/cs2;
188 #elif EOS == ISOTHERMAL
189  EXPAND( ,
190  LL[2][VXt] = 1.0; ,
191  LL[3][VXb] = 1.0;)
192 #endif
193 
194 /* ------------------------------------------------------------
195  3. If required, verify eigenvectors consistency by
196 
197  1) checking that A = L.Lambda.R, where A is
198  the Jacobian dF/dU
199  2) verify orthonormality, L.R = R.L = I
200  ------------------------------------------------------------ */
201 
202 #if CHECK_EIGENVECTORS == YES
203 {
204  static double **A, **ALR;
205  double dA;
206 
207  if (A == NULL){
208  A = ARRAY_2D(NFLX, NFLX, double);
209  ALR = ARRAY_2D(NFLX, NFLX, double);
210  #if COMPONENTS != 3
211  print1 ("! PrimEigenvectors(): eigenvector check requires 3 components\n");
212  #endif
213  }
214  #if COMPONENTS != 3
215  return;
216  #endif
217 
218  /* --------------------------------------
219  Construct the Jacobian analytically
220  -------------------------------------- */
221 
222  for (i = 0; i < NFLX; i++){
223  for (j = 0; j < NFLX; j++){
224  A[i][j] = ALR[i][j] = 0.0;
225  }}
226 
227  #if HAVE_ENERGY
228  A[RHO][RHO] = q[VXn] ; A[RHO][VXn] = q[RHO];
229  A[VXn][VXn] = q[VXn] ; A[VXn][PRS] = 1.0/q[RHO];
230  A[VXt][VXt] = q[VXn] ;
231  A[VXb][VXb] = q[VXn] ;
232  A[PRS][VXn] = cs2*q[RHO]; A[PRS][PRS] = q[VXn];
233  #elif EOS == ISOTHERMAL
234  A[RHO][RHO] = q[VXn]; A[RHO][VXn] = q[RHO];
235  A[VXn][VXn] = q[VXn];
236  A[VXn][RHO] = cs2/q[RHO];
237  A[VXt][VXt] = q[VXn];
238  A[VXb][VXb] = q[VXn];
239  #endif
240 
241  for (i = 0; i < NFLX; i++){
242  for (j = 0; j < NFLX; j++){
243  ALR[i][j] = 0.0;
244  for (k = 0; k < NFLX; k++) ALR[i][j] += RR[i][k]*lambda[k]*LL[k][j];
245  }}
246 
247  for (i = 0; i < NFLX; i++){
248  for (j = 0; j < NFLX; j++){
249  dA = ALR[i][j] - A[i][j];
250  if (fabs(dA) > 1.e-8){
251  print ("! PrimEigenvectors: eigenvectors not consistent\n");
252  print ("! A[%d][%d] = %16.9e, R.Lambda.L[%d][%d] = %16.9e\n",
253  i,j, A[i][j], i,j,ALR[i][j]);
254  print ("! cs2 = %12.6e\n",cs2);
255  print ("\n\n A = \n"); ShowMatrix(A, NFLX, 1.e-8);
256  print ("\n\n R.Lambda.L = \n"); ShowMatrix(ALR, NFLX, 1.e-8);
257  QUIT_PLUTO(1);
258  }
259  }}
260 
261 /* -- check orthornomality -- */
262 
263  for (i = 0; i < NFLX; i++){
264  for (j = 0; j < NFLX; j++){
265  a = 0.0;
266  for (k = 0; k < NFLX; k++) a += LL[i][k]*RR[k][j];
267  if ( (i == j && fabs(a-1.0) > 1.e-8) ||
268  (i != j && fabs(a)>1.e-8) ) {
269  print ("! PrimEigenvectors: Eigenvectors not orthogonal\n");
270  print ("! i,j = %d, %d %12.6e \n",i,j,a);
271  print ("! g_dir: %d\n",g_dir);
272  QUIT_PLUTO(1);
273  }
274  }}
275 }
276 #endif
277 }
278 /* ********************************************************************* */
279 void ConsEigenvectors (double *u, double *v, double a2,
280  double **LL, double **RR, double *lambda)
281 /*!
282  * Provide conservative eigenvectors for HD equations.
283  *
284  * \param [in] u array of conservative variables
285  * \param [in] v array of primitive variables
286  * \param [in] a2 square of sound speed
287  * \param [out] LL left conservative eigenvectors
288  * \param [out] RR right conservative eigenvectors
289  * \param [out] lambda eigenvalues
290  *
291  * \b References:
292  * - "Riemann Solvers and Numerical Methods for Fluid Dynamics", \n
293  * Toro, 1997 Springer-Verlag (Page 107)
294  *
295  *********************************************************************** */
296 {
297  int i, j, k, nv;
298  double H, a, gt_1, vmag2;
299 
300  #if EOS == PVTE_LAW
301  print1( "! ConsEigenvectors: cannot be used presently\n");
302  QUIT_PLUTO(1);
303  #endif
304 
305 /* --------------------------------------------------------------
306  If eigenvector check is required, we make sure that
307  U = U(V) is exactly converted from V.
308  This is not the case, with the present version of the code,
309  with FD scheme where V = 0.5*(VL + VR) and U = 0.5*(UL + UR)
310  hence U != U(V).
311  --------------------------------------------------------------- */
312 
313  #if CHECK_EIGENVECTORS == YES
314  u[RHO] = v[RHO];
315  u[MX1] = v[RHO]*v[VX1];
316  u[MX2] = v[RHO]*v[VX2];
317  u[MX3] = v[RHO]*v[VX3];
318  #if EOS == IDEAL
319  u[ENG] = v[PRS]/(g_gamma-1.0)
320  + 0.5*v[RHO]*(v[VX1]*v[VX1] + v[VX2]*v[VX2] + v[VX3]*v[VX3]);
321  #endif
322  #endif
323 
324  #if EOS == IDEAL
325  gt_1 = 1.0/(g_gamma - 1.0);
326 /*
327  a2 = g_gamma*v[PRS]/v[RHO];
328  a = sqrt(a2);
329 */
330  #elif EOS == ISOTHERMAL
331 /*
332  a2 = g_isoSoundSpeed2;
333  a = g_isoSoundSpeed;
334 */
335  #endif
336  a = sqrt(a2);
337 
338  vmag2 = EXPAND(v[VXn]*v[VXn], + v[VXt]*v[VXt], + v[VXb]*v[VXb]);
339 
340 /* -----------------------------------------------
341  to compute H we use primitive variable since
342  U and V may have been averaged in WENO_RF and
343  are not the map of each other.
344  Otherwise left and right eigenvectors wouldn't
345  be orthonormal.
346  ----------------------------------------------- */
347 /*
348  H = (u[ENG] + v[PRS])/v[RHO];
349 */
350  #if EOS == IDEAL
351  H = 0.5*vmag2 + a2*gt_1;
352  #endif
353 
354 /* ======================================================
355  RIGHT EIGENVECTORS
356  ====================================================== */
357 
358 /* lambda = u - c */
359 
360  k = 0;
361  lambda[k] = v[VXn] - a;
362 
363  RR[RHO][k] = 1.0;
364  EXPAND(RR[MXn][k] = lambda[k]; ,
365  RR[MXt][k] = v[VXt]; ,
366  RR[MXb][k] = v[VXb]; )
367  #if EOS == IDEAL
368  RR[ENG][k] = H - a*v[VXn];
369  #endif
370 
371 /* lambda = u + c */
372 
373  k = 1;
374  lambda[k] = v[VXn] + a;
375 
376  RR[RHO][k] = 1.0;
377  EXPAND(RR[MXn][k] = lambda[k]; ,
378  RR[MXt][k] = v[VXt]; ,
379  RR[MXb][k] = v[VXb];)
380  #if EOS == IDEAL
381  RR[ENG][k] = H + a*v[VXn];
382  #endif
383 
384 /* lambda = u */
385 
386  #if EOS == IDEAL
387  k = 2;
388  lambda[k] = v[VXn];
389 
390  RR[RHO][k] = 1.0;
391  EXPAND(RR[MXn][k] = v[VXn]; ,
392  RR[MXt][k] = v[VXt]; ,
393  RR[MXb][k] = v[VXb];)
394  RR[ENG][k] = 0.5*vmag2;
395 
396  #if COMPONENTS > 1
397  k = 3;
398  lambda[k] = v[VXn];
399  RR[MXt][k] = 1.0;
400  RR[ENG][k] = v[VXt];
401  #endif
402 
403  #if COMPONENTS == 3
404  k = 4;
405  lambda[k] = v[VXn];
406  RR[MXb][k] = 1.0;
407  RR[ENG][k] = v[VXb];
408  #endif
409  #elif EOS == ISOTHERMAL
410  EXPAND( ,
411  lambda[2] = v[VXn]; RR[MXt][2] = 1.0; ,
412  lambda[3] = v[VXn]; RR[MXb][3] = 1.0;)
413  #endif
414 
415 /* ======================================================
416  LEFT EIGENVECTORS
417  ====================================================== */
418 
419 /* lambda = u - c */
420 
421  k = 0;
422  #if EOS == IDEAL
423  LL[k][RHO] = H + a*gt_1*(v[VXn] - a);
424  EXPAND(LL[k][MXn] = -(v[VXn] + a*gt_1); ,
425  LL[k][MXt] = -v[VXt]; ,
426  LL[k][MXb] = -v[VXb];)
427  LL[k][ENG] = 1.0;
428  #elif EOS == ISOTHERMAL
429  LL[k][RHO] = 0.5*(1.0 + v[VXn]/a);
430  LL[k][MXn] = -0.5/a;
431  #endif
432 
433 /* lambda = u + c */
434 
435  k = 1;
436  #if EOS == IDEAL
437  LL[k][RHO] = H - a*gt_1*(v[VXn] + a);
438  EXPAND(LL[k][MXn] = -v[VXn] + a*gt_1; ,
439  LL[k][MXt] = -v[VXt]; ,
440  LL[k][MXb] = -v[VXb];)
441  LL[k][ENG] = 1.0;
442  #elif EOS == ISOTHERMAL
443  LL[k][RHO] = 0.5*(1.0 - v[VXn]/a);
444  LL[k][MXn] = 0.5/a;
445  #endif
446 
447 /* lambda = u */
448 
449  #if EOS == IDEAL
450  k = 2;
451  LL[k][RHO] = -2.0*H + 4.0*gt_1*a2;
452  EXPAND(LL[k][MXn] = 2.0*v[VXn]; ,
453  LL[k][MXt] = 2.0*v[VXt]; ,
454  LL[k][MXb] = 2.0*v[VXb];)
455  LL[k][ENG] = -2.0;
456 
457  #if COMPONENTS > 1
458  k = 3;
459  LL[k][RHO] = -2.0*v[VXt]*a2*gt_1;
460  LL[k][MXt] = 2.0*a2*gt_1;
461  LL[k][ENG] = 0.0;
462  #endif
463 
464  #if COMPONENTS == 3
465  k = 4;
466  LL[k][RHO] = -2.0*v[VXb]*a2*gt_1;
467  LL[k][MXb] = 2.0*a2*gt_1;
468  #endif
469 
470  for (k = 0; k < NFLX; k++){ /* normalization */
471  for (nv = 0; nv < NFLX; nv++){
472  LL[k][nv] *= (g_gamma - 1.0)/(2.0*a2);
473  }}
474 
475  #elif EOS == ISOTHERMAL
476  EXPAND( ,
477  k = 2; LL[k][RHO] = -v[VXt]; LL[k][VXt] = 1.0; ,
478  k = 3; LL[k][RHO] = -v[VXb]; LL[k][VXb] = 1.0; )
479  #endif
480 
481 /* -----------------------------------------
482  Check eigenvectors consistency
483  ----------------------------------------- */
484 
485 #if CHECK_EIGENVECTORS == YES
486 {
487  static double **A, **ALR;
488  double dA, vel2, Bmag2, vB;
489 
490  if (A == NULL){
491  A = ARRAY_2D(NFLX, NFLX, double);
492  ALR = ARRAY_2D(NFLX, NFLX, double);
493  #if COMPONENTS != 3
494  print ("! ConsEigenvectors: eigenvector check requires 3 components\n");
495  #endif
496  }
497  #if COMPONENTS != 3
498  return;
499  #endif
500 
501  /* --------------------------------------
502  Construct the Jacobian analytically
503  -------------------------------------- */
504 
505  for (i = 0; i < NFLX; i++){
506  for (j = 0; j < NFLX; j++){
507  A[i][j] = ALR[i][j] = 0.0;
508  }}
509 
510  vel2 = v[VXn]*v[VXn] + v[VXt]*v[VXt] + v[VXb]*v[VXb];
511  #if EOS == IDEAL
512  A[RHO][MXn] = 1.0;
513 
514  A[MXn][RHO] = -v[VXn]*v[VXn] + (g_gamma - 1.0)*vel2*0.5;
515  A[MXn][MXn] = (3.0 - g_gamma)*v[VXn];
516  A[MXn][MXt] = (1.0 - g_gamma)*v[VXt];
517  A[MXn][MXb] = (1.0 - g_gamma)*v[VXb];
518  A[MXn][ENG] = g_gamma - 1.0;
519 
520  A[MXt][RHO] = - v[VXn]*v[VXt];
521  A[MXt][MXn] = v[VXt];
522  A[MXt][MXt] = v[VXn];
523 
524  A[MXb][RHO] = - v[VXn]*v[VXb];
525  A[MXb][MXn] = v[VXb];
526  A[MXb][MXb] = v[VXn];
527 
528  A[ENG][RHO] = 0.5*(g_gamma - 1.0)*vel2*v[VXn] - (u[ENG] + v[PRS])/v[RHO]*v[VXn];
529 
530  A[ENG][MXn] = (1.0 - g_gamma)*v[VXn]*v[VXn] + (u[ENG] + v[PRS])/v[RHO];
531 
532  A[ENG][MXt] = (1.0 - g_gamma)*v[VXn]*v[VXt];
533  A[ENG][MXb] = (1.0 - g_gamma)*v[VXn]*v[VXb];
534 
535  A[ENG][ENG] = g_gamma*v[VXn];
536  #elif EOS == ISOTHERMAL
537  A[RHO][MXn] = 1.0;
538 
539  A[MXn][RHO] = -v[VXn]*v[VXn] + a2;
540  A[MXn][MXn] = 2.0*v[VXn];
541 
542  A[MXt][RHO] = - v[VXn]*v[VXt];
543  A[MXt][MXn] = v[VXt];
544  A[MXt][MXt] = v[VXn];
545 
546  A[MXb][RHO] = - v[VXn]*v[VXb];
547  A[MXb][MXn] = v[VXb];
548  A[MXb][MXb] = v[VXn];
549  #endif
550 
551  for (i = 0; i < NFLX; i++){
552  for (j = 0; j < NFLX; j++){
553  ALR[i][j] = 0.0;
554  for (k = 0; k < NFLX; k++){
555  ALR[i][j] += RR[i][k]*lambda[k]*LL[k][j];
556  }
557  }}
558 
559  for (i = 0; i < NFLX; i++){
560  for (j = 0; j < NFLX; j++){
561  if (fabs(ALR[i][j] - A[i][j]) > 1.e-6){
562  print ("! ConsEigenvectors: eigenvectors not consistent\n");
563  print ("! g_dir = %d\n",g_dir);
564  print ("! A[%d][%d] = %16.9e, R.Lambda.L[%d][%d] = %16.9e\n",
565  i,j, A[i][j], i,j,ALR[i][j]);
566  print ("\n\n A = \n"); ShowMatrix(A, NFLX, 1.e-8);
567  print ("\n\n R.Lambda.L = \n"); ShowMatrix(ALR, NFLX, 1.e-8);
568 
569  print ("\n\n RR = \n"); ShowMatrix(RR, NFLX, 1.e-8);
570  print ("\n\n LL = \n"); ShowMatrix(LL, NFLX, 1.e-8);
571  QUIT_PLUTO(1);
572  }
573  }}
574 
575 /* -- check orthornomality -- */
576 
577  for (i = 0; i < NFLX; i++){
578  for (j = 0; j < NFLX; j++){
579  dA = 0.0;
580  for (k = 0; k < NFLX; k++) dA += LL[i][k]*RR[k][j];
581  if ( (i == j && fabs(dA-1.0) > 1.e-8) ||
582  (i != j && fabs(dA)>1.e-8) ) {
583  print ("! ConsEigenvectors: Eigenvectors not orthogonal\n");
584  print ("! i,j = %d, %d %12.6e \n",i,j,dA);
585  print ("! g_dir: %d\n",g_dir);
586  QUIT_PLUTO(1);
587  }
588  }}
589 }
590 #endif
591 }
592 /* ********************************************************************* */
593 void PrimToChar (double **L, double *v, double *w)
594 /*!
595  * Compute the matrix-vector multiplcation between the
596  * the L matrix (containing primitive left eigenvectors)
597  * and the vector v. The result containing the characteristic
598  * variables is stored in the vector w = L.v
599  *
600  * For efficiency purpose, multiplication is done
601  * explicitly, so that only nonzero entries
602  * of the left primitive eigenvectors are considered.
603  *
604  * \param [in] L Left eigenvectors
605  * \param [in] v (difference of) primitive variables
606  * \param [out] w (difference of) characteristic variables
607  *
608  *********************************************************************** */
609 {
610  int nv;
611 
612  #if HAVE_ENERGY
613  w[0] = L[0][VXn]*v[VXn] + L[0][PRS]*v[PRS];
614  w[1] = L[1][VXn]*v[VXn] + L[1][PRS]*v[PRS];
615  EXPAND( w[2] = v[RHO] + L[2][PRS]*v[PRS]; ,
616  w[3] = v[VXt]; ,
617  w[4] = v[VXb];)
618  #elif EOS == ISOTHERMAL
619  w[0] = L[0][RHO]*v[RHO] + L[0][VXn]*v[VXn];
620  w[1] = L[1][RHO]*v[RHO] + L[1][VXn]*v[VXn];
621  EXPAND( ,
622  w[2] = v[VXt]; ,
623  w[3] = v[VXb];)
624  #endif
625 
626 /* -----------------------------------------------
627  For passive scalars, the characteristic
628  variable is equal to the primitive one,
629  since l = r = (0,..., 1 , 0 ,....)
630  ----------------------------------------------- */
631 
632 #if NSCL > 0
633  NSCL_LOOP(nv) w[nv] = v[nv];
634 #endif
635 }
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
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
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
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
NSCL
#define NSCL
Definition: pluto.h:572
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
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MXn
int MXn
Definition: globals.h:74
PrimToChar
void PrimToChar(double **L, double *v, double *w)
Definition: eigenv.c:593
VXn
int VXn
Definition: globals.h:73
CHECK_EIGENVECTORS
#define CHECK_EIGENVECTORS
Definition: pluto.h:411
ISOTHERMAL
#define ISOTHERMAL
Definition: pluto.h:49
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
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
q
static Runtime q
Definition: runtime_setup.c:490
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
MXb
int MXb
Definition: globals.h:74