PLUTO
Classes | Macros | Typedefs | Functions | Variables
mod_defs.h File Reference

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

Go to the source code of this file.

Classes

struct  MAP_PARAM
 

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   (2 + COMPONENTS)
 
#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
 

Typedefs

typedef struct MAP_PARAM Map_param
 

Functions

int ConsToPrim (double **, double **, int, int, unsigned char *)
 
void PrimEigenvectors (double *, double, double, double *, double **, double **)
 
int EnergySolve (Map_param *)
 
int EntropySolve (Map_param *)
 
int PressureFix (Map_param *)
 
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 *)
 
void VelocityLimiter (double *, double *, double *)
 
void ConvertTo4vel (double **, int, int)
 
void ConvertTo3vel (double **, int, int)
 

Variables

Riemann_Solver TwoShock_Solver
 
Riemann_Solver LF_Solver
 
Riemann_Solver HLL_Solver
 
Riemann_Solver HLLC_Solver
 

Detailed Description

Set labels, indexes and prototypes for the RHD module.

Contains variable names and prototypes for the RHD 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   (2 + COMPONENTS)

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.

Typedef Documentation

typedef struct MAP_PARAM Map_param

The Map_param structure is used to pass input/output arguments during the conversion from conservative to primitive variables operated by the ConsToPrim() function. The output parameter, rho, W, lor and p, must be set at the end of every root-finder routine (EnergySolve(), EntropySolve() and PressureFix()). Additionally, some of the input parameters must be re-computed in EntropySolve() and PressureFix().

Function Documentation

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 ConvertTo3vel ( double **  ,
int  ,
int   
)

Definition at line 25 of file four_vel.c.

30 {
31  int i;
32  double lor;
33 
34  for (i = beg; i <= end; i++){
35  lor = EXPAND( v[i][VX1]*v[i][VX1],
36  + v[i][VX2]*v[i][VX2],
37  + v[i][VX3]*v[i][VX3]);
38  lor = sqrt(1.0 + lor);
39 
40  EXPAND(v[i][VX1] /= lor; ,
41  v[i][VX2] /= lor; ,
42  v[i][VX3] /= lor;)
43  }
44 }
VX2
#define VX2
Definition: mod_defs.h:29
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
i
int i
Definition: analysis.c:2
VX3
#define VX3
Definition: mod_defs.h:30

Here is the caller graph for this function:

void ConvertTo4vel ( double **  ,
int  ,
int   
)

Definition at line 4 of file four_vel.c.

9 {
10  int i;
11  double lor;
12 
13  for (i = beg; i <= end; i++){
14  lor = EXPAND( v[i][VX1]*v[i][VX1],
15  + v[i][VX2]*v[i][VX2],
16  + v[i][VX3]*v[i][VX3]);
17  lor = 1.0/sqrt(1.0 - lor);
18 
19  EXPAND(v[i][VX1] *= lor; ,
20  v[i][VX2] *= lor; ,
21  v[i][VX3] *= lor;)
22  }
23 }
VX2
#define VX2
Definition: mod_defs.h:29
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
i
int i
Definition: analysis.c:2
VX3
#define VX3
Definition: mod_defs.h:30

Here is the caller graph for this function:

int EnergySolve ( Map_param *  par)

Use Newton algorithm to recover pressure p0 from conservative quantities D, m and E

Solve f(W) = 0, where f(W) is Equation (A4) or (A6).

Returns
Return (0) is successful, (1) otherwise.

Definition at line 23 of file energy_solve.c.

29 {
30  int iter;
31  double p, D_1, alpha, alpha2, lor2, lor, m, tol=1.e-14;
32  double tau, theta, h, dh_dp, dh_dtau, gmmr;
33  double yp, dyp, dp, scrh;
34  double D, E, m2;
35 
36  D = par->D;
37  E = par->E;
38  m2 = par->m2;
39 
40  #if EOS == IDEAL
41  gmmr = g_gamma/(g_gamma - 1.0);
42  #endif
43  m = sqrt(m2);
44  p = m - E;
45  p = MAX(p, 1.e-18);
46  D_1 = 1.0/D;
47  for (iter = 0; iter < MAX_ITER; iter++) {
48 
49  alpha = E + p;
50  alpha2 = alpha*alpha;
51  lor2 = 1.0 - m2/alpha2;
52 
53  lor2 = MAX(lor2,1.e-9);
54  lor2 = 1.0/lor2;
55  lor = sqrt(lor2);
56 
57  tau = lor*D_1;
58  theta = p*tau;
59 
60  #if EOS == IDEAL
61  h = 1.0 + gmmr*theta;
62  dh_dp = gmmr*tau;
63  dh_dtau = gmmr*p;
64  #elif EOS == TAUB
65  h = 2.5*theta + sqrt(2.25*theta*theta + 1.0);
66  scrh = (5.0*h - 8.0*theta)/(2.0*h - 5.0*theta);
67  dh_dp = tau*scrh;
68  dh_dtau = p*scrh;
69  #endif
70 
71  yp = D*h*lor - E - p;
72  dyp = D*lor*dh_dp - m2*lor2*lor/(alpha2*alpha)*(lor*dh_dtau + D*h) - 1.0;
73  dp = yp/dyp;
74  p -= dp;
75  if (fabs (dp) < tol*p) break;
76  }
77 
78 /* ----------------------------------------------------------
79  check if solution is consistent
80  ---------------------------------------------------------- */
81 
82  if (p != p) {
83  WARNING(
84  print("! EnergySolve: NaN found while recovering pressure, ");
85  )
86  return 1;
87  }
88 
89  if (p < 0.0){
90  WARNING(
91  print("! EnergySolve: negative pressure (%8.2e), ", p);
92  )
93  return 1;
94  }
95 
96  par->rho = par->D/lor;
97  par->W = D*h*lor;
98  par->prs = p;
99  par->lor = lor;
100 
101 /* -- Recompute entropy consistently -- */
102 
103 #if ENTROPY_SWITCH
104  #if EOS == IDEAL
105  par->sigma_c = p*lor/pow(par->rho,g_gamma-1);
106  #elif EOS == TAUB
107  theta = p/rho;
108  par->sigma_c = p*lor/pow(par->rho,2.0/3.0)*(1.5*theta + sqrt(2.25*theta*theta + 1.0);
109  #endif
110 #endif
111 
112  return 0; /* -- success -- */
113 }
alpha
static double alpha
Definition: init.c:3
MAX
#define MAX(a, b)
Definition: macros.h:101
g_gamma
double g_gamma
Definition: globals.h:112
MAP_PARAM::m2
double m2
Square of total momentum (input).
Definition: mod_defs.h:89
MAP_PARAM::D
double D
Lab density (input).
Definition: mod_defs.h:86
configure.scrh
tuple scrh
Definition: configure.py:200
WARNING
#define WARNING(a)
Definition: macros.h:217
gmmr
static double gmmr
Definition: two_shock.c:42
flow_past_cyc.D
tuple D
Definition: flow_past_cyc.py:11
MAP_PARAM::rho
double rho
proper density (output)
Definition: mod_defs.h:91
MAP_PARAM::W
double W
D*h*lor (output).
Definition: mod_defs.h:92
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MAP_PARAM::lor
double lor
Lorentz factor (output).
Definition: mod_defs.h:93
MAP_PARAM::prs
double prs
Thermal pressure (output).
Definition: mod_defs.h:94
MAX_ITER
#define MAX_ITER
Definition: energy_solve.c:21
MAP_PARAM::sigma_c
double sigma_c
Conserved entropy (input).
Definition: mod_defs.h:87
MAP_PARAM::E
double E
Total energy (input).
Definition: mod_defs.h:88
int EntropySolve ( Map_param *  par)

Convert the conservative variables {D, m, sigma_c} (where sigma_c = D*sigma is the conserved entropy) to primitive variable using a Newton-Raphson/Bisection scheme.

LAST MODIFIED

February 5th (2012) by C. Zanni & A. Mignone (zanni.nosp@m.@oat.nosp@m.o.ina.nosp@m.f.it, migno.nosp@m.ne@o.nosp@m.ato.i.nosp@m.naf..nosp@m.it)

Definition at line 21 of file entropy_solve.c.

33 {
34  int k, done;
35  double v2, v2max, v2min, acc = 1.e-13;
36  double h, W, W2, W3;
37  double scrh, x, temp;
38  double lor, lor2, rho, sigma, th;
39  double dW_dlor;
40  double fv2, dfv2, dv2, dv2old;
41  double m2, D;
42 
43  D = par->D;
44  m2 = par->m2;
45 
46  sigma = par->sigma_c/D;
47  v2max = 1.0 - 1.e-8;
48  v2min = 0.0;
49 
50  done = 0;
51 
52  v2 = 0.5*(v2max+v2min);
53  dv2old = v2max-v2min;
54  dv2 = dv2old;
55 
56  lor2 = 1.0/(1.0 - v2);
57  lor = sqrt(lor2);
58  rho = D/lor;
59 
60  #if EOS == IDEAL
61  x = pow(rho,g_gamma - 1.0);
62  h = 1.0 + g_gamma/(g_gamma - 1.0)*sigma*x;
63  W = D*h*lor;
64  dW_dlor = (2 - g_gamma)*W/lor + (g_gamma - 1.0)*D;
65  #elif EOS == TAUB
66  x = pow(rho,2.0/3.0);
67  th = sigma*x/sqrt(1.0 + 3.0*sigma*x);
68  h = 2.5*th + sqrt(2.25*th*th + 1.0);
69  W = D*h*lor;
70 
71  scrh = -2.0*x/(3.0*lor)*sigma*(2.0*sigma*x - 3.0*th*th);
72  scrh /= 2.0*th*(1.0 + 3.0*sigma*x);
73  dW_dlor = W/lor + D*lor*(5.0*h - 8.0*th)/(2.0*h - 5.0*th)*scrh;
74  #endif
75 
76  W2 = W*W;
77  W3 = W2*W;
78 
79  fv2 = v2 - m2/W2;
80  dfv2 = 1.0 + m2/W3*dW_dlor*lor2*lor;
81 
82  for (k = 1; k < MAX_ITER; k++) {
83  if ((((v2-v2max)*dfv2-fv2)*((v2-v2min)*dfv2-fv2) > 0.0)
84  || (fabs(2.*fv2) > fabs(dv2old*dfv2))) {
85  dv2old = dv2;
86  dv2 = 0.5*(v2max-v2min);
87  v2 = v2min+dv2;
88  if (v2min == v2) done = 1;
89  } else {
90  dv2old = dv2;
91  dv2 = fv2/dfv2;
92  temp = v2;
93  v2 -= dv2;
94  if (temp == v2) done = 1;
95  }
96 
97  if (fabs(dv2) < acc*v2 || fabs(fv2) < acc) done = 1;
98 
99  lor2 = 1.0/(1.0 - v2);
100  lor = sqrt(lor2);
101 
102  rho = D/lor;
103 
104  #if EOS == IDEAL
105  x = pow(rho,g_gamma - 1.0);
106  h = 1.0 + g_gamma/(g_gamma - 1.0)*sigma*x;
107  W = D*h*lor;
108  dW_dlor = (2 - g_gamma)*W/lor + (g_gamma - 1.0)*D;
109  #elif EOS == TAUB
110  x = pow(rho,2.0/3.0);
111  th = sigma*x/sqrt(1.0 + 3.0*sigma*x);
112  h = 2.5*th + sqrt(2.25*th*th + 1.0);
113  W = D*h*lor;
114 
115  scrh = -2.0*x/(3.0*lor)*sigma*(2.0*sigma*x - 3.0*th*th);
116  scrh /= 2.0*th*(1.0 + 3.0*sigma*x);
117  dW_dlor = W/lor + D*lor*(5.0*h-8.0*th)/(2.0*h-5.0*th)*scrh;
118  #endif
119 
120  W2 = W*W;
121 
122  if (done) break;
123 
124  W3 = W2*W;
125  fv2 = v2 - m2/W2;
126 
127  if (fv2 < 0.0) v2min = v2;
128  else v2max = v2;
129 
130  dfv2 = 1.0 + m2/W3*dW_dlor*lor2*lor;
131  }
132 
133  #if EOS == IDEAL
134  par->prs = sigma*x*rho;
135  #elif EOS == TAUB
136  par->prs = th*rho;
137  #endif
138 
139  if (k == MAX_ITER) {
140  WARNING(
141  print ("! EntropySolve: too many iterations,%d, ",k);
142  )
143  par->prs = g_smallPressure;
144  par->E = W - par->prs;
145  return(1);
146  }
147 
148  if (par->prs < 0.0 || sigma < 0.0 || W < 0.0) {
149  WARNING(
150  print ("! EntropySolve: negative pressure, p = %12.6e\n",par->prs);
151  )
152  par->prs = g_smallPressure;
153  par->E = W - par->prs;
154  return(1);
155  }
156 
157  par->rho = par->D/lor;
158  par->lor = lor;
159  par->E = W - par->prs; /* redefine energy */
160  return(0); /* -- success -- */
161 }
g_gamma
double g_gamma
Definition: globals.h:112
MAP_PARAM::m2
double m2
Square of total momentum (input).
Definition: mod_defs.h:89
MAP_PARAM::D
double D
Lab density (input).
Definition: mod_defs.h:86
configure.scrh
tuple scrh
Definition: configure.py:200
WARNING
#define WARNING(a)
Definition: macros.h:217
flow_past_cyc.D
tuple D
Definition: flow_past_cyc.py:11
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MAX_ITER
#define MAX_ITER
Definition: entropy_solve.c:19
MAP_PARAM::prs
double prs
Thermal pressure (output).
Definition: mod_defs.h:94
MAP_PARAM::sigma_c
double sigma_c
Conserved entropy (input).
Definition: mod_defs.h:87
MAP_PARAM::E
double E
Total energy (input).
Definition: mod_defs.h:88
void Flux ( double **  ucons,
double **  uprim,
double *  h,
double **  fx,
real *  pr,
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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the caller graph for this function:

int PressureFix ( Map_param *  par)

Fix p to a small value, solve for the square of velocity by using secant algorithm applied to Eq (9) of Mignone, Plewa & Bodo (2005). This step involved re-computing W at each step of the iteration. Once the root has been found, we recompute total energy E. Return 0 if succesful, 1 otherwise.

Definition at line 22 of file pressure_fix.c.

31 {
32  int k, done=0;
33  double v2, v2c, fc, f, dW, S2_W2;
34  double fmin, fmax, v2min, v2max;
35 
36  par->prs = g_smallPressure;
37 
38  v2max = 1.0-1.e-8;
39  v2c = 0.95;
40  fc = VelocitySquareFunc(v2c, par);
41  v2 = 0.96;
42  for (k = 1; k < MAX_ITER; k++){
43  f = VelocitySquareFunc(v2, par);
44  if (done == 1) break;
45  dW = (v2 - v2c)/(f - fc)*f;
46  v2c = v2; fc = f;
47  v2 -= dW;
48  v2 = MIN(v2max,v2);
49  v2 = MAX(v2, 0.0);
50  if (fabs(f) < 1.e-9) done = 1;
51  }
52  if (v2 >= 1.0 || k >= MAX_ITER) {
53  print ("! PressureFix: too many iter while fixing p , v^2 = %f\n", v2);
54  return (1);
55  }
56 
57 /* -- Redefine energy, density and entropy -- */
58 
59  par->E = par->W - par->prs;
60  par->rho = par->D/par->lor;
61 
62 #if ENTROPY_SWITCH
63 {
64  double rho = par->rho;
65  double th = par->prs/rho;
66  #if EOS == IDEAL
67  par->sigma_c = par->prs*par->lor/pow(rho,g_gamma-1);
68  #elif EOS == TAUB
69  th = par->prs/rho;
70  par->sigma_c = par->prs*par->lor/pow(rho,2.0/3.0)*(1.5*th + sqrt(2.25*th*th + 1.0);
71  #endif
72 }
73 #endif
74 
75  return(0); /* -- success -- */
76 }
MAX
#define MAX(a, b)
Definition: macros.h:101
g_gamma
double g_gamma
Definition: globals.h:112
MAP_PARAM::D
double D
Lab density (input).
Definition: mod_defs.h:86
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
MIN
#define MIN(a, b)
Definition: macros.h:104
MAP_PARAM::rho
double rho
proper density (output)
Definition: mod_defs.h:91
MAP_PARAM::W
double W
D*h*lor (output).
Definition: mod_defs.h:92
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
VelocitySquareFunc
static double VelocitySquareFunc(double, Map_param *)
Definition: pressure_fix.c:79
MAP_PARAM::lor
double lor
Lorentz factor (output).
Definition: mod_defs.h:93
MAP_PARAM::prs
double prs
Thermal pressure (output).
Definition: mod_defs.h:94
MAP_PARAM::sigma_c
double sigma_c
Conserved entropy (input).
Definition: mod_defs.h:87
MAX_ITER
#define MAX_ITER
Definition: pressure_fix.c:19
MAP_PARAM::E
double E
Total energy (input).
Definition: mod_defs.h:88

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.

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

Here is the call graph for this function:

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 call 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
void VelocityLimiter ( double *  v,
double *  vp,
double *  vm 
)

Check whether the total reconstructed velocity is > 1 If a superluminal value occurs, flatten distribution.

Definition at line 16 of file vel_limiter.c.

22 {
23 #if RECONSTRUCT_4VEL == NO && COMPONENTS > 1
24  int nv;
25  double v2m, v2p;
26 
27  v2m = EXPAND(vm[VX1]*vm[VX1], + vm[VX2]*vm[VX2], + vm[VX3]*vm[VX3]);
28  v2p = EXPAND(vp[VX1]*vp[VX1], + vp[VX2]*vp[VX2], + vp[VX3]*vp[VX3]);
29  if (v2m >= 1.0 || v2p >= 1.0){
30  for (nv = NVAR; nv--; ) vm[nv] = vp[nv] = v[nv];
31  }
32 #endif
33 }
VX2
#define VX2
Definition: mod_defs.h:29
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
VX3
#define VX3
Definition: mod_defs.h:30
NVAR
#define NVAR
Definition: pluto.h:609

Here is the caller graph for this function:

Variable Documentation

Riemann_Solver HLL_Solver

Definition at line 119 of file mod_defs.h.

Riemann_Solver HLLC_Solver

Definition at line 119 of file mod_defs.h.

Riemann_Solver LF_Solver

Definition at line 119 of file mod_defs.h.

Riemann_Solver TwoShock_Solver

Definition at line 119 of file mod_defs.h.