PLUTO
Functions
prim_eqn.c File Reference

Compute the right hand side of the primitive MHD equations. More...

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

Go to the source code of this file.

Functions

void PrimRHS (double *v, double *dv, double cs2, double h, double *Adv)
 
void PrimSource (const State_1D *state, int beg, int end, double *a2, double *h, double **src, Grid *grid)
 

Detailed Description

Compute the right hand side of the primitive MHD equations.

Implements the right hand side of the quasi-linear form of the MHD equations. In 1D this may be written as

\[ \partial_t{\mathbf{V}} = - A\cdot\partial_x\mathbf{V} + \mathbf{S} \]

where $ A $ is the matrix of the primitive form of the equations, $ S $ is the source term.

Reference

  • "A solution adaptive upwind scheme for ideal MHD", Powell et al., JCP (1999) 154, 284.

The function PrimRHS() implements the first term while PrimSource() implements the source term part.

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

Definition in file prim_eqn.c.

Function Documentation

void PrimRHS ( double *  v,
double *  dv,
double  cs2,
double  h,
double *  Adv 
)

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.
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 30 of file prim_eqn.c.

52 {
53  int nv;
54  double tau, scrh;
55  double ch2;
56 
57  tau = 1.0/v[RHO];
58 
59 /* ---------------------------------------------
60  Adv[k] Contains A[k][*]*dv[*]
61  --------------------------------------------- */
62 
63  Adv[RHO] = v[VXn]*dv[RHO] + v[RHO]*dv[VXn];
64  scrh = EXPAND(0.0, + v[BXt]*dv[BXt], + v[BXb]*dv[BXb]);
65 
66  #if EOS == IDEAL || EOS == PVTE_LAW
67  Adv[VXn] = v[VXn]*dv[VXn] + tau*(dv[PRS] + scrh);
68  #elif EOS == ISOTHERMAL
69  Adv[VXn] = v[VXn]*dv[VXn] + tau*(cs2*dv[RHO] + scrh);
70  #else
71  print ("! PRIM_RHS: not defined for this EoS\n");
72  QUIT_PLUTO(1);
73  #endif
74 
75  EXPAND( ; ,
76  Adv[VXt] = v[VXn]*dv[VXt] - tau*v[BXn]*dv[BXt]; ,
77  Adv[VXb] = v[VXn]*dv[VXb] - tau*v[BXn]*dv[BXb]; )
78 
79  #if DIVB_CONTROL == EIGHT_WAVES
80  Adv[BXn] = v[VXn]*dv[BXn];
81  #elif DIVB_CONTROL == DIV_CLEANING
82  ch2 = glm_ch*glm_ch;
83  Adv[BXn] = dv[PSI_GLM];
84  Adv[PSI_GLM] = dv[BXn]*ch2;
85  #else
86  Adv[BXn] = 0.0;
87  #endif
88 
89 /* ------------------------------------------------------------------- */
90 /*! \note
91  In the 7-wave and 8-wave formulations we use the the same matrix
92  being decomposed into right and left eigenvectors during
93  the Characteristic Tracing step.
94  Note, however, that it DOES NOT include two additional terms
95  (-vy*dV[BX] for By, -vz*dv[BX] for Bz) that are needed in the
96  7-wave form and are added using source terms.
97  --------------------------------------------------------------------- */
98 
99  EXPAND( ; ,
100  Adv[BXt] = v[BXt]*dv[VXn] - v[BXn]*dv[VXt] + v[VXn]*dv[BXt]; ,
101  Adv[BXb] = v[BXb]*dv[VXn] - v[BXn]*dv[VXb] + v[VXn]*dv[BXb];)
102 
103  #if EOS == IDEAL
104  Adv[PRS] = g_gamma*v[PRS]*dv[VXn] + v[VXn]*dv[PRS];
105  #elif EOS == PVTE_LAW
106  Adv[PRS] = cs2*v[RHO]*dv[VXn] + v[VXn]*dv[PRS];
107  #endif
108 
109 /* -------------------------------------------------------------
110  Now Define Tracers
111  ------------------------------------------------------------- */
112 
113 #if NSCL > 0
114  NSCL_LOOP(nv) Adv[nv] = v[VXn]*dv[nv];
115 #endif
116 
117 }
IDEAL
#define IDEAL
Definition: pluto.h:45
EOS
#define EOS
Definition: pluto.h:341
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
g_gamma
double g_gamma
Definition: globals.h:112
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
func_param::v
double v[NVAR]
Definition: eos.h:106
DIVB_CONTROL
#define DIVB_CONTROL
Definition: definitions_01.h:17
EIGHT_WAVES
#define EIGHT_WAVES
Definition: pluto.h:120
VXb
int VXb
Definition: globals.h:73
BXn
int BXn
Definition: globals.h:75
VXt
int VXt
Definition: globals.h:73
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
VXn
int VXn
Definition: globals.h:73
glm_ch
double glm_ch
The propagation speed of divergence error.
Definition: glm.c:21
BXt
int BXt
Definition: globals.h:75
BXb
int BXb
Definition: globals.h:75
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

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 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) \]

Definition at line 120 of file prim_eqn.c.

162 {
163  int nv, i, j, k;
164  double tau, dA_dV, th;
165  double hscale; /* scale factor */
166  double *v, *vp, *A, *dV, r_inv, ct;
167  double *x1, *x2, *x3;
168  double *x1p, *x2p, *x3p;
169  double *dx1, *dx2, *dx3;
170  static double *phi_p;
171  double g[3], ch2, db, scrh;
172 
173  #if ROTATING_FRAME == YES
174  print1 ("! PrimSource: does not work with rotations\n");
175  QUIT_PLUTO(1);
176  #endif
177 
178 /* ----------------------------------------------------------
179  1. Memory allocation and pointer shortcuts
180  ---------------------------------------------------------- */
181 
182  if (phi_p == NULL) phi_p = ARRAY_1D(NMAX_POINT, double);
183 
184  #if GEOMETRY == CYLINDRICAL
185  x1 = grid[IDIR].xgc; x1p = grid[IDIR].xr; dx1 = grid[IDIR].dx;
186  x2 = grid[JDIR].xgc; x2p = grid[JDIR].xr; dx2 = grid[JDIR].dx;
187  x3 = grid[KDIR].xgc; x3p = grid[KDIR].xr; dx3 = grid[KDIR].dx;
188  #else
189  x1 = grid[IDIR].x; x1p = grid[IDIR].xr; dx1 = grid[IDIR].dx;
190  x2 = grid[JDIR].x; x2p = grid[JDIR].xr; dx2 = grid[JDIR].dx;
191  x3 = grid[KDIR].x; x3p = grid[KDIR].xr; dx3 = grid[KDIR].dx;
192  #endif
193 
194  #ifdef GLM_MHD
195  ch2 = glm_ch*glm_ch;
196  #endif
197 
198  A = grid[g_dir].A;
199  dV = grid[g_dir].dV;
200  hscale = 1.0;
201 
202  i = g_i; j = g_j; k = g_k;
203 
204 /* ----------------------------------------------------------
205  initialize all elements of src to zero
206  ---------------------------------------------------------- */
207 
208  memset((void *)src[0], '\0',NMAX_POINT*NVAR*sizeof(double));
209 
210 /* ----------------------------------------------------------
211  2. Compute geometrical source terms
212  ---------------------------------------------------------- */
213 
214 #if GEOMETRY == CYLINDRICAL
215 
216  if (g_dir == IDIR) {
217  for (i = beg; i <= end; i++){
218  v = state->v[i];
219 
220  tau = 1.0/v[RHO];
221  dA_dV = 1.0/x1[i];
222 
223  src[i][RHO] = -v[RHO]*v[VXn]*dA_dV;
224  EXPAND( ,
225  src[i][iBZ] = (v[iBR]*v[iVZ] - v[iVR]*v[iBZ])*dA_dV; ,
226  src[i][iVR] = (v[iVPHI]*v[iVPHI] - v[iBPHI]*v[iBPHI]*tau)*dA_dV;
227  src[i][iVPHI] = (-v[iVR]*v[iVPHI] + v[iBR]*v[iBPHI]*tau)*dA_dV;)
228 
229  #if EOS == IDEAL
230  src[i][PRS] = a2[i]*src[i][RHO];
231  #endif
232 
233  #ifdef GLM_MHD
234  src[i][PSI_GLM] = -v[iBR]*dA_dV*ch2;
235  #endif
236  }
237  }
238 
239 #elif GEOMETRY == POLAR
240 
241  if (g_dir == IDIR) {
242  for (i = beg; i <= end; i++){
243  v = state->v[i];
244 
245  tau = 1.0/v[RHO];
246  dA_dV = 1.0/x1[i];
247  src[i][RHO] = -v[RHO]*v[VXn]*dA_dV;
248 
249  EXPAND( ,
250  src[i][iVR] = (v[iVPHI]*v[iVPHI] - v[iBPHI]*v[iBPHI]*tau)*dA_dV;
251  src[i][iVPHI] = (-v[iVR]*v[iVPHI] + v[iBR]*v[iBPHI]*tau)*dA_dV; ,
252  src[i][iBZ] = ( v[iBR]*v[iVZ] - v[iVR]*v[iBZ])*dA_dV;)
253 
254  #if EOS == IDEAL
255  src[i][PRS] = a2[i]*src[i][RHO];
256  #endif
257  #ifdef GLM_MHD
258  src[i][PSI_GLM] = -v[iBR]*dA_dV*ch2;
259  #endif
260 
261  }
262  }
263 
264 #elif GEOMETRY == SPHERICAL
265 
266  print1 ("! PrimSource: not implemented in Spherical geometry\n");
267  QUIT_PLUTO(1);
268 
269 #endif
270 
271 /* ----------------------------------------------------------
272  3. Add body forces. This includes:
273  - Coriolis terms for the shearing box module
274  - Body forces
275  ---------------------------------------------------------- */
276 
277 #ifdef SHEARINGBOX
278 
279  if (g_dir == IDIR){
280  for (i = beg; i <= end; i++) {
281  src[i][VX1] = 2.0*state->v[i][VX2]*SB_OMEGA;
282  }
283  }else if (g_dir == JDIR){
284  for (j = beg; j <= end; j++) {
285  #ifdef FARGO
286  src[j][VX2] = (SB_Q - 2.0)*state->v[j][VX1]*SB_OMEGA;
287  #else
288  src[j][VX2] = -2.0*state->v[j][VX1]*SB_OMEGA;
289  #endif
290  }
291  }
292 
293 #endif
294 
295 
296 #if (BODY_FORCE != NO)
297  if (g_dir == IDIR) {
298 
299  i = beg-1;
300  j = g_j;
301  k = g_k;
302  #if BODY_FORCE & POTENTIAL
303  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
304  #endif
305  for (i = beg; i <= end; i++){
306  #if BODY_FORCE & VECTOR
307  v = state->v[i];
308  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
309  src[i][VX1] += g[IDIR];
310  #endif
311  #if BODY_FORCE & POTENTIAL
312  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
313  src[i][VX1] -= (phi_p[i] - phi_p[i-1])/(hscale*dx1[i]);
314  #endif
315 
316  /* -- Add tangential components in 1D -- */
317 
318  #if DIMENSIONS == 1
319  EXPAND( ,
320  src[i][VX2] += g[JDIR]; ,
321  src[i][VX3] += g[KDIR];)
322  #endif
323  }
324 
325  }else if (g_dir == JDIR){
326 
327  i = g_i;
328  j = beg - 1;
329  k = g_k;
330  #if BODY_FORCE & POTENTIAL
331  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
332  #endif
333  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
334  hscale = x1[i];
335  #endif
336  for (j = beg; j <= end; j++){
337  #if BODY_FORCE & VECTOR
338  v = state->v[j];
339  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
340  src[j][VX2] += g[JDIR];
341  #endif
342  #if BODY_FORCE & POTENTIAL
343  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
344  src[j][VX2] -= (phi_p[j] - phi_p[j-1])/(hscale*dx2[j]);
345  #endif
346 
347  /* -- Add 3rd component in 2D -- */
348 
349  #if DIMENSIONS == 2 && COMPONENTS == 3
350  src[j][VX3] += g[KDIR];
351  #endif
352  }
353 
354  }else if (g_dir == KDIR){
355 
356  i = g_i;
357  j = g_j;
358  k = beg - 1;
359  #if BODY_FORCE & POTENTIAL
360  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
361  #endif
362  #if GEOMETRY == SPHERICAL
363  th = x2[j];
364  hscale = x1[i]*sin(th);
365  #endif
366  for (k = beg; k <= end; k++){
367  #if BODY_FORCE & VECTOR
368  v = state->v[k];
369  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
370  src[k][VX3] += g[KDIR];
371  #endif
372  #if BODY_FORCE & POTENTIAL
373  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
374  src[k][VX3] -= (phi_p[k] - phi_p[k-1])/(hscale*dx3[k]);
375  #endif
376  }
377  }
378 #endif
379 
380 /* -----------------------------------------------------------
381  4. MHD, div.B related source terms
382  ----------------------------------------------------------- */
383 
384  #ifdef GLM_MHD
385  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
386  print1 ("! Error: Div. Cleaning does not work in this configuration.\n");
387  print1 ("! Try RK integrator instead\n");
388  QUIT_PLUTO(1);
389  #endif
390  for (i = beg; i <= end; i++){
391  v = state->v[i];
392 
393  tau = 1.0/v[RHO];
394  db = 0.5*( A[i] *(state->v[i+1][BXn] + state->v[i][BXn])
395  - A[i-1]*(state->v[i-1][BXn] + state->v[i][BXn]))/dV[i];
396  #if GLM_EXTENDED == NO
397  EXPAND(src[i][VXn] += v[BXn]*tau*db; ,
398  src[i][VXt] += v[BXt]*tau*db; ,
399  src[i][VXb] += v[BXb]*tau*db;)
400  #endif
401  EXPAND( ,
402  src[i][BXt] += v[VXt]*db; ,
403  src[i][BXb] += v[VXb]*db;)
404 
405  #if EOS == IDEAL
406  scrh = EXPAND(v[VXn]*v[BXn], + v[VXt]*v[BXt], + v[VXb]*v[BXb]);
407  src[i][PRS] += (1.0 - g_gamma)*scrh*db;
408  #if GLM_EXTENDED == NO
409  scrh = 0.5*(state->v[i+1][PSI_GLM] - state->v[i-1][PSI_GLM])/grid[g_dir].dx[i];
410  src[i][PRS] += (g_gamma - 1.0)*v[BXn]*scrh;
411  #endif
412  #endif
413  }
414  #endif
415 
416 /* ---------------------------------------------------------------
417  5. Add FARGO source terms (except for SHEARINGBOX).
418  When SHEARINGBOX is also enabled, we do not include
419  these source terms since they're all provided by body_force)
420  --------------------------------------------------------------- */
421 
422  #if (defined FARGO && !defined SHEARINGBOX)
423  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
424  print1 ("! Time Stepping works only in Cartesian or cylindrical coords\n");
425  print1 ("! Use RK instead\n");
426  QUIT_PLUTO(1);
427  #endif
428 
429  double **wA, *dx, *dz;
430  wA = FARGO_GetVelocity();
431  if (g_dir == IDIR){
432  k = g_k;
433  dx = grid[IDIR].dx;
434  for (i = beg; i <= end; i++){
435  v = state->v[i];
436  src[i][VX2] -= 0.5*v[VX1]*(wA[k][i+1] - wA[k][i-1])/dx[i];
437  }
438  }else if (g_dir == KDIR){
439  i = g_i;
440  dz = grid[KDIR].dx;
441  for (k = beg; k <= end; k++){
442  v = state->v[k];
443  src[k][VX2] -= 0.5*v[VX3]*(wA[k+1][i] - wA[k-1][i])/dz[k];
444  }
445  }
446  #endif
447 }
IDEAL
#define IDEAL
Definition: pluto.h:45
EOS
#define EOS
Definition: pluto.h:341
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
g_gamma
double g_gamma
Definition: globals.h:112
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
iBZ
#define iBZ
Definition: mod_defs.h:138
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
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
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
iBR
#define iBR
Definition: mod_defs.h:150
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
VXb
int VXb
Definition: globals.h:73
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
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
iBPHI
#define iBPHI
Definition: mod_defs.h:152
VXt
int VXt
Definition: globals.h:73
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
glm_ch
double glm_ch
The propagation speed of divergence error.
Definition: glm.c:21
iVZ
#define iVZ
Definition: mod_defs.h:55
JDIR
#define JDIR
Definition: pluto.h:194
BXt
int BXt
Definition: globals.h:75
NVAR
#define NVAR
Definition: pluto.h:609
BXb
int BXb
Definition: globals.h:75
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: