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

Module header file for relativistic MHD (RMHD). More...

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

struct  MAP_PARAM
 

Macros

#define RESISTIVE_RMHD   NO
 
#define RHO   0
 
#define MX1   1
 
#define MX2   (COMPONENTS >= 2 ? 2: 255)
 
#define MX3   (COMPONENTS == 3 ? 3: 255)
 
#define BX1   (COMPONENTS + 1)
 
#define BX2   (COMPONENTS >= 2 ? (BX1+1): 255)
 
#define BX3   (COMPONENTS == 3 ? (BX1+2): 255)
 
#define PSI_GLM   (2*COMPONENTS + 1 + HAVE_ENERGY)
 
#define VX1   MX1
 
#define VX2   MX2
 
#define VX3   MX3
 
#define NFLX   (2 + 2*COMPONENTS + (DIVB_CONTROL == DIV_CLEANING))
 
#define AX1   (NVAR + 1)
 
#define AX2   (NVAR + 2)
 
#define AX3   (NVAR + 3)
 
#define AX   AX1 /* backward compatibility */
 
#define AY   AX2
 
#define AZ   AX3
 
#define iVR   VX1
 
#define iVZ   VX2
 
#define iVPHI   VX3
 
#define iMR   MX1
 
#define iMZ   MX2
 
#define iMPHI   MX3
 
#define iBR   BX1
 
#define iBZ   BX2
 
#define iBPHI   BX3
 
#define iVR   VX1
 
#define iVPHI   VX2
 
#define iVZ   VX3
 
#define iMR   MX1
 
#define iMPHI   MX2
 
#define iMZ   MX3
 
#define iBR   BX1
 
#define iBPHI   BX2
 
#define iBZ   BX3
 
#define iVR   VX1
 
#define iVTH   VX2
 
#define iVPHI   VX3
 
#define iMR   MX1
 
#define iMTH   MX2
 
#define iMPHI   MX3
 
#define iBR   BX1
 
#define iBTH   BX2
 
#define iBPHI   BX3
 
#define RMHD_FAST_EIGENVALUES   NO
 If set to YES, use approximate (and faster) expressions when computing the fast magnetosonic speed, see Sect. More...
 
#define RMHD_REDUCED_ENERGY   YES
 By turning RMHD_REDUCED_ENERGY to YES, we let PLUTO evolve the total energy minus the mass density contribution. More...
 

Typedefs

typedef struct MAP_PARAM Map_param
 

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 PRIM_EIGENVECTORS (double *, double, double, double *, double **, double **)
 
int Eigenvalues (double *, double, double, double *)
 
int EntropySolve (Map_param *)
 
int EnergySolve (Map_param *)
 
int PressureFix (Map_param *)
 
void Flux (double **, double **, double *, double **, double *, int, int)
 
void HLL_Speed (double **, double **, double *, double *, double *, double *, double *, double *, int, int)
 
int MaxSignalSpeed (double **, double *, double *, double *, double *, int, int)
 
void PrimToCons (double **, double **, int, int)
 
void VelocityLimiter (double *, double *, double *)
 
int Magnetosonic (double *vp, double cs2, double h, double *lambda)
 
int QuarticSolve (double, double, double, double, double *)
 
int CubicSolve (double, double, double, double *)
 
void POWELL_DIVB_SOURCE (const State_1D *, int, int, Grid *)
 
void HLL_DIVB_SOURCE (const State_1D *, double **, int, int, Grid *)
 

Variables

Riemann_Solver LF_Solver
 
Riemann_Solver HLL_Solver
 
Riemann_Solver HLLC_Solver
 
Riemann_Solver HLLD_Solver
 
Riemann_Solver HLL_Linde_Solver
 

Detailed Description

Module header file for relativistic MHD (RMHD).

Set label, indexes and basic prototyping for the relativistic MHD module.

Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
C. Zanni (zanni.nosp@m.@oat.nosp@m.o.ina.nosp@m.f.it)
Date
April 2, 2015

Definition in file mod_defs.h.

Macro Definition Documentation

#define AX   AX1 /* backward compatibility */

Definition at line 115 of file mod_defs.h.

#define AX1   (NVAR + 1)

Definition at line 111 of file mod_defs.h.

#define AX2   (NVAR + 2)

Definition at line 112 of file mod_defs.h.

#define AX3   (NVAR + 3)

Definition at line 113 of file mod_defs.h.

#define AY   AX2

Definition at line 116 of file mod_defs.h.

#define AZ   AX3

Definition at line 117 of file mod_defs.h.

#define BX1   (COMPONENTS + 1)

Definition at line 29 of file mod_defs.h.

#define BX2   (COMPONENTS >= 2 ? (BX1+1): 255)

Definition at line 30 of file mod_defs.h.

#define BX3   (COMPONENTS == 3 ? (BX1+2): 255)

Definition at line 31 of file mod_defs.h.

#define iBPHI   BX3

Definition at line 168 of file mod_defs.h.

#define iBPHI   BX2

Definition at line 168 of file mod_defs.h.

#define iBPHI   BX3

Definition at line 168 of file mod_defs.h.

#define iBR   BX1

Definition at line 166 of file mod_defs.h.

#define iBR   BX1

Definition at line 166 of file mod_defs.h.

#define iBR   BX1

Definition at line 166 of file mod_defs.h.

#define iBTH   BX2

Definition at line 167 of file mod_defs.h.

#define iBZ   BX2

Definition at line 152 of file mod_defs.h.

#define iBZ   BX3

Definition at line 152 of file mod_defs.h.

#define iMPHI   MX3

Definition at line 164 of file mod_defs.h.

#define iMPHI   MX2

Definition at line 164 of file mod_defs.h.

#define iMPHI   MX3

Definition at line 164 of file mod_defs.h.

#define iMR   MX1

Definition at line 162 of file mod_defs.h.

#define iMR   MX1

Definition at line 162 of file mod_defs.h.

#define iMR   MX1

Definition at line 162 of file mod_defs.h.

#define iMTH   MX2

Definition at line 163 of file mod_defs.h.

#define iMZ   MX2

Definition at line 148 of file mod_defs.h.

#define iMZ   MX3

Definition at line 148 of file mod_defs.h.

#define iVPHI   VX3

Definition at line 160 of file mod_defs.h.

#define iVPHI   VX2

Definition at line 160 of file mod_defs.h.

#define iVPHI   VX3

Definition at line 160 of file mod_defs.h.

#define iVR   VX1

Definition at line 158 of file mod_defs.h.

#define iVR   VX1

Definition at line 158 of file mod_defs.h.

#define iVR   VX1

Definition at line 158 of file mod_defs.h.

#define iVTH   VX2

Definition at line 159 of file mod_defs.h.

#define iVZ   VX2

Definition at line 144 of file mod_defs.h.

#define iVZ   VX3

Definition at line 144 of file mod_defs.h.

#define MX1   1

Definition at line 26 of file mod_defs.h.

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

Definition at line 27 of file mod_defs.h.

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

Definition at line 28 of file mod_defs.h.

#define NFLX   (2 + 2*COMPONENTS + (DIVB_CONTROL == DIV_CLEANING))

Definition at line 45 of file mod_defs.h.

#define PSI_GLM   (2*COMPONENTS + 1 + HAVE_ENERGY)

Definition at line 38 of file mod_defs.h.

#define RESISTIVE_RMHD   NO

Definition at line 16 of file mod_defs.h.

#define RHO   0

Definition at line 25 of file mod_defs.h.

#define RMHD_FAST_EIGENVALUES   NO

If set to YES, use approximate (and faster) expressions when computing the fast magnetosonic speed, see Sect.

3.3 of Del Zanna, A&A (2006), 473. Solutions of quartic equation is avoided and replace with upper bounds provided by quadratic equation.

Definition at line 177 of file mod_defs.h.

#define RMHD_REDUCED_ENERGY   YES

By turning RMHD_REDUCED_ENERGY to YES, we let PLUTO evolve the total energy minus the mass density contribution.

Definition at line 193 of file mod_defs.h.

#define VX1   MX1

Definition at line 41 of file mod_defs.h.

#define VX2   MX2

Definition at line 42 of file mod_defs.h.

#define VX3   MX3

Definition at line 43 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 in the relativistic modules (RHD and RMHD). 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().

Enumeration Type Documentation

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

Definition at line 58 of file mod_defs.h.

58  {
59  KFASTM, KFASTP, KENTRP
60 
61  #if DIVB_CONTROL != DIV_CLEANING
62  , KDIVB
63  #endif
64 
65  #if COMPONENTS >= 2
66  , KSLOWM, KSLOWP
67  #if COMPONENTS == 3
68  , KALFVM, KALFVP
69  #endif
70  #endif
71 
72  #if DIVB_CONTROL == DIV_CLEANING
73  , KPSI_GLMM, KPSI_GLMP
74  #endif
75 };
KPSI_GLMP
Definition: mod_defs.h:75
KFASTP
Definition: mod_defs.h:58
KFASTM
Definition: mod_defs.h:58
KENTRP
Definition: mod_defs.h:59
KPSI_GLMM
Definition: mod_defs.h:75

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

Here is the call graph for this function:

Here is the caller graph for this function:

int CubicSolve ( double  ,
double  ,
double  ,
double *   
)
int Eigenvalues ( double *  v,
double  cs2,
double  h,
double *  lambda 
)

Compute an approximate expression for the fast magnetosonic speed using the upper-bound estimated outlined by Leismann et al. (A&A 2005, 436, 503), Eq. [27].

Definition at line 277 of file eigenv.c.

284 {
285  double vel2, lor2, Bmag2, b2, ca2, om2, vB2;
286  double vB, scrh, vl, vx, vx2;
287 
288  vel2 = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
289  vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
290  Bmag2 = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]);
291 
292  if (vel2 >= 1.0){
293  print ("! Eigenavalues(): |v|^2 = %f > 1\n",vel2);
294  return 1;
295  }
296 
297  vx = v[VXn];
298  vx2 = vx*vx;
299  vB2 = vB*vB;
300  b2 = Bmag2*(1.0 - vel2) + vB2;
301  ca2 = b2/(v[RHO]*h + b2);
302 
303  om2 = cs2 + ca2 - cs2*ca2;
304  vl = vx*(1.0 - om2)/(1.0 - vel2*om2);
305  scrh = om2*(1.0 - vel2)*((1.0 - vel2*om2) - vx2*(1.0 - om2));
306  scrh = sqrt(scrh)/(1.0 - vel2*om2);
307  lambda[KFASTM] = vl - scrh;
308  lambda[KFASTP] = vl + scrh;
309 
310  if (fabs(lambda[KFASTM])>1.0 || fabs(lambda[KFASTP]) > 1.0){
311  print ("! Eigenvalues(): vm, vp = %8.3e, %8.3e\n",lambda[KFASTM],lambda[KFASTP]);
312  QUIT_PLUTO(1);
313  }
314  if (lambda[KFASTM] != lambda[KFASTM] || lambda[KFASTP] != lambda[KFASTP]){
315  print ("! Eigenvalues(): nan, vm, vp = %8.3e, %8.3e\n",lambda[KFASTM],lambda[KFASTP]);
316  QUIT_PLUTO(1);
317  }
318 
319  scrh = fabs(vx)/sqrt(cs2);
320  g_maxMach = MAX(scrh, g_maxMach);
321 }
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
func_param::v
double v[NVAR]
Definition: eos.h:106
VX1
#define VX1
Definition: mod_defs.h:28
g_maxMach
double g_maxMach
The maximum Mach number computed during integration.
Definition: globals.h:119
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
KFASTP
Definition: mod_defs.h:58
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
KFASTM
Definition: mod_defs.h:58
VXn
int VXn
Definition: globals.h:73
BX3
#define BX3
Definition: mod_defs.h:27
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call 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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Here is the caller graph for this function:

void HLL_DIVB_SOURCE ( const State_1D *  ,
double **  ,
int  ,
int  ,
Grid *   
)

Definition at line 53 of file source.c.

68 {
69  int i, nv;
70  double vc[NVAR], *A, *src, *vm;
71  double r, s, vB;
72  static double *divB, *vp;
73  Grid *GG;
74 
75  if (divB == NULL){
76  divB = ARRAY_1D(NMAX_POINT, double);
77  vp = ARRAY_1D(NMAX_POINT, double);
78  }
79 
80  GG = grid + g_dir;
81  vm = vp - 1;
82  A = grid[g_dir].A;
83 
84 /* --------------------------------------------
85  Compute normal component of the field
86  -------------------------------------------- */
87 
88  for (i = beg - 1; i <= end; i++) {
89  vp[i] = state->flux[i][RHO] < 0.0 ? state->vR[i][BXn]: state->vL[i][BXn];
90  }
91 
92 /* --------------------------------------------
93  Compute div.B contribution from the normal
94  direction (1) in different geometries
95  -------------------------------------------- */
96 
97 
98  #if GEOMETRY == CARTESIAN
99 
100  for (i = beg; i <= end; i++) {
101  divB[i] = (vp[i] - vm[i])/GG->dx[i];
102  }
103 
104  #elif GEOMETRY == CYLINDRICAL
105 
106  if (g_dir == IDIR){ /* -- r -- */
107  for (i = beg; i <= end; i++) {
108  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/GG->dV[i];
109  }
110  }else if (g_dir == JDIR){ /* -- z -- */
111  for (i = beg; i <= end; i++) {
112  divB[i] = (vp[i] - vm[i])/GG->dx[i];
113  }
114  }
115 
116  #elif GEOMETRY == POLAR
117 
118  if (g_dir == IDIR){ /* -- r -- */
119  Ar = grid[IDIR].A;
120  for (i = beg; i <= end; i++) {
121  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/GG->dV[i];
122  }
123  }else if (g_dir == JDIR){ /* -- phi -- */
124  r = grid[IDIR].x[*g_i];
125  for (i = beg; i <= end; i++) {
126  divB[i] = (vp[i] - vm[i])/(r*GG->dx[i]);
127  }
128  }else if (g_dir == KDIR){ /* -- z -- */
129  for (i = beg; i <= end; i++) {
130  divB[i] = (vp[i] - vm[i])/GG->dx[i];
131  }
132  }
133 
134  #elif GEOMETRY == SPHERICAL
135 
136  if (g_dir == IDIR){ /* -- r -- */
137  for (i = beg; i <= end; i++) {
138  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/GG->dV[i];
139  }
140  }else if (g_dir == JDIR){ /* -- theta -- */
141  r = grid[IDIR].x[*g_i];
142  for (i = beg; i <= end; i++) {
143  divB[i] = (vp[i]*A[i] - vm[i]*A[i - 1])/(r*GG->dV[i]);
144  }
145  }else if (g_dir == KDIR){ /* -- phi -- */
146  r = grid[IDIR].x[*g_i];
147  s = sin(grid[JDIR].x[*g_j]);
148  for (i = beg; i <= end; i++) {
149  divB[i] = (vp[i] - vm[i])/(r*s*GG->dx[i]);
150  }
151  }
152 
153  #endif
154 
155  /* -----------------------------------------
156  compute total source terms
157  ----------------------------------------- */
158 
159  for (i = beg; i <= end; i++) {
160 
161  src = state->src[i];
162  for (nv = NFLX; nv--; ) vc[nv] = state->vh[i][nv];
163 
164  vB = EXPAND(vc[VX1]*vc[BX1], + vc[VX2]*vc[BX2], + vc[VX3]*vc[BX3]);
165  src[RHO] = 0.0;
166  EXPAND(src[MX1] = 0.0; ,
167  src[MX2] = 0.0; ,
168  src[MX3] = 0.0;)
169 
170  EXPAND(src[BX1] = -vc[VX1]*divB[i]; ,
171  src[BX2] = -vc[VX2]*divB[i]; ,
172  src[BX3] = -vc[VX3]*divB[i];)
173 
174  #if EOS != ISOTHERMAL
175  src[ENG] = 0.0;
176  #endif
177 
178  }
179 }
MX3
#define MX3
Definition: mod_defs.h:22
EOS
#define EOS
Definition: pluto.h:341
MX1
#define MX1
Definition: mod_defs.h:20
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
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
BXn
int BXn
Definition: globals.h:75
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
NFLX
#define NFLX
Definition: mod_defs.h:32
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
GRID
Definition: structs.h:78
MX2
#define MX2
Definition: mod_defs.h:21
divB
double *** divB
Definition: analysis.c:4
ISOTHERMAL
#define ISOTHERMAL
Definition: pluto.h:49
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
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
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609

Here is the caller graph for this function:

void HLL_Speed ( double **  vL,
double **  vR,
double *  a2L,
double *  a2R,
double *  hL,
double *  hR,
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

Definition at line 18 of file hll_speed.c.

35 {
36  int i, err;
37  static real *sl_min, *sl_max;
38  static real *sr_min, *sr_max;
39 
40  if (sl_min == NULL){
41  sl_min = ARRAY_1D(NMAX_POINT, double);
42  sl_max = ARRAY_1D(NMAX_POINT, double);
43 
44  sr_min = ARRAY_1D(NMAX_POINT, double);
45  sr_max = ARRAY_1D(NMAX_POINT, double);
46  }
47 
48 /* ----------------------------------------------
49  DAVIS Estimate
50  ---------------------------------------------- */
51 
52  MaxSignalSpeed (vL, a2L, hL, sl_min, sl_max, beg, end);
53  MaxSignalSpeed (vR, a2R, hR, sr_min, sr_max, beg, end);
54 /*
55  err = MAX_CH_SPEED (vL, sl_min, sl_max, beg, end);
56  if (err != 0) return err;
57  err = MAX_CH_SPEED (vR, sr_min, sr_max, beg, end);
58  if (err != 0) return err;
59 */
60  for (i = beg; i <= end; i++) {
61  SL[i] = MIN(sl_min[i], sr_min[i]);
62  SR[i] = MAX(sl_max[i], sr_max[i]);
63  }
64 /* return 0; */
65 }
MAX
#define MAX(a, b)
Definition: macros.h:101
real
double real
Definition: pluto.h:488
MIN
#define MIN(a, b)
Definition: macros.h:104
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
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:

int Magnetosonic ( double *  vp,
double  cs2,
double  h,
double *  lambda 
)

Find the four magneto-sonic speeds (fast and slow) for the relativistic MHD equations (RMHD). We follow the notations introduced in Eq. (16) in

Del Zanna, Bucciantini and Londrillo, A&A, 400, 397–413, 2003

and also Eq. (57-58) of Mignone & Bodo, MNRAS, 2006 for the degenerate cases.

The quartic equation is solved analytically and the solver (function 'QuarticSolve') assumes all roots are real (which should be always the case here, since eigenvalues are recovered from primitive variables). The coefficients of the quartic were found through the following MAPLE script: 

 *   ------------------------------------------------
restart;

u[0] := gamma;
u[x] := gamma*v[x];
b[0] := gamma*vB;
b[x] := B[x]/gamma + b[0]*v[x];
wt   := w + b^2;

#############################################################
"fdZ will be equation (16) of Del Zanna 2003 times w_{tot}";

e2  := c[s]^2 + b^2*(1 - c[s]^2)/wt;
fdZ := (1-e2)*(u[0]*lambda - u[x])^4 + (1-lambda^2)*
       (c[s]^2*(b[0]*lambda - b[x])^2/wt - e2*(u[0]*lambda - u[x])^2);

########################################################
"print the coefficients of the quartic polynomial fdZ";

coeff(fMB,lambda,4);
coeff(fMB,lambda,3);
coeff(fMB,lambda,2);
coeff(fMB,lambda,1);
coeff(fMB,lambda,0);

fdZ := fdZ*wt;  

######################################################
"fMB will be equation (56) of Mignone & Bodo (2006)";

a  := gamma*(lambda-v[x]);
BB := b[x] - lambda*b[0];
fMB := w*(1-c[s]^2)*a^4 - (1-lambda^2)*((b^2+w*c[s]^2)*a^2-c[s]^2*BB^2);

######################################
"check that fdZ = fMB";

simplify(fdZ - fMB);

########################################################
"print the coefficients of the quartic polynomial fMB";

coeff(fMB,lambda,4);
coeff(fMB,lambda,3);
coeff(fMB,lambda,2);
coeff(fMB,lambda,1);
coeff(fMB,lambda,0);


###############################################
"Degenerate case 2: Bx = 0, Eq (58) in MB06";

B[x] := 0;

a2 := w*(c[s]^2 + gamma^2*(1-c[s]^2)) + Q; 
a1 := -2*w*gamma^2*v[x]*(1-c[s]^2);
a0 := w*(-c[s]^2 + gamma^2*v[x]^2*(1-c[s]^2))-Q;

"delc is a good way to evaluate the determinant so that it is > 0";
del  := a1^2 - 4*a2*a0;
delc := 4*(w*c[s]^2 + Q)*(w*(1-c[s]^2)*gamma^2*(1-v[x]^2) + (w*c[s]^2+Q));
simplify(del-delc);

############################################################################
"check the correctness of the coefficients of Eq. (58) in MB06";

Q  := b^2 - c[s]^2*vB^2;
divide(fMB, (lambda-v[x])^2,'fMB2');
simplify(a2 - coeff(fMB2,lambda,2)/gamma^2);
simplify(a1 - coeff(fMB2,lambda,1)/gamma^2);
simplify(a0 - coeff(fMB2,lambda,0)/gamma^2);

Definition at line 34 of file eigenv.c.

127 {
128  int iflag;
129  double vB, vB2, Bm2, vm2, w, w_1;
130  double eps2, one_m_eps2, a2_w;
131  double vx, vx2, u0, u02;
132  double a4, a3, a2, a1, a0;
133  double scrh1, scrh2, scrh;
134  double b2, del, z[4];
135 
136 #if RMHD_FAST_EIGENVALUES
137  Eigenvalues (vp, cs2, h, lambda);
138  return 0;
139 #endif
140 
141  scrh = fabs(vp[VXn])/sqrt(cs2);
142  g_maxMach = MAX(scrh, g_maxMach);
143 
144  vB = EXPAND(vp[VX1]*vp[BX1], + vp[VX2]*vp[BX2], + vp[VX3]*vp[BX3]);
145  u02 = EXPAND(vp[VX1]*vp[VX1], + vp[VX2]*vp[VX2], + vp[VX3]*vp[VX3]);
146  Bm2 = EXPAND(vp[BX1]*vp[BX1], + vp[BX2]*vp[BX2], + vp[BX3]*vp[BX3]);
147  vm2 = u02;
148 
149  if (u02 >= 1.0){
150  print ("! MAGNETOSONIC: |v|= %f > 1\n",u02);
151  return 1;
152  }
153 
154  if (u02 < 1.e-14) {
155 
156  /* -----------------------------------------------------
157  if total velocity is = 0, the eigenvalue equation
158  reduces to a biquadratic:
159 
160  x^4 + a1*x^2 + a0 = 0
161 
162  with
163  a0 = cs^2*bx*bx/W > 0,
164  a1 = - a0 - eps^2 < 0.
165  ----------------------------------------------------- */
166 
167  w_1 = 1.0/(vp[RHO]*h + Bm2);
168  eps2 = cs2 + Bm2*w_1*(1.0 - cs2);
169  a0 = cs2*vp[BXn]*vp[BXn]*w_1;
170  a1 = - a0 - eps2;
171  del = a1*a1 - 4.0*a0;
172  del = MAX(del, 0.0);
173  del = sqrt(del);
174  lambda[KFASTP] = sqrt(0.5*(-a1 + del));
175  lambda[KFASTM] = -lambda[KFASTP];
176  #if COMPONENTS > 1
177  lambda[KSLOWP] = sqrt(0.5*(-a1 - del));
178  lambda[KSLOWM] = -lambda[KSLOWP];
179  #endif
180  return 0;
181  }
182 
183  vB2 = vB*vB;
184  u02 = 1.0/(1.0 - u02);
185  b2 = Bm2/u02 + vB2;
186  u0 = sqrt(u02);
187  w_1 = 1.0/(vp[RHO]*h + b2);
188  vx = vp[VXn];
189  vx2 = vx*vx;
190 
191  if (fabs(vp[BXn]) < 1.e-14){
192 
193  w = vp[RHO]*h;
194  scrh1 = b2 + cs2*(w - vB2); /* wc_s^2 + Q */
195  scrh2 = w*(1.0 - cs2)*u02;
196  del = 4.0*scrh1*(scrh2*(1.0 - vx2) + scrh1);
197 
198  a2 = scrh1 + scrh2;
199  a1 = -2.0*vx*scrh2;
200 
201  lambda[KFASTP] = 0.5*(-a1 + sqrt(del))/a2;
202  lambda[KFASTM] = 0.5*(-a1 - sqrt(del))/a2;
203  #if COMPONENTS > 1
204  lambda[KSLOWP] = lambda[KSLOWM] = vp[VXn];
205  #endif
206 
207  return 0;
208 
209  }else{
210 
211  scrh1 = vp[BXn]/u02 + vB*vx; /* -- this is bx/u0 -- */
212  scrh2 = scrh1*scrh1;
213 
214  a2_w = cs2*w_1;
215  eps2 = (cs2*vp[RHO]*h + b2)*w_1;
216  one_m_eps2 = u02*vp[RHO]*h*(1.0 - cs2)*w_1;
217 
218  /* ---------------------------------------
219  Define coefficients for the quartic
220  --------------------------------------- */
221 
222  scrh = 2.0*(a2_w*vB*scrh1 - eps2*vx);
223  a4 = one_m_eps2 - a2_w*vB2 + eps2;
224  a3 = - 4.0*vx*one_m_eps2 + scrh;
225  a2 = 6.0*vx2*one_m_eps2 + a2_w*(vB2 - scrh2) + eps2*(vx2 - 1.0);
226  a1 = - 4.0*vx*vx2*one_m_eps2 - scrh;
227  a0 = vx2*vx2*one_m_eps2 + a2_w*scrh2 - eps2*vx2;
228 
229  if (a4 < 1.e-12){
230  print ("! Magnetosonic: cannot divide by a4\n");
231  return 1;
232  }
233 
234  scrh = 1.0/a4;
235 
236  a3 *= scrh;
237  a2 *= scrh;
238  a1 *= scrh;
239  a0 *= scrh;
240  iflag = QuarticSolve(a3, a2, a1, a0, z);
241 
242  if (iflag){
243  print ("! Magnetosonic: Cannot find max speed [dir = %d]\n", g_dir);
244  print ("! rho = %12.6e\n",vp[RHO]);
245  print ("! prs = %12.6e\n",vp[PRS]);
246  EXPAND(print ("! vx = %12.6e\n",vp[VX1]); ,
247  print ("! vy = %12.6e\n",vp[VX2]); ,
248  print ("! vz = %12.6e\n",vp[VX3]);)
249  EXPAND(print ("! Bx = %12.6e\n",vp[BX1]); ,
250  print ("! By = %12.6e\n",vp[BX2]); ,
251  print ("! Bz = %12.6e\n",vp[BX3]);)
252  print ("! f(x) = %12.6e + x*(%12.6e + x*(%12.6e ",
253  a0*a4, a1*a4, a2*a4);
254  print (" + x*(%12.6e + x*%12.6e)))\n", a3*a4, a4);
255  return 1;
256  }
257 /*
258 if (z[3] < z[2] || z[3] < z[1] || z[3] < z[0] ||
259  z[2] < z[1] || z[2] < z[0] || z[1] < z[0]){
260  printf ("!Sorting not correct\n");
261  printf ("z = %f %f %f %f\n",z[0],z[1],z[2],z[3]);
262  QUIT_PLUTO(1);
263 }
264 */
265 
266  lambda[KFASTM] = z[0];
267  lambda[KFASTP] = z[3];
268  #if COMPONENTS > 1
269  lambda[KSLOWP] = z[2];
270  lambda[KSLOWM] = z[1];
271  #endif
272  }
273  return 0;
274 }
MAX
#define MAX(a, b)
Definition: macros.h:101
QuarticSolve
int QuarticSolve(double, double, double, double, double *)
Definition: quartic.c:23
VX2
#define VX2
Definition: mod_defs.h:29
Eigenvalues
void Eigenvalues(double **v, double *csound2, double **lambda, int beg, int end)
Definition: eigenv.c:67
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
VX1
#define VX1
Definition: mod_defs.h:28
BXn
int BXn
Definition: globals.h:75
g_maxMach
double g_maxMach
The maximum Mach number computed during integration.
Definition: globals.h:119
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
KFASTP
Definition: mod_defs.h:58
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
KFASTM
Definition: mod_defs.h:58
VXn
int VXn
Definition: globals.h:73
BX3
#define BX3
Definition: mod_defs.h:27
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26

Here is the call graph for this function:

Here is the caller graph for this function:

int MaxSignalSpeed ( double **  v,
double *  a2,
double *  h,
double *  cmin,
double *  cmax,
int  beg,
int  end 
)

Return the rightmost (cmax) and leftmost (cmin) wave propagation speed in the Riemann fan

Definition at line 13 of file eigenv.c.

20 {
21  int i, err;
22  double lambda[NFLX];
23 
24  for (i = beg; i <= end; i++) {
25  err = Magnetosonic (v[i], a2[i], h[i], lambda);
26  if (err != 0) return i;
27  cmax[i] = lambda[KFASTP];
28  cmin[i] = lambda[KFASTM];
29  }
30  return 0;
31 }
func_param::v
double v[NVAR]
Definition: eos.h:106
NFLX
#define NFLX
Definition: mod_defs.h:32
KFASTP
Definition: mod_defs.h:58
KFASTM
Definition: mod_defs.h:58
i
int i
Definition: analysis.c:2
Magnetosonic
int Magnetosonic(double *vp, double cs2, double h, double *lambda)
Definition: eigenv.c:34

Here is the call graph for this function:

void POWELL_DIVB_SOURCE ( const State_1D *  ,
int  ,
int  ,
Grid *   
)

Definition at line 5 of file source.c.

10 {
11  int i;
12  double divB;
13  double g_2, vB, *vc;
14  Grid *GG;
15 
16  GG = grid + g_dir;
17 
18 /* ------------------------------------------------------------
19  POWELL's divB monopole source term
20  ------------------------------------------------------------ */
21 
22  for (i = is; i <= ie; i++) {
23 
24  vc = state->vh[i];
25 
26 /* Make (divB) contribution from the normal direction (1) */
27 
28  divB = (state->vR[i][BXn] + state->vL[i][BXn])*GG->A[i] -
29  (state->vR[i - 1][BXn] + state->vL[i - 1][BXn])*GG->A[i - 1];
30 /*
31  divB = state->bn[i]*GG->A[i] - state->bn[i - 1]*GG->A[i - 1];
32 */
33  divB /= 2.0*GG->dV[i];
34 
35  vB = EXPAND(vc[VX1]*vc[BX1], + vc[VX2]*vc[BX2], + vc[VX3]*vc[BX3]);
36  g_2 = EXPAND(vc[VX1]*vc[VX1], + vc[VX2]*vc[VX2], + vc[VX3]*vc[VX3]);
37  g_2 = 1.0 - g_2;
38 
39  EXPAND (state->src[i][MX1] = -divB*(vc[BX1]*g_2 + vB*vc[VX1]); ,
40  state->src[i][MX2] = -divB*(vc[BX2]*g_2 + vB*vc[VX2]); ,
41  state->src[i][MX3] = -divB*(vc[BX3]*g_2 + vB*vc[VX3]);)
42 
43  state->src[i][ENG] = -divB*vB;
44 
45  EXPAND (state->src[i][BX1] = -divB*vc[VX1]; ,
46  state->src[i][BX2] = -divB*vc[VX2]; ,
47  state->src[i][BX3] = -divB*vc[VX3];)
48 
49  }
50 }
MX3
#define MX3
Definition: mod_defs.h:22
MX1
#define MX1
Definition: mod_defs.h:20
VX2
#define VX2
Definition: mod_defs.h:29
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
VX1
#define VX1
Definition: mod_defs.h:28
BXn
int BXn
Definition: globals.h:75
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
GRID
Definition: structs.h:78
MX2
#define MX2
Definition: mod_defs.h:21
divB
double *** divB
Definition: analysis.c:4
BX3
#define BX3
Definition: mod_defs.h:27
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

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 call graph for this function:

void PRIM_EIGENVECTORS ( double *  ,
double  ,
double  ,
double *  ,
double **  ,
double **   
)
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

Here is the call graph for this function:

Here is the caller graph for this function:

int QuarticSolve ( double  b,
double  c,
double  d,
double  e,
double *  z 
)

Solve a quartic equation in the form

\[ z^4 + bz^3 + cz^2 + dz + e = 0 \]

For its purpose, it is assumed that ALL roots are double. This makes things faster.

Parameters
[in]bcoefficient of the quartic
[in]ccoefficient of the quartic
[in]dcoefficient of the quartic
[in]ecoefficient of the quartic
[out]zvector containing the (double) roots of the quartic

Reference:

http://www.1728.com/quartic2.htm

Returns
Return 0 on success, 1 if cubic solver fails, 2 if NaN has been found and 3 if quartic is not satisfied.

Definition at line 23 of file quartic.c.

47 {
48  int status;
49  status = QuarticSolveNew(b,c,d,e,z);
50 
51  if (status != 0){
52  debug_print = 1;
53  QuarticSolveNew(b,c,d,e,z);
54  QUIT_PLUTO(1);
55  }
56 }
QuarticSolveNew
int QuarticSolveNew(double b, double c, double d, double e, double *z)
Definition: quartic.c:175
menu.c
tuple c
Definition: menu.py:375
debug_print
static int debug_print
Definition: quartic.c:20
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 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

Variable Documentation

Riemann_Solver HLL_Linde_Solver

Definition at line 221 of file mod_defs.h.

Riemann_Solver HLL_Solver

Definition at line 221 of file mod_defs.h.

Riemann_Solver HLLC_Solver

Definition at line 221 of file mod_defs.h.

void HLLD_Solver

Solve Riemann problem for the isothermal MHD equations using the three-state HLLD Riemann solver of Mignone (2007).

Parameters
[in,out]statepointer to State_1D structure
[in]beginitial grid index
[out]endfinal grid index
[out]cmax1D array of maximum characteristic speeds
[in]gridpointer to array of Grid structures.

Definition at line 221 of file mod_defs.h.

Riemann_Solver LF_Solver

Definition at line 221 of file mod_defs.h.