PLUTO
Macros | Enumerations | Functions | Variables
mod_defs.h File Reference

Set labels, indexes and prototypes for the HD module. More...

Go to the source code of this file.

Macros

#define RHO   0
 
#define MX1   1
 
#define MX2   (COMPONENTS >= 2 ? 2: 255)
 
#define MX3   (COMPONENTS == 3 ? 3: 255)
 
#define VX1   MX1
 
#define VX2   MX2
 
#define VX3   MX3
 
#define NFLX   (1 + COMPONENTS + HAVE_ENERGY)
 
#define iVR   VX1
 
#define iVZ   VX2
 
#define iVPHI   VX3
 
#define iMR   MX1
 
#define iMZ   MX2
 
#define iMPHI   MX3
 
#define iVR   VX1
 
#define iVPHI   VX2
 
#define iVZ   VX3
 
#define iMR   MX1
 
#define iMPHI   MX2
 
#define iMZ   MX3
 
#define iVR   VX1
 
#define iVTH   VX2
 
#define iVPHI   VX3
 
#define iMR   MX1
 
#define iMTH   MX2
 
#define iMPHI   MX3
 

Enumerations

enum  KWAVES {
  KSOUNDM, KSOUNDP, KFASTM, KFASTP,
  KPSI_GLMM, KPSI_GLMP, KFASTM, KFASTP,
  KENTRP, KPSI_GLMM, KPSI_GLMP
}
 

Functions

int ConsToPrim (double **, double **, int, int, unsigned char *)
 
void Eigenvalues (double **, double *, double **, int, int)
 
void PrimEigenvectors (double *, double, double, double *, double **, double **)
 
void ConsEigenvectors (double *, double *, double, double **, double **, double *)
 
void Flux (double **, double **, double *, double **, double *, int, int)
 
void HLL_Speed (double **, double **, double *, double *, double *, double *, int, int)
 
void MaxSignalSpeed (double **, double *, double *, double *, int, int)
 
void PrimToCons (double **, double **, int, int)
 
void PrimRHS (double *, double *, double, double, double *)
 
void PrimSource (const State_1D *, int, int, double *, double *, double **, Grid *)
 

Variables

Riemann_Solver TwoShock_Solver
 
Riemann_Solver LF_Solver
 
Riemann_Solver Roe_Solver
 
Riemann_Solver HLL_Solver
 
Riemann_Solver HLLC_Solver
 
Riemann_Solver RusanovDW_Solver
 
Riemann_Solver AUSMp_Solver
 

Detailed Description

Set labels, indexes and prototypes for the HD module.

Contains variable names and prototypes for the HD module

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
April, 2, 2015

Definition in file mod_defs.h.

Macro Definition Documentation

#define iMPHI   MX3

Definition at line 71 of file mod_defs.h.

#define iMPHI   MX2

Definition at line 71 of file mod_defs.h.

#define iMPHI   MX3

Definition at line 71 of file mod_defs.h.

#define iMR   MX1

Definition at line 69 of file mod_defs.h.

#define iMR   MX1

Definition at line 69 of file mod_defs.h.

#define iMR   MX1

Definition at line 69 of file mod_defs.h.

#define iMTH   MX2

Definition at line 70 of file mod_defs.h.

#define iMZ   MX2

Definition at line 59 of file mod_defs.h.

#define iMZ   MX3

Definition at line 59 of file mod_defs.h.

#define iVPHI   VX3

Definition at line 67 of file mod_defs.h.

#define iVPHI   VX2

Definition at line 67 of file mod_defs.h.

#define iVPHI   VX3

Definition at line 67 of file mod_defs.h.

#define iVR   VX1

Definition at line 65 of file mod_defs.h.

#define iVR   VX1

Definition at line 65 of file mod_defs.h.

#define iVR   VX1

Definition at line 65 of file mod_defs.h.

#define iVTH   VX2

Definition at line 66 of file mod_defs.h.

#define iVZ   VX2

Definition at line 55 of file mod_defs.h.

#define iVZ   VX3

Definition at line 55 of file mod_defs.h.

#define MX1   1

Definition at line 20 of file mod_defs.h.

#define MX2   (COMPONENTS >= 2 ? 2: 255)

Definition at line 21 of file mod_defs.h.

#define MX3   (COMPONENTS == 3 ? 3: 255)

Definition at line 22 of file mod_defs.h.

#define NFLX   (1 + COMPONENTS + HAVE_ENERGY)

Definition at line 32 of file mod_defs.h.

#define RHO   0

Definition at line 19 of file mod_defs.h.

#define VX1   MX1

Definition at line 28 of file mod_defs.h.

#define VX2   MX2

Definition at line 29 of file mod_defs.h.

#define VX3   MX3

Definition at line 30 of file mod_defs.h.

Enumeration Type Documentation

enum KWAVES
Enumerator
KSOUNDM 
KSOUNDP 
KFASTM 
KFASTP 
KPSI_GLMM 
KPSI_GLMP 
KFASTM 
KFASTP 
KENTRP 
KPSI_GLMM 
KPSI_GLMP 

Definition at line 80 of file mod_defs.h.

80  {
81  KSOUNDM, KSOUNDP
82  #if HAVE_ENERGY
83  , KENTRP
84  #endif
85 };
KSOUNDP
Definition: mod_defs.h:81
KSOUNDM
Definition: mod_defs.h:81
KENTRP
Definition: mod_defs.h:59

Function Documentation

void ConsEigenvectors ( double *  u,
double *  v,
double  a2,
double **  Lc,
double **  Rc,
double *  lambda 
)

Provide conservative eigenvectors for HD equations.

Parameters
[in]uarray of conservative variables
[in]varray of primitive variables
[in]a2square of sound speed
[out]LLleft conservative eigenvectors
[out]RRright conservative eigenvectors
[out]lambdaeigenvalues

References:

  • "Riemann Solvers and Numerical Methods for Fluid Dynamics",
    Toro, 1997 Springer-Verlag (Page 107)

Provide conservative eigenvectors for MHD equations.

Parameters
[in]uarray of conservative variables
[in]varray of primitive variables
[in]a2square of sound speed
[out]Lcleft conservative eigenvectors
[out]Rcright conservative eigenvectors
[out]lambdaeigenvalues

Reference

  • "High-order conservative finite difference GLM-MHD schemes for cell-centered MHD"
    Mignone, Tzeferacos & Bodo, JCP (2010) 229, 5896
  • "A High-order WENO Finite Difference Scheme for the Equations of Ideal MHD"
    Jiang,Wu, JCP 150,561 (1999)

With the following corrections:

The components (By, Bz) of L_{1,7} page 571 should be +{rho}a instead of -. Also, since Bx is not considered as a parameter, one must also add a component in the fast and slow waves. This can be seen by forming the conservative eigenvectors from the primitive ones, see the paper from Powell, Roe Linde.

The Lc_{1,3,5,7}[Bx] components, according to Serna 09 are -(1-g_gamma)*a_{f,s}*Bx and not (1-g_gamma)*a_{f,s}*Bx. Both are orthonormal though. She is WRONG! -Petros-

Definition at line 279 of file eigenv.c.

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 }
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
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
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
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
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
VXn
int VXn
Definition: globals.h:73
CHECK_EIGENVECTORS
#define CHECK_EIGENVECTORS
Definition: pluto.h:411
i
int i
Definition: analysis.c:2
VX3
#define VX3
Definition: mod_defs.h:30
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
MXb
int MXb
Definition: globals.h:74

Here is the call graph for this function:

Here is the caller graph for this function:

int ConsToPrim ( double **  ucons,
double **  uprim,
int  beg,
int  end,
unsigned char *  flag 
)

Convert from conservative to primitive variables.

Parameters
[in]uconsarray of conservative variables
[out]uprimarray of primitive variables
[in]begstarting index of computation
[in]endfinal index of computation
[in,out]flagarray of flags tagging, in input, zones where entropy must be used to recover pressure and, on output, zones where conversion was not successful.
Returns
Return 0 if conversion was successful in all zones in the range [ibeg,iend]. Return 1 if one or more zones could not be converted correctly and either pressure, density or energy took non-physical values.

Definition at line 89 of file mappers.c.

109 {
110  int i, nv, err, ifail;
111  int use_entropy, use_energy=1;
112  double tau, rho, gmm1, rhoe, T;
113  double kin, m2, rhog1;
114  double *u, *v;
115 
116 #if EOS == IDEAL
117  gmm1 = g_gamma - 1.0;
118 #endif
119 
120  ifail = 0;
121  for (i = ibeg; i <= iend; i++) {
122 
123  u = ucons[i];
124  v = uprim[i];
125 
126  m2 = EXPAND(u[MX1]*u[MX1], + u[MX2]*u[MX2], + u[MX3]*u[MX3]);
127 
128  /* -- Check density positivity -- */
129 
130  if (u[RHO] < 0.0) {
131  print("! ConsToPrim: rho < 0 (%8.2e), ", u[RHO]);
132  Where (i, NULL);
133  u[RHO] = g_smallDensity;
134  flag[i] |= FLAG_CONS2PRIM_FAIL;
135  ifail = 1;
136  }
137 
138  /* -- Compute density, velocity and scalars -- */
139 
140  v[RHO] = rho = u[RHO];
141  tau = 1.0/u[RHO];
142  EXPAND(v[VX1] = u[MX1]*tau; ,
143  v[VX2] = u[MX2]*tau; ,
144  v[VX3] = u[MX3]*tau;)
145 
146  kin = 0.5*m2/u[RHO];
147 
148  /* -- Check energy positivity -- */
149 
150 #if HAVE_ENERGY
151  if (u[ENG] < 0.0) {
152  WARNING(
153  print("! ConsToPrim: E < 0 (%8.2e), ", u[ENG]);
154  Where (i, NULL);
155  )
156  u[ENG] = g_smallPressure/gmm1 + kin;
157  flag[i] |= FLAG_CONS2PRIM_FAIL;
158  ifail = 1;
159  }
160 #endif
161 
162  /* -- Compute pressure from total energy or entropy -- */
163 
164 #if EOS == IDEAL
165  #if ENTROPY_SWITCH
166  use_entropy = (flag[i] & FLAG_ENTROPY);
167  use_energy = !use_entropy;
168  if (use_entropy){
169  rhog1 = pow(rho, gmm1);
170  v[PRS] = u[ENTR]*rhog1;
171  if (v[PRS] < 0.0){
172  WARNING(
173  print("! ConsToPrim: p(S) < 0 (%8.2e, %8.2e), ", v[PRS], u[ENTR]);
174  Where (i, NULL);
175  )
176  v[PRS] = g_smallPressure;
177  flag[i] |= FLAG_CONS2PRIM_FAIL;
178  ifail = 1;
179  }
180  u[ENG] = v[PRS]/gmm1 + kin; /* -- redefine energy -- */
181  }
182  #endif /* ENTROPY_SWITCH */
183 
184  if (use_energy){
185  v[PRS] = gmm1*(u[ENG] - kin);
186  if (v[PRS] < 0.0){
187  WARNING(
188  print("! ConsToPrim: p(E) < 0 (%8.2e), ", v[PRS]);
189  Where (i, NULL);
190  )
191  v[PRS] = g_smallPressure;
192  u[ENG] = v[PRS]/gmm1 + kin; /* -- redefine energy -- */
193  flag[i] |= FLAG_CONS2PRIM_FAIL;
194  ifail = 1;
195  }
196  #if ENTROPY_SWITCH
197  u[ENTR] = v[PRS]/pow(rho,gmm1);
198  #endif
199  }
200 
201  #if NSCL > 0
202  NSCL_LOOP(nv) v[nv] = u[nv]*tau;
203  #endif
204 
205 #elif EOS == PVTE_LAW
206 
207  /* -- Convert scalars here since EoS may need ion fractions -- */
208 
209  #if NSCL > 0
210  NSCL_LOOP(nv) v[nv] = u[nv]*tau;
211  #endif
212 
213  if (u[ENG] != u[ENG]){
214  print("! ConsToPrim: NaN found\n");
215  Show(ucons,i);
216  QUIT_PLUTO(1);
217  }
218  rhoe = u[ENG] - kin;
219 
220  err = GetEV_Temperature (rhoe, v, &T);
221  if (err){ /* If something went wrong while retrieving the */
222  /* temperature, we floor \c T to \c T_CUT_RHOE, */
223  /* recompute internal and total energies. */
224  T = T_CUT_RHOE;
225  WARNING(
226  print ("! ConsToPrim: rhoe < 0 or T < T_CUT_RHOE; "); Where(i,NULL);
227  )
228  rhoe = InternalEnergy(v, T);
229  u[ENG] = rhoe + kin; /* -- redefine total energy -- */
230  flag[i] |= FLAG_CONS2PRIM_FAIL;
231  ifail = 1;
232  }
233  v[PRS] = Pressure(v, T);
234 
235 #endif /* EOS == PVTE_LAW */
236 
237  /* -- Dust-- */
238 
239 #if DUST == YES
240  u[RHO_D] = MAX(u[RHO_D], 1.e-20);
241  v[RHO_D] = u[RHO_D];
242  EXPAND(v[MX1_D] = u[MX1_D]/v[RHO_D]; ,
243  v[MX2_D] = u[MX2_D]/v[RHO_D]; ,
244  v[MX3_D] = u[MX3_D]/v[RHO_D];)
245 #endif
246 
247  }
248  return ifail;
249 }
MX3
#define MX3
Definition: mod_defs.h:22
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
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
T_CUT_RHOE
#define T_CUT_RHOE
Sets the lowest cut-off temperature used in the PVTE_LAW equation of state.
Definition: eos.h:69
g_smallDensity
double g_smallDensity
Small value for density fix.
Definition: globals.h:109
InternalEnergy
double InternalEnergy(double *v, double T)
Definition: internal_energy.c:255
RHO
#define RHO
Definition: mod_defs.h:19
WARNING
#define WARNING(a)
Definition: macros.h:217
func_param::rhoe
double rhoe
Definition: eos.h:107
func_param::v
double v[NVAR]
Definition: eos.h:106
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
VX1
#define VX1
Definition: mod_defs.h:28
FLAG_CONS2PRIM_FAIL
#define FLAG_CONS2PRIM_FAIL
Definition: pluto.h:189
FLAG_ENTROPY
#define FLAG_ENTROPY
Update pressure using entropy equation.
Definition: pluto.h:182
GetEV_Temperature
int GetEV_Temperature(double rhoe, double *v, double *T)
Definition: internal_energy.c:290
MX2
#define MX2
Definition: mod_defs.h:21
Where
void Where(int, Grid *)
Definition: tools.c:235
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
Show
void Show(double **, int)
Definition: tools.c:132
i
int i
Definition: analysis.c:2
Pressure
double Pressure(double *v, double T)
Definition: thermal_eos.c:179
VX3
#define VX3
Definition: mod_defs.h:30
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
void Eigenvalues ( double **  v,
double *  csound2,
double **  lambda,
int  beg,
int  end 
)

Compute eigenvalues for the HD equations

Parameters
[in]v1-D array of primitive variables
[out]csound21-D array containing the square of sound speed
[out]lambda1-D array [i][nv] containing the eigenvalues
[in]begstarting index of computation
[in]endfinal index of computation

Compute eigenvalues for the MHD equations

Parameters
[in]v1-D array of primitive variables
[out]csound21-D array containing the square of sound speed
[out]lambda1-D array [i][nv] containing the eigenvalues
[in]begstarting index of computation
[in]endfinal index of computation

Definition at line 67 of file eigenv.c.

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 }
func_param::v
double v[NVAR]
Definition: eos.h:106
NFLX
#define NFLX
Definition: mod_defs.h:32
k
int k
Definition: analysis.c:2
VXn
int VXn
Definition: globals.h:73
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

void Flux ( double **  u,
double **  w,
double *  a2,
double **  fx,
double *  p,
int  beg,
int  end 
)
Parameters
[in]u1D array of conserved quantities
[in]w1D array of primitive quantities
[in]a21D array of sound speeds
[out]fx1D array of fluxes (total pressure excluded)
[out]p1D array of pressure values
[in]beginitial index of computation
[in]endfinal index of computation
Returns
This function has no return value.

Definition at line 23 of file fluxes.c.

36 {
37  int nv, ii;
38 
39  for (ii = beg; ii <= end; ii++) {
40  fx[ii][RHO] = u[ii][MXn];
41  EXPAND(fx[ii][MX1] = u[ii][MX1]*w[ii][VXn]; ,
42  fx[ii][MX2] = u[ii][MX2]*w[ii][VXn]; ,
43  fx[ii][MX3] = u[ii][MX3]*w[ii][VXn];)
44 #if HAVE_ENERGY
45  p[ii] = w[ii][PRS];
46  fx[ii][ENG] = (u[ii][ENG] + w[ii][PRS])*w[ii][VXn];
47 #elif EOS == ISOTHERMAL
48  p[ii] = a2[ii]*w[ii][RHO];
49 #endif
50 /*
51 #if DUST == YES
52  fx[ii][RHO_D] = u[ii][MXn_D];
53  EXPAND(fx[ii][MX1_D] = u[ii][MX1_D]*w[ii][VXn_D]; ,
54  fx[ii][MX2_D] = u[ii][MX2_D]*w[ii][VXn_D]; ,
55  fx[ii][MX3_D] = u[ii][MX3_D]*w[ii][VXn_D];)
56 #endif
57 */
58  }
59 }
MX3
#define MX3
Definition: mod_defs.h:22
MX1
#define MX1
Definition: mod_defs.h:20
RHO
#define RHO
Definition: mod_defs.h:19
MX2
#define MX2
Definition: mod_defs.h:21
MXn
int MXn
Definition: globals.h:74
VXn
int VXn
Definition: globals.h:73
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395
void HLL_Speed ( double **  vL,
double **  vR,
double *  a2L,
double *  a2R,
double *  SL,
double *  SR,
int  beg,
int  end 
)

Compute leftmost (SL) and rightmost (SR) speed for the Riemann fan.

Parameters
[in]vLleft state for the Riemann solver
[in]vRright state for the Riemann solver
[in]a2L1-D array containing the square of the sound speed for the left state
[in]a2R1-D array containing the square of the sound speed for the right state
[out]SLthe (estimated) leftmost speed of the Riemann fan
[out]SRthe (estimated) rightmost speed of the Riemann fan
[in]begstarting index of computation
[in]endfinal index of computation

Switches:

ROE_ESTIMATE (YES/NO), DAVIS_ESTIMATE (YES/NO). TVD_ESTIMATE (YES/NO) JAN_HLL (YES/NO)

These switches set how the wave speed estimates are going to be computed. Only one can be set to 'YES', and the rest of them must be set to 'NO'

ROE_ESTIMATE: b_m = (0, (u_R - c_R, u_L - c_L, u_{roe} - c_{roe})) b_m = (0, (u_R + c_R, u_L + c_L, u_{roe} + c_{roe}))

where u_{roe} and c_{roe} are computed using Roe averages.

DAVIS_ESTIMATE: b_m = (0, (u_R - c_R, u_L - c_L)) b_m = (0, (u_R + c_R, u_L + c_L))

Compute leftmost (SL) and rightmost (SR) speed for the Riemann fan.

Parameters
[in]vLleft state for the Riemann solver
[in]vRright state for the Riemann solver
[in]a2L1-D array containing the square of the sound speed for the left state
[in]a2R1-D array containing the square of the sound speed for the right state
[out]SLthe (estimated) leftmost speed of the Riemann fan
[out]SRthe (estimated) rightmost speed of the Riemann fan
[in]begstarting index of computation
[in]endfinal index of computation

Definition at line 24 of file hll_speed.c.

58 {
59  int i;
60  double scrh, s, c;
61  double aL, dL;
62  double aR, dR;
63  double dvx, dvy, dvz;
64  double a_av, du, vx;
65  static double *sl_min, *sl_max;
66  static double *sr_min, *sr_max;
67 
68  if (sl_min == NULL){
69  sl_min = ARRAY_1D(NMAX_POINT, double);
70  sl_max = ARRAY_1D(NMAX_POINT, double);
71 
72  sr_min = ARRAY_1D(NMAX_POINT, double);
73  sr_max = ARRAY_1D(NMAX_POINT, double);
74  }
75 
76 /* ----------------------------------------------
77  DAVIS Estimate
78  ---------------------------------------------- */
79 
80  #if DAVIS_ESTIMATE == YES
81 
82  MaxSignalSpeed (vL, a2L, sl_min, sl_max, beg, end);
83  MaxSignalSpeed (vR, a2R, sr_min, sr_max, beg, end);
84  for (i = beg; i <= end; i++) {
85 
86  SL[i] = MIN(sl_min[i], sr_min[i]);
87  SR[i] = MAX(sl_max[i], sr_max[i]);
88 
89  aL = sqrt(a2L[i]);
90  aR = sqrt(a2R[i]);
91 
92  scrh = fabs(vL[i][VXn]) + fabs(vR[i][VXn]);
93  scrh /= aL + aR;
94 
95  g_maxMach = MAX(scrh, g_maxMach);
96  }
97 
98  #endif
99 
100 /* ----------------------------------------------
101  Einfeldt Estimate for wave speed
102  ---------------------------------------------- */
103 
104  #if EINFELDT_ESTIMATE == YES
105 
106  for (i = beg; i <= end; i++) {
107 
108  aL = sqrt(a2L[i]);
109  aR = sqrt(a2R[i]);
110 
111  dL = sqrt(vL[i][RHO]);
112  dR = sqrt(vR[i][RHO]);
113  a_av = 0.5*dL*dR/( (dL + dR)*(dL + dR));
114  du = vR[i][VXn] - vL[i][VXn];
115  scrh = (dR*aL*aL + dR*aR*aR)/(dL + dR);
116  scrh += a_av*du*du;
117 
118  SL[i] = (dl*ql[VXn] + dr*qr[VXn])/(dl + dr);
119 
120  bmin = MIN(0.0, um[VXn] - sqrt(scrh));
121  bmax = MAX(0.0, um[VXn] + sqrt(scrh));
122  }
123  #endif
124 
125 /* ----------------------------------------------
126  Roe-like Estimate
127  ---------------------------------------------- */
128 
129  #if ROE_ESTIMATE == YES
130 
131  for (i = beg; i <= end; i++) {
132 
133  aL = sqrt(a2L[i]);
134  aR = sqrt(a2R[i]);
135 
136  s = 1.0/(1.0 + sqrt(vR[i][RHO]/vL[i][RHO]));
137  c = 1.0 - s;
138 
139  vx = s*vL[i][VXn] + c*vR[i][VXn];
140 
141 /* ************************************************************
142  The following definition of the averaged sound speed
143  is equivalent to original formulation given by Roe;
144  here we just simplify the expression without explicitly
145  computing the Enthalpy:
146  ************************************************************ */
147 
148  EXPAND(dvx = vL[i][VX1] - vR[i][VX1]; ,
149  dvy = vL[i][VX2] - vR[i][VX2]; ,
150  dvz = vL[i][VX3] - vR[i][VX3];)
151 
152  scrh = EXPAND(dvx*dvx, + dvy*dvy, + dvz*dvz);
153 
154  a_av = sqrt(s*aL2 + c*aR2 + 0.5*s*c*gmm1*scrh);
155 
156  SL[i] = MIN(vL[i][VXn] - aL, vx - a_av);
157  SR[i] = MAX(vR[i][VXn] + aR, vx + a_av);
158 
159  scrh = (fabs(vL[i][VXn]) + fabs(vR[i][VXn]))/(aL + aR);
160  g_maxMach = MAX(scrh, g_maxMach);
161 
162  }
163  #endif
164 
165 }
MAX
#define MAX(a, b)
Definition: macros.h:101
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
VX1
#define VX1
Definition: mod_defs.h:28
MIN
#define MIN(a, b)
Definition: macros.h:104
g_maxMach
double g_maxMach
The maximum Mach number computed during integration.
Definition: globals.h:119
menu.c
tuple c
Definition: menu.py:375
VXn
int VXn
Definition: globals.h:73
s
#define s
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
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
void MaxSignalSpeed ( double **  v,
double *  cs2,
double *  cmin,
double *  cmax,
int  beg,
int  end 
)

Compute the maximum and minimum characteristic velocities for the HD equation from the vector of primitive variables v.

Parameters
[in]v1-D array of primitive variables
[in]cs21-D array containing the square of the sound speed
[out]cmin1-D array containing the leftmost characteristic
[out]cmin1-D array containing the rightmost characteristic
[in]begstarting index of computation
[in]endfinal index of computation

Definition at line 34 of file eigenv.c.

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 }
a
static double a
Definition: init.c:135
func_param::v
double v[NVAR]
Definition: eos.h:106
VXn
int VXn
Definition: globals.h:73
i
int i
Definition: analysis.c:2
void PrimEigenvectors ( double *  q,
double  a2,
double  h,
double *  lambda,
double **  LL,
double **  RR 
)

Provide left and right eigenvectors and corresponding eigenvalues for the primitive form of the HD equations (adiabatic, pvte & isothermal cases).

Parameters
[in]qvector of primitive variables
[in]cs2sound speed squared
[in]henthalpy
[out]lambdaeigenvalues
[out]LLmatrix containing left primitive eigenvectors (rows)
[out]RRmatrix containing right primitive eigenvectors (columns)
Note
It is highly recommended that LL and RR be initialized to zero BEFORE since only non-zero entries are treated here.

Wave names and their order are defined as enumeration constants in mod_defs.h.

Advection modes associated with passive scalars are simple cases for which lambda = u (entropy mode) and l = r = (0, ... , 1, 0, ...). For this reason they are NOT defined here and must be treated seperately.

References:

  • "Riemann Solvers and Numerical Methods for Fluid Dynamics",
    Toro, 1997 Springer-Verlag, Eq. [3.18]

Provide left and right eigenvectors and corresponding eigenvalues for the PRIMITIVE form of the MHD equations (adiabatic & isothermal cases).

Parameters
[in]qvector of primitive variables
[in]a2sound speed squared
[in]henthalpy
[out]lambdaeigenvalues
[out]LLleft primitive eigenvectors
[out]RRright primitive eigenvectors
Note
It is highly recommended that LL and RR be initialized to zero before since only non-zero entries are treated here.

Wave names and their order are defined as enumeration constants in mod_defs.h. Notice that, the characteristic decomposition may differ depending on the way div.B is treated.

Advection modes associated with passive scalars are simple cases for which lambda = u (entropy mode) and l = r = (0, ... , 1, 0, ...). For this reason they are NOT defined here and must be treated seperately.

References:

  • "Notes on the Eigensystem of Magnetohydrodynamics",
    P.L. Roe, D.S. Balsara, SIAM Journal on Applied Mathematics, 56, 57 (1996)
  • "A solution adaptive upwind scheme for ideal MHD",
    K. G. Powell, P. L. Roe, and T. J. Linde, Journal of Computational Physics, 154, 284-309, (1999).
  • "A second-order unsplit Godunov scheme for cell-centered MHD: the CTU-GLM scheme"
    Mignone & Tzeferacos, JCP (2010) 229, 2117
  • "ATHENA: A new code for astrophysical MHD", J. Stone, T. Gardiner, ApJS, 178, 137 (2008)

The derivation of the isothermal eigenvectors follows the consideration given in roe.c

Definition at line 91 of file eigenv.c.

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 }
EOS
#define EOS
Definition: pluto.h:341
a
static double a
Definition: init.c:135
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
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
VXb
int VXb
Definition: globals.h:73
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
j
int j
Definition: analysis.c:2
VXt
int VXt
Definition: globals.h:73
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
VXn
int VXn
Definition: globals.h:73
CHECK_EIGENVECTORS
#define CHECK_EIGENVECTORS
Definition: pluto.h:411
ISOTHERMAL
#define ISOTHERMAL
Definition: pluto.h:49
i
int i
Definition: analysis.c:2
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
void PrimRHS ( double *  w,
double *  dw,
double  cs2,
double  h,
double *  Adw 
)

Compute the matrix-vector multiplication $ A(\mathbf{v})\cdot \Delta\mathbf{v} $ where A is the matrix of the quasi-linear form of the HD equations.

Parameters
[in]wvector of primitive variables;
[in]dwlimited (linear) slopes;
[in]cs2local sound speed;
[in]hlocal enthalpy;
[out]Adwmatrix-vector product.
Returns
This function has no return value.

Compute the matrix-vector multiplication $ A(\mathbf{v})\cdot d\mathbf{v} $ where A is the matrix of the quasi-linear form of the MHD equations.

References

  • "A solution adaptive upwind scheme for ideal MHD", Powell et al., JCP (1999) 154, 284
  • "An unsplit Godunov method for ideal MHD via constrained transport" Gardiner & Stone, JCP (2005) 205, 509
Parameters
[in]vvector of primitive variables
[in]dvlimited (linear) slopes
[in]cs2local sound speed
[in]hlocal enthalpy
[out]AdVmatrix-vector product
Returns
This function has no return value.

Compute the matrix-vector multiplication $ A(\mathbf{v})\cdot \Delta\mathbf{v} $ where A is the matrix of the quasi-linear form of the RHD equations.

Parameters
[in]wvector of primitive variables;
[in]dwlimited (linear) slopes;
[in]cs2local sound speed;
[in]hlocal enthalpy;
[out]Adwmatrix-vector product.
Note
Returns
This function has no return value.
Note
In the 7-wave and 8-wave formulations we use the the same matrix being decomposed into right and left eigenvectors during the Characteristic Tracing step. Note, however, that it DOES NOT include two additional terms (-vy*dV[BX] for By, -vz*dv[BX] for Bz) that are needed in the

7-wave form and are added using source terms.

Definition at line 31 of file prim_eqn.c.

44 {
45  int nv;
46  double u;
47 
48  u = v[VXn];
49 
50  Adv[RHO] = u*dv[RHO] + v[RHO]*dv[VXn];
51  #if EOS == IDEAL || EOS == PVTE_LAW
52  EXPAND(Adv[VXn] = u*dv[VXn] + dv[PRS]/v[RHO]; ,
53  Adv[VXt] = u*dv[VXt]; ,
54  Adv[VXb] = u*dv[VXb];)
55  Adv[PRS] = cs2*v[RHO]*dv[VXn] + u*dv[PRS];
56  #elif EOS == ISOTHERMAL
57  EXPAND(Adv[VXn] = u*dv[VXn] + cs2*dv[RHO]/v[RHO]; ,
58  Adv[VXt] = u*dv[VXt]; ,
59  Adv[VXb] = u*dv[VXb];)
60  #endif
61 
62  /* ---- scalars ---- */
63 
64 #if NSCL > 0
65  NSCL_LOOP(nv) Adv[nv] = u*dv[nv];
66 #endif
67 
68 }
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
RHO
#define RHO
Definition: mod_defs.h:19
func_param::v
double v[NVAR]
Definition: eos.h:106
VXb
int VXb
Definition: globals.h:73
NSCL
#define NSCL
Definition: pluto.h:572
VXt
int VXt
Definition: globals.h:73
VXn
int VXn
Definition: globals.h:73
void PrimSource ( const State_1D *  state,
int  beg,
int  end,
double *  a2,
double *  h,
double **  src,
Grid *  grid 
)

Compute source terms of the HD equations in primitive variables.

  • Geometrical sources;
  • Gravity;
  • Fargo source terms.

The rationale for choosing during which sweep a particular source term has to be incorporated should match the same criterion used during the conservative update. For instance, in polar or cylindrical coordinates, curvilinear source terms are included during the radial sweep only.

Parameters
[in]statepointer to a State_1D structure;
[in]beginitial index of computation;
[in]endfinal index of computation;
[in]a2array of sound speed;
[in]harray of enthalpies (not needed in MHD);
[out]srcarray of source terms;
[in]gridpointer to a Grid structure.
Returns
This function has no return value.
Note
This function does not work in spherical coordinates yet.

Compute source terms of the MHD equations in primitive variables. These include:

  • Geometrical sources;
  • Shearing-box terms
  • Gravity;
  • terms related to divergence of B control (Powell eight wave and GLM);
  • FARGO source terms.

The rationale for choosing during which sweep a particular source term has to be incorporated should match the same criterion used during the conservative update. For instance, in polar or cylindrical coordinates, curvilinear source terms are included during the radial sweep only.

Parameters
[in]statepointer to a State_1D structure
[in]beginitial index of computation
[in]endfinal index of computation
[in]a2array of sound speed
[in]harray of enthalpies (not needed in MHD)
[out]srcarray of source terms
[in]gridpointer to a Grid structure
Note
This function does not work in spherical coordinates yet. For future implementations we annotate hereafter the induction equation in spherical coordinates:

\[ \partial_tB_r + \frac{1}{r}\partial_\theta E_\phi - \frac{1}{r\sin\theta}\partial_\phi E_\theta = -E_\phi\cot\theta/r \]

\[ \partial_t B_\theta + \frac{1}{r\sin\theta}\partial_\phi E_r - \partial_rE_\phi = E_\phi/r \]

\[ \partial_t B_\phi + \partial_r E_\theta - \frac{1}{r}\partial_\theta E_r = - E_\theta/r\]

where

\[ E_\phi = -(v \times B)_\phi = - (v_r B_\theta - v_\theta B_r) \,,\qquad E_\theta = -(v \times B)_\theta = - (v_\phi B_r - v_r B_\phi) \]

Compute source terms of the RHD equations in primitive variables.

  • Geometrical sources;
  • Gravity;

The rationale for choosing during which sweep a particular source term has to be incorporated should match the same criterion used during the conservative update. For instance, in polar or cylindrical coordinates, curvilinear source terms are included during the radial sweep only.

Parameters
[in]statepointer to a State_1D structure;
[in]beginitial index of computation;
[in]endfinal index of computation;
[in]a2array of sound speed;
[in]harray of enthalpies (not needed in MHD);
[out]srcarray of source terms;
[in]gridpointer to a Grid structure.
Returns
This function has no return value.

Definition at line 71 of file prim_eqn.c.

97 {
98  int nv, i, j, k;
99  double tau, dA_dV, th;
100  double hscale; /* scale factor */
101  double *v, *vp, *A, *dV, r_inv, ct;
102  double *x1, *x2, *x3;
103  double *x1p, *x2p, *x3p;
104  double *dx1, *dx2, *dx3;
105  static double *phi_p;
106  double g[3], scrh;
107 
108 #if ROTATING_FRAME == YES
109  print1 ("! PrimSource(): does not work with rotations\n");
110  QUIT_PLUTO(1);
111 #endif
112 
113 /* ----------------------------------------------------------
114  1. Memory allocation and pointer shortcuts
115  ---------------------------------------------------------- */
116 
117  if (phi_p == NULL) phi_p = ARRAY_1D(NMAX_POINT, double);
118 
119  #if GEOMETRY == CYLINDRICAL
120  x1 = grid[IDIR].xgc; x1p = grid[IDIR].xr; dx1 = grid[IDIR].dx;
121  x2 = grid[JDIR].xgc; x2p = grid[JDIR].xr; dx2 = grid[JDIR].dx;
122  x3 = grid[KDIR].xgc; x3p = grid[KDIR].xr; dx3 = grid[KDIR].dx;
123  #else
124  x1 = grid[IDIR].x; x1p = grid[IDIR].xr; dx1 = grid[IDIR].dx;
125  x2 = grid[JDIR].x; x2p = grid[JDIR].xr; dx2 = grid[JDIR].dx;
126  x3 = grid[KDIR].x; x3p = grid[KDIR].xr; dx3 = grid[KDIR].dx;
127  #endif
128 
129  A = grid[g_dir].A;
130  dV = grid[g_dir].dV;
131  hscale = 1.0;
132 
133  i = g_i; j = g_j; k = g_k;
134 
135 /* ----------------------------------------------------------
136  initialize all elements of src to zero
137  ---------------------------------------------------------- */
138 
139  memset((void *)src[0], '\0',NMAX_POINT*NVAR*sizeof(double));
140 
141 /* ----------------------------------------------------------
142  2. Compute geometrical source terms
143  ---------------------------------------------------------- */
144 
145 #if GEOMETRY == CYLINDRICAL
146 
147  if (g_dir == IDIR) {
148  for (i = beg; i <= end; i++){
149  v = state->v[i];
150  tau = 1.0/v[RHO];
151  dA_dV = 1.0/x1[i];
152 
153  src[i][RHO] = -v[RHO]*v[VXn]*dA_dV;
154  #if COMPONENTS == 3
155  src[i][iVR] = v[iVPHI]*v[iVPHI]*dA_dV;
156  src[i][iVPHI] = -v[iVR]*v[iVPHI]*dA_dV;
157  #endif
158  #if EOS == IDEAL
159  src[i][PRS] = a2[i]*src[i][RHO];
160  #endif
161  }
162  }
163 
164 #elif GEOMETRY == POLAR
165 
166  if (g_dir == IDIR) {
167  for (i = beg; i <= end; i++){
168  v = state->v[i];
169 
170  tau = 1.0/v[RHO];
171  dA_dV = 1.0/x1[i];
172  src[i][RHO] = -v[RHO]*v[VXn]*dA_dV;
173 
174  #if COMPONENTS >= 2
175  src[i][iVR] = v[iVPHI]*v[iVPHI]*dA_dV;
176  #endif
177 
178  /* -- in 1D, all sources should be included during this sweep -- */
179 
180  #if DIMENSIONS == 1 && COMPONENTS >= 2
181  src[i][iVPHI] = -v[iVR]*v[iVPHI]*dA_dV;
182  #endif
183 
184  #if EOS == IDEAL
185  src[i][PRS] = a2[i]*src[i][RHO];
186  #endif
187  }
188  } else if (g_dir == JDIR) {
189  dA_dV = 1.0/x1[i];
190  hscale = x1[i];
191  for (j = beg; j <= end; j++){
192  v = state->v[j];
193  tau = 1.0/v[RHO];
194  src[j][iVPHI] = -v[iVR]*v[iVPHI]*dA_dV;
195  }
196  }
197 
198 #elif GEOMETRY == SPHERICAL
199 
200  print1 ("! PrimSource: not implemented in Spherical geometry\n");
201  QUIT_PLUTO(1);
202 
203 #endif
204 
205 /* ----------------------------------------------------------
206  3. Add body forces. This includes:
207  - Coriolis terms for the shearing box module
208  - Body forces
209  ---------------------------------------------------------- */
210 
211 #ifdef SHEARINGBOX
212 
213  if (g_dir == IDIR){
214  for (i = beg; i <= end; i++) {
215  src[i][VX1] = 2.0*state->v[i][VX2]*SB_OMEGA;
216  }
217  }else if (g_dir == JDIR){
218  for (j = beg; j <= end; j++) {
219  #ifdef FARGO
220  src[j][VX2] = (SB_Q - 2.0)*state->v[j][VX1]*SB_OMEGA;
221  #else
222  src[j][VX2] = -2.0*state->v[j][VX1]*SB_OMEGA;
223  #endif
224  }
225  }
226 
227 #endif
228 
229  #if (BODY_FORCE != NO)
230  if (g_dir == IDIR) {
231  i = beg-1;
232  #if BODY_FORCE & POTENTIAL
233  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
234  #endif
235  for (i = beg; i <= end; i++){
236  #if BODY_FORCE & VECTOR
237  v = state->v[i];
238  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
239  src[i][VX1] += g[IDIR];
240  #endif
241  #if BODY_FORCE & POTENTIAL
242  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
243  src[i][VX1] -= (phi_p[i] - phi_p[i-1])/(hscale*dx1[i]);
244  #endif
245 
246  /* -- Add tangential components in 1D -- */
247 
248  #if DIMENSIONS == 1
249  EXPAND( ,
250  src[i][VX2] += g[JDIR]; ,
251  src[i][VX3] += g[KDIR];)
252  #endif
253  }
254  }else if (g_dir == JDIR){
255  j = beg - 1;
256  #if BODY_FORCE & POTENTIAL
257  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
258  #endif
259  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
260  hscale = x1[i];
261  #endif
262  for (j = beg; j <= end; j++){
263  #if BODY_FORCE & VECTOR
264  v = state->v[j];
265  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
266  src[j][VX2] += g[JDIR];
267  #endif
268  #if BODY_FORCE & POTENTIAL
269  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
270  src[j][VX2] -= (phi_p[j] - phi_p[j-1])/(hscale*dx2[j]);
271  #endif
272 
273  /* -- Add 3rd component in 2D -- */
274 
275  #if DIMENSIONS == 2 && COMPONENTS == 3
276  src[j][VX3] += g[KDIR];
277  #endif
278  }
279  }else if (g_dir == KDIR){
280  k = beg - 1;
281  #if BODY_FORCE & POTENTIAL
282  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
283  #endif
284  #if GEOMETRY == SPHERICAL
285  th = x2[j];
286  hscale = x1[i]*sin(th);
287  #endif
288  for (k = beg; k <= end; k++){
289  #if BODY_FORCE & VECTOR
290  v = state->v[k];
291  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
292  src[k][VX3] += g[KDIR];
293  #endif
294  #if BODY_FORCE & POTENTIAL
295  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
296  src[k][VX3] -= (phi_p[k] - phi_p[k-1])/(hscale*dx3[k]);
297  #endif
298  }
299  }
300  #endif
301 
302 /* ---------------------------------------------------------------
303  5. Add FARGO source terms (except for SHEARINGBOX).
304  When SHEARINGBOX is also enabled, we do not include
305  these source terms since they're all provided by body_force)
306  --------------------------------------------------------------- */
307 
308  #if (defined FARGO) && !(defined SHEARINGBOX)
309  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
310  print1 ("! Time Stepping works only in Cartesian or cylindrical coords\n");
311  print1 ("! Use RK instead\n");
312  QUIT_PLUTO(1);
313  #endif
314 
315  double **wA, *dx, *dz;
316  wA = FARGO_GetVelocity();
317  if (g_dir == IDIR){
318  dx = grid[IDIR].dx;
319  for (i = beg; i <= end; i++){
320  v = state->v[i];
321  src[i][VX2] -= 0.5*v[VX1]*(wA[k][i+1] - wA[k][i-1])/dx[i];
322  }
323  }else if (g_dir == KDIR){
324  dz = grid[KDIR].dx;
325  for (k = beg; k <= end; k++){
326  v = state->v[k];
327  src[k][VX2] -= 0.5*v[VX3]*(wA[k+1][i] - wA[k-1][i])/dz[k];
328  }
329  }
330  #endif
331 }
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
VX2
#define VX2
Definition: mod_defs.h:29
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
GRID::xr
double * xr
Definition: structs.h:81
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
wA
static double ** wA
Definition: fargo_velocity.c:24
iVPHI
#define iVPHI
Definition: mod_defs.h:67
func_param::v
double v[NVAR]
Definition: eos.h:106
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GRID::dx
double * dx
Definition: structs.h:83
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
iVR
#define iVR
Definition: mod_defs.h:65
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
SB_OMEGA
#define SB_OMEGA
Disk local orbital frequency .
Definition: shearingbox.h:81
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
GRID::xgc
double * xgc
Cell volumetric centroid (!= x when geometry != CARTESIAN).
Definition: structs.h:84
BodyForcePotential
double BodyForcePotential(double, double, double)
Definition: init.c:479
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
VXn
int VXn
Definition: globals.h:73
SB_Q
#define SB_Q
The shear parameter, .
Definition: shearingbox.h:76
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
FARGO_GetVelocity
double ** FARGO_GetVelocity(void)
Definition: fargo_velocity.c:270
BodyForceVector
void BodyForceVector(double *, double *, double, double, double)
Definition: init.c:441
VX3
#define VX3
Definition: mod_defs.h:30
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the caller graph for this function:

void PrimToCons ( double **  uprim,
double **  ucons,
int  beg,
int  end 
)

Convert primitive variables to conservative variables.

Parameters
[in]uprimarray of primitive variables
[out]uconsarray of conservative variables
[in]begstarting index of computation
[in]endfinal index of computation

Convert primitive variables in conservative variables.

Parameters
[in]uprimarray of primitive variables
[out]uconsarray of conservative variables
[in]begstarting index of computation
[in]endfinal index of computation

Definition at line 26 of file mappers.c.

36 {
37  int i, nv, status;
38  double *v, *u;
39  double rho, rhoe, T, gmm1;
40 
41  #if EOS == IDEAL
42  gmm1 = g_gamma - 1.0;
43  #endif
44  for (i = ibeg; i <= iend; i++) {
45 
46  v = uprim[i];
47  u = ucons[i];
48 
49  u[RHO] = rho = v[RHO];
50  EXPAND (u[MX1] = rho*v[VX1]; ,
51  u[MX2] = rho*v[VX2]; ,
52  u[MX3] = rho*v[VX3];)
53 
54  #if EOS == IDEAL
55  u[ENG] = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
56  u[ENG] = 0.5*rho*u[ENG] + v[PRS]/gmm1;
57 
58  #elif EOS == PVTE_LAW
59  status = GetPV_Temperature(v, &T);
60  if (status != 0){
61  T = T_CUT_RHOE;
62  v[PRS] = Pressure(v, T);
63  }
64  rhoe = InternalEnergy(v, T);
65 
66  u[ENG] = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
67  u[ENG] = 0.5*rho*u[ENG] + rhoe;
68 
69  if (u[ENG] != u[ENG]){
70  print("! PrimToCons(): E = NaN, rhoe = %8.3e, T = %8.3e\n",rhoe, T);
71  QUIT_PLUTO(1);
72  }
73  #endif
74 
75 #if DUST == YES
76  u[RHO_D] = v[RHO_D];
77  EXPAND(u[MX1_D] = v[RHO_D]*v[VX1_D]; ,
78  u[MX2_D] = v[RHO_D]*v[VX2_D]; ,
79  u[MX3_D] = v[RHO_D]*v[VX3_D];)
80 #endif
81 
82 #if NSCL > 0
83  NSCL_LOOP(nv) u[nv] = rho*v[nv];
84 #endif
85  }
86 
87 }
IDEAL
#define IDEAL
Definition: pluto.h:45
MX3
#define MX3
Definition: mod_defs.h:22
EOS
#define EOS
Definition: pluto.h:341
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
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
T_CUT_RHOE
#define T_CUT_RHOE
Sets the lowest cut-off temperature used in the PVTE_LAW equation of state.
Definition: eos.h:69
InternalEnergy
double InternalEnergy(double *v, double T)
Definition: internal_energy.c:255
RHO
#define RHO
Definition: mod_defs.h:19
func_param::rhoe
double rhoe
Definition: eos.h:107
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
NSCL
#define NSCL
Definition: pluto.h:572
MX2
#define MX2
Definition: mod_defs.h:21
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
i
int i
Definition: analysis.c:2
Pressure
double Pressure(double *v, double T)
Definition: thermal_eos.c:179
VX3
#define VX3
Definition: mod_defs.h:30
GetPV_Temperature
int GetPV_Temperature(double *v, double *T)
Definition: thermal_eos.c:221
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Variable Documentation

Riemann_Solver AUSMp_Solver

Definition at line 108 of file mod_defs.h.

Riemann_Solver HLL_Solver

Definition at line 106 of file mod_defs.h.

Riemann_Solver HLLC_Solver

Definition at line 106 of file mod_defs.h.

Riemann_Solver LF_Solver

Definition at line 106 of file mod_defs.h.

Riemann_Solver Roe_Solver

Definition at line 106 of file mod_defs.h.

Riemann_Solver RusanovDW_Solver

Definition at line 106 of file mod_defs.h.

Riemann_Solver TwoShock_Solver

Definition at line 106 of file mod_defs.h.