PLUTO
Functions
eigenv.c File Reference

Compute the eigenvalues for the relativisitc MHD equations. More...

#include "pluto.h"
Include dependency graph for eigenv.c:

Go to the source code of this file.

Functions

int MaxSignalSpeed (double **v, double *a2, double *h, double *cmin, double *cmax, int beg, int end)
 
int Magnetosonic (double *vp, double cs2, double h, double *lambda)
 
int Eigenvalues (double *v, double cs2, double h, double *lambda)
 

Detailed Description

Compute the eigenvalues for the relativisitc MHD equations.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Sep 10, 2012

Definition in file eigenv.c.

Function Documentation

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 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: