PLUTO
rhs.c
Go to the documentation of this file.
1 #include "pluto.h"
2 
3 #ifdef CH_SPACEDIM /* implies Chombo is being used */
4  #define USE_PR_GRADIENT YES
5 #else
6  #ifdef FINITE_DIFFERENCE
7  #define USE_PR_GRADIENT NO /* -- default for Finite Difference schemes -- */
8  #else
9  #define USE_PR_GRADIENT YES /* -- default, do not change!! -- */
10  #endif
11 #endif
12 
13 /* *********************************************************** */
14 void RightHandSide (const State_1D *state, Time_Step *Dts,
15  int beg, int end, double dt, Grid *grid)
16 /*
17  * Compute right hand side of the MHD equations in different geometries,
18  * taking contributions from one direction at a time.
19  * The right hand side (rhs) consists of the following contributions:
20  *
21  * rhs = dt/dV * (Ap*Fp - Am*Fm) + dt*S
22  *
23  * where
24  *
25  * Ap, Am: interface areas,
26  * Fp, Fm: interface fluxes,
27  * dt: time step
28  * S: source term including
29  * * geometrical source terms
30  * * gravity
31  *
32  * In order to compute rhs, this function takes the following steps:
33  *
34  * - initialize rhs with flux differences (#1)
35  * - add geometrical source terms (#2)
36  * - add gravity (#4)
37  *
38  * LAST MODIFIED
39  *
40  * July, 31 2012 by A. Mignone (mignone@ph.unito.it)
41  *
42  ************************************************************* */
43 {
44  int i, j, k, nv;
45  double dtdx, dtdV, scrh;
46  double *x1, *x1p, *dx1, *dV1;
47  double *x2, *x2p, *dx2, *dV2;
48  double *x3, *x3p, *dx3, *dV3;
49  double **flux, **rhs, *p, *v;
50  double cl;
51  double g[3];
52  static double **fA, *h;
53 
54  #if GEOMETRY != CARTESIAN
55  if (fA == NULL) {
56  fA = ARRAY_2D(NMAX_POINT, NVAR, double);
57  h = ARRAY_1D(NMAX_POINT, double);
58  }
59  #endif
60 
61 /* --------------------------------------------------
62  Compute passive scalar fluxes
63  -------------------------------------------------- */
64 
65  #if NSCL > 0
66  AdvectFlux (state, beg - 1, end, grid);
67  #endif
68 
69 /* --------------------------
70  pointer shortcuts
71  -------------------------- */
72 
73  rhs = state->rhs;
74  flux = state->flux;
75  p = state->press;
76 
77  x1 = grid[IDIR].x; x2 = grid[JDIR].x; x3 = grid[KDIR].x;
78  x1p = grid[IDIR].xr; x2p = grid[JDIR].xr; x3p = grid[KDIR].xr;
79  dx1 = grid[IDIR].dx; dx2 = grid[JDIR].dx; dx3 = grid[KDIR].dx;
80  dV1 = grid[IDIR].dV; dV2 = grid[JDIR].dV; dV3 = grid[KDIR].dV;
81 
82  i = g_i; /* will be redefined during x1-sweep */
83  j = g_j; /* will be redefined during x2-sweep */
84  k = g_k; /* will be redefined during x3-sweep */
85 
86 /* ------------------------------------------------
87  Add pressure to normal component of
88  momentum flux if necessary.
89  ------------------------------------------------ */
90 
91  #if USE_PR_GRADIENT == NO
92  for (i = beg - 1; i <= end; i++) flux[i][MXn] += p[i];
93  #endif
94 
95 #if GEOMETRY == CARTESIAN
96 /* ***********************************************************
97 
98  Compute right-hand side of the RMHD/RHD
99  equations in CARTESIAN geometry.
100 
101  *********************************************************** */
102 {
103  double x, y, z;
104 
105  if (g_dir == IDIR){
106 
107  /* ****************************************************
108  Cartesian x-direction,
109 
110  - initialize rhs with flux differences (I1)
111  - enforce conservation of total angular
112  momentum and/or energy (I3)
113  - add gravity (I4)
114  **************************************************** */
115 
116  y = x2[j];
117  z = x3[k];
118  for (i = beg; i <= end; i++) {
119  x = x1[i];
120  dtdx = dt/dx1[i];
121 
122  /* -----------------------------------------------
123  I1. initialize rhs with flux difference
124  ----------------------------------------------- */
125 
126  NVAR_LOOP(nv) rhs[i][nv] = -dtdx*(flux[i][nv] - flux[i-1][nv]);
127  #if USE_PR_GRADIENT == YES
128  rhs[i][MX1] -= dtdx*(p[i] - p[i-1]);
129  #endif
130 
131  /* ----------------------------------------------------
132  I4. Include gravity
133  ---------------------------------------------------- */
134 
135  v = state->vh[i];
136  #if (BODY_FORCE & VECTOR)
137  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
138  rhs[i][MX1] += dt*v[RHO]*g[IDIR];
139  #if DIMENSIONS == 1 && COMPONENTS > 1 /* In 1D and when COMPONENTS == 2 */
140  /* or 3, add gravity contributions */
141  /* along non-existing dimensions. */
142  EXPAND( ,
143  rhs[i][MX2] += dt*vc[RHO]*g[JDIR]; ,
144  rhs[i][MX3] += dt*vc[RHO]*g[KDIR];)
145  #endif
146 
147  #if HAVE_ENERGY
148  rhs[i][ENG] += dt*0.5*(flux[i][RHO] + flux[i-1][RHO])*g[IDIR];
149  #if DIMENSIONS == 1 && COMPONENTS > 1
150  rhs[i][ENG] += dt*(EXPAND(0.0, + vc[RHO]*vc[VX2]*g[JDIR],
151  + vc[RHO]*vc[VX3]*g[KDIR]));
152  #endif
153  #endif
154  #endif
155  }
156  } else if (g_dir == JDIR){
157 
158  /* ****************************************************
159  Cartesian y-direction,
160 
161  - initialize rhs with flux differences (J1)
162  - add gravity (J4)
163  **************************************************** */
164 
165  x = x1[i];
166  z = x3[k];
167  for (j = beg; j <= end; j++) {
168  y = x2[j];
169  dtdx = dt/dx2[j];
170 
171  /* -----------------------------------------------
172  J1. initialize rhs with flux difference
173  ----------------------------------------------- */
174 
175  NVAR_LOOP(nv) rhs[j][nv] = -dtdx*(flux[j][nv] - flux[j-1][nv]);
176  #if USE_PR_GRADIENT == YES
177  rhs[j][MX2] -= dtdx*(p[j] - p[j-1]);
178  #endif
179 
180  /* ----------------------------------------------------
181  J4. Include gravity
182  ---------------------------------------------------- */
183 
184  v = state->vh[j];
185  #if (BODY_FORCE & VECTOR)
186  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
187  rhs[j][MX2] += dt*v[RHO]*g[JDIR];
188  #if DIMENSIONS == 2 && COMPONENTS == 3
189  rhs[j][MX3] += dt*vc[RHO]*g[KDIR];
190  #endif
191  #if HAVE_ENERGY
192  rhs[j][ENG] += dt*0.5*(flux[j][RHO] + flux[j-1][RHO])*g[JDIR];
193  #if DIMENSIONS == 2 && COMPONENTS == 3
194  rhs[j][ENG] += dt*vc[RHO]*vc[VX3]*g[KDIR];
195  #endif
196  #endif
197  #endif
198  }
199 
200  }else if (g_dir == KDIR){
201 
202  /* ****************************************************
203  Cartesian z-direction,
204 
205  - initialize rhs with flux differences (K1)
206  - enforce conservation of total angular
207  momentum and/or energy (K3)
208  - add gravity (K4)
209  **************************************************** */
210 
211  x = x1[i];
212  y = x2[j];
213  for (k = beg; k <= end; k++) {
214  z = x3[k];
215  dtdx = dt/dx3[k];
216 
217  /* -----------------------------------------------
218  K1. initialize rhs with flux difference
219  ----------------------------------------------- */
220 
221  NVAR_LOOP(nv) rhs[k][nv] = -dtdx*(flux[k][nv] - flux[k-1][nv]);
222  #if USE_PR_GRADIENT == YES
223  rhs[k][MX3] -= dtdx*(p[k] - p[k-1]);
224  #endif
225 
226  /* ----------------------------------------------------
227  K4. Include gravity
228  ---------------------------------------------------- */
229 
230  v = state->vh[k];
231  #if (BODY_FORCE & VECTOR)
232  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
233  rhs[k][MX3] += dt*v[RHO]*g[KDIR];
234  #if HAVE_ENERGY
235  rhs[k][ENG] += dt*0.5*(flux[k][RHO] + flux[k-1][RHO])*g[KDIR];
236  #endif
237  #endif
238  }
239  }
240 }
241 #elif GEOMETRY == CYLINDRICAL
242 
243 /* ***********************************************************
244 
245  Compute right-hand side of the RMHD/RHD
246  equations in CYLINDRICAL geometry.
247 
248  *********************************************************** */
249 {
250  double R, z, phi;
251  double R1, R_1;
252  double vB, vel2, lor2, wt, mphi, B2;
253 
254  if (g_dir == IDIR) {
255 
256  /* ****************************************************
257  Cylindrical radial direction:
258  multiply fluxes times interface area
259  **************************************************** */
260 
261  z = x2[g_j];
262  phi = 0.0;
263  for (i = beg - 1; i <= end; i++){
264  R = grid[IDIR].A[i];
265 
266  fA[i][RHO] = flux[i][RHO]*R;
267  EXPAND(fA[i][iMR] = flux[i][iMR]*R; ,
268  fA[i][iMZ] = flux[i][iMZ]*R; ,
269  fA[i][iMPHI] = flux[i][iMPHI]*R*R;)
270  #if PHYSICS == RMHD
271  EXPAND(fA[i][iBR] = flux[i][iBR]*R; ,
272  fA[i][iBZ] = flux[i][iBZ]*R; ,
273  fA[i][iBPHI] = flux[i][iBPHI]*R;)
274  #endif
275  #if HAVE_ENERGY
276  fA[i][ENG] = flux[i][ENG]*R;
277  #endif
278  #ifdef GLM_MHD
279  fA[i][PSI_GLM] = flux[i][PSI_GLM]*R;
280  #endif
281  NSCL_LOOP(nv) fA[i][nv] = flux[i][nv]*R;
282  }
283 
284  /* ****************************************************
285  Cylindrical radial direction,
286 
287  - initialize rhs with flux differences (I1)
288  - add source terms (I2)
289  - add gravity (I4)
290  **************************************************** */
291 
292  #if COMPONENTS == 3
293  Enthalpy (state->v, h, beg, end);
294  #endif
295  for (i = beg; i <= end; i++){
296  R = x1[i];
297  dtdV = dt/dV1[i];
298  dtdx = dt/dx1[i];
299  R_1 = 1.0/R;
300 
301  /* -----------------------------------------------
302  I1. initialize rhs with flux difference
303  ----------------------------------------------- */
304 
305  rhs[i][RHO] = -dtdV*(fA[i][RHO] - fA[i-1][RHO]);
306  EXPAND(rhs[i][iMR] = - dtdV*(fA[i][iMR] - fA[i-1][iMR]); ,
307  rhs[i][iMZ] = - dtdV*(fA[i][iMZ] - fA[i-1][iMZ]); ,
308  rhs[i][iMPHI] = - dtdV*(fA[i][iMPHI] - fA[i-1][iMPHI])*fabs(R_1);)
309  #if USE_PR_GRADIENT == YES
310  rhs[i][iMR] -= dtdx*(p[i] - p[i-1]);
311  #endif
312  #if PHYSICS == RMHD
313  EXPAND(rhs[i][iBR] = - dtdV*(fA[i][iBR] - fA[i-1][iBR]); ,
314  rhs[i][iBZ] = - dtdV*(fA[i][iBZ] - fA[i-1][iBZ]); ,
315  rhs[i][iBPHI] = - dtdV*(fA[i][iBPHI] - fA[i-1][iBPHI]);)
316  #ifdef GLM_MHD
317  rhs[i][iBR] = - dtdx*(flux[i][iBR] - flux[i-1][iBR]);
318  rhs[i][PSI_GLM] = - dtdV*(fA[i][PSI_GLM] - fA[i-1][PSI_GLM]);
319  #endif
320  #endif
321  #if HAVE_ENERGY
322  rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);
323  #endif
324  NSCL_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
325 
326 
327  /* -------------------------------------------------------
328  I2. Add source terms
329  [for Bphi we use the non-conservative formulation
330  with the source term since it has been found to
331  be more stable at low resolution in toroidal jet
332  simulations]
333  ------------------------------------------------------- */
334 
335  v = state->vh[i];
336  #if COMPONENTS == 3
337  vel2 = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
338  lor2 = 1.0/(1.0 - vel2);
339  #if PHYSICS == RMHD
340  vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
341  B2 = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]);
342  wt = v[RHO]*h[i]*lor2 + B2;
343  mphi = wt*v[iVPHI] - vB*v[iBPHI];
344  #elif PHYSICS == RHD
345  wt = v[RHO]*h[i]*lor2;
346  mphi = wt*v[iVPHI];
347  #endif
348 
349  rhs[i][iMR] += dt*mphi*v[iVPHI]*R_1;
350  #if PHYSICS == RMHD
351  rhs[i][iMR] -= dt*(v[iBPHI]/lor2 + vB*v[iVPHI])*v[iBPHI]*R_1;
352  rhs[i][iBPHI] -= dt*(v[iVPHI]*v[iBR] - v[iBPHI]*v[iVR])*R_1;
353  #endif
354  #endif /* COMPONENTS == 3 */
355 
356  /* ----------------------------------------------------
357  I4. Include gravity
358  ---------------------------------------------------- */
359 
360  #if (BODY_FORCE & VECTOR)
361  BodyForceVector(v, g, x1[i], x2[g_j], x3[g_k]);
362  rhs[i][iMR] += dt*v[RHO]*g[g_dir];
363  #if HAVE_ENERGY
364  rhs[i][ENG] += dt*0.5*(flux[i][RHO] + flux[i-1][RHO])*g[g_dir];
365  #endif
366  #endif
367  }
368 
369  } else if (g_dir == JDIR) {
370 
371  /* ****************************************************
372  Cylindrical vertical direction:
373 
374  - initialize rhs with flux differences (J1)
375  - add gravity (J4)
376  **************************************************** */
377 
378  R = x1[i];
379  phi = 0.0;
380  for (j = beg; j <= end; j++){
381  z = x2[j];
382  dtdx = dt/dx2[j];
383 
384  /* -----------------------------------------------
385  J1. initialize rhs with flux difference
386  ----------------------------------------------- */
387 
388  NVAR_LOOP(nv) rhs[j][nv] = -dtdx*(flux[j][nv] - flux[j-1][nv]);
389  #if USE_PR_GRADIENT == YES
390  rhs[j][iMZ] += - dtdx*(p[j] - p[j-1]);
391  #endif
392 
393  /* ----------------------------------------------------
394  J4. Include gravity
395  ---------------------------------------------------- */
396 
397  v = state->vh[j];
398  #if (BODY_FORCE & VECTOR)
399  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
400  rhs[j][iMZ] += dt*v[RHO]*g[JDIR];
401  #if HAVE_ENERGY
402  rhs[j][ENG] += dt*0.5*(flux[j][RHO] + flux[j-1][RHO])*g[JDIR];
403  #endif
404  #endif
405  }
406  }
407 }
408 
409 #elif GEOMETRY == POLAR
410 
411 /* ***********************************************************
412 
413  Compute right-hand side of the RMHD/RHD
414  equations in POLAR geometry.
415 
416  *********************************************************** */
417 {
418  double R, phi, z;
419  double R_1;
420  double vB, vel2, g_2, wt, mphi, B2;
421 
422  if (g_dir == IDIR) {
423 
424  /* ****************************************************
425  Polar radial direction:
426  multiply fluxes times interface area
427  **************************************************** */
428 
429  phi = x2[j];
430  z = x3[k];
431  for (i = beg - 1; i <= end; i++) {
432  R = grid[IDIR].A[i];
433 
434  fA[i][RHO] = flux[i][RHO]*R;
435  EXPAND(fA[i][iMR] = flux[i][iMR]*R; ,
436  fA[i][iMPHI] = flux[i][iMPHI]*R*R; ,
437  fA[i][iMZ] = flux[i][iMZ]*R;)
438  #if PHYSICS == RMHD
439  EXPAND(fA[i][iBR] = flux[i][iBR]*R; ,
440  fA[i][iBPHI] = flux[i][iBPHI]; ,
441  fA[i][iBZ] = flux[i][iBZ]*R;)
442  #endif
443  #if HAVE_ENERGY
444  fA[i][ENG] = flux[i][ENG]*R;
445  #endif
446  #ifdef GLM_MHD
447  fA[i][PSI_GLM] = flux[i][PSI_GLM]*R;
448  #endif
449  NSCL_LOOP(nv) fA[i][nv] = flux[i][nv]*R;
450  }
451 
452  /* ****************************************************
453  Polar radial direction,
454 
455  - initialize rhs with flux differences (I1)
456  - add source terms (I2)
457  - enforce conservation of total angular
458  momentum and/or energy (I3)
459  - add gravity (I4)
460  **************************************************** */
461 
462  #if COMPONENTS >= 2 /* -- need enthalpy for source term computation -- */
463  Enthalpy (state->v, h, beg, end);
464  #endif
465  for (i = beg; i <= end; i++) {
466  R = x1[i];
467  dtdV = dt/dV1[i];
468  dtdx = dt/dx1[i];
469  R_1 = grid[IDIR].r_1[i];
470 
471  /* -----------------------------------------------
472  I1. initialize rhs with flux difference
473  ----------------------------------------------- */
474 
475  rhs[i][RHO] = -dtdV*(fA[i][RHO] - fA[i-1][RHO]);
476  EXPAND(rhs[i][iMR] = - dtdV*(fA[i][iMR] - fA[i-1][iMR])
477  - dtdx*(p[i] - p[i-1]); ,
478  rhs[i][iMPHI] = - dtdV*(fA[i][iMPHI] - fA[i-1][iMPHI])*R_1; ,
479  rhs[i][iMZ] = - dtdV*(fA[i][iMZ] - fA[i-1][iMZ]);)
480  #if PHYSICS == RMHD
481  EXPAND(rhs[i][iBR] = -dtdV*(fA[i][iBR] - fA[i-1][iBR]); ,
482  rhs[i][iBPHI] = -dtdx*(fA[i][iBPHI] - fA[i-1][iBPHI]); ,
483  rhs[i][iBZ] = -dtdV*(fA[i][iBZ] - fA[i-1][iBZ]);)
484  #ifdef GLM_MHD
485  rhs[i][iBR] = -dtdx*(flux[i][iBR] - flux[i-1][iBR]);
486  rhs[i][PSI_GLM] = -dtdV*(fA[i][PSI_GLM] - fA[i-1][PSI_GLM]);
487  #endif
488  #endif
489  #if HAVE_ENERGY
490  rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);
491  #endif
492 
493  NSCL_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
494 
495 
496  /* ----------------------------------------------------
497  I2. Add source terms to the radial momentum eqn.
498  S = w*gamma^2*v(phi)^2 - B(phi)^2 - E(phi)^2
499  ---------------------------------------------------- */
500 
501  v = state->vh[i];
502  #if COMPONENTS >= 2
503  vel2 = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
504  g_2 = 1.0 - vel2;
505  #if PHYSICS == RMHD
506  vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
507  B2 = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]);
508  wt = v[RHO]*h[i]/g_2 + B2;
509  mphi = wt*v[iVPHI] - vB*v[iBPHI];
510  #elif PHYSICS == RHD
511  wt = v[RHO]*h[i]/g_2;
512  mphi = wt*v[iVPHI];
513  #endif
514 
515  rhs[i][iMR] += dt*mphi*v[iVPHI]*R_1;
516  #if PHYSICS == RMHD
517  rhs[i][iMR] -= dt*(v[iBPHI]*g_2 + vB*v[iVPHI])*v[iBPHI]*R_1;
518  #endif
519  #endif /* COMPONENTS >= 2 */
520 
521  /* ----------------------------------------------------
522  I4. Include gravity
523  ---------------------------------------------------- */
524 
525  #if (BODY_FORCE & VECTOR)
526  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
527  rhs[i][iMR] += dt*v[RHO]*g[IDIR];
528  #if HAVE_ENERGY
529  rhs[i][ENG] += dt*0.5*(flux[i][RHO] + flux[i-1][RHO])*g[JDIR];
530  #endif
531  #endif
532  }
533 
534  } else if (g_dir == JDIR) {
535 
536  /* ****************************************************
537  Polar azimuthal direction:
538 
539  - initialize rhs with flux differences (J1)
540  - add gravity (J4)
541  **************************************************** */
542 
543  R = x1[i];
544  z = x3[k];
545  scrh = dt/R;
546  for (j = beg; j <= end; j++){
547  phi = x2[j];
548  dtdx = scrh/dx2[j];
549 
550  /* ------------------------------------------------
551  J1. Compute equations rhs for phi-contributions
552  ------------------------------------------------ */
553 
554  NVAR_LOOP(nv) rhs[j][nv] = -dtdx*(flux[j][nv] - flux[j-1][nv]);
555  rhs[j][iMPHI] -= dtdx*(p[j] - p[j-1]);
556 
557  /* -------------------------------------------------------
558  J4. Include gravity
559  ------------------------------------------------------- */
560 
561  v = state->vh[j];
562  #if (BODY_FORCE & VECTOR)
563  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
564  rhs[j][iMPHI] += dt*v[RHO]*g[JDIR];
565  #if HAVE_ENERGY
566  rhs[j][ENG] += dt*0.5*(flux[j][RHO] + flux[j-1][RHO])*g[JDIR];
567  #endif
568  #endif
569  }
570 
571  } else if (g_dir == KDIR) {
572 
573  /* ****************************************************
574  Polar vertical direction:
575 
576  - initialize rhs with flux differences (K1)
577  - enforce conservation of total angular
578  momentum and/or energy (K3)
579  - add gravity (K4)
580  **************************************************** */
581 
582  R = x1[i];
583  phi = x2[j];
584  for (k = beg; k <= end; k++){
585  z = x3[k];
586  dtdx = dt/dx3[k];
587 
588  /* -----------------------------------------------
589  K1. initialize rhs with flux difference
590  ----------------------------------------------- */
591 
592  NVAR_LOOP(nv) rhs[k][nv] = -dtdx*(flux[k][nv] - flux[k-1][nv]);
593  rhs[k][iMZ] -= dtdx*(p[k] - p[k-1]);
594 
595  /* ----------------------------------------------------
596  K4. Include gravity
597  ---------------------------------------------------- */
598 
599  v = state->vh[k];
600  #if (BODY_FORCE & VECTOR)
601  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
602  rhs[k][iMZ] += dt*v[RHO]*g[KDIR];
603  #if HAVE_ENERGY
604  rhs[k][ENG] += dt*0.5*(flux[k][RHO] + flux[k-1][RHO])*g[KDIR];
605  #endif
606  #endif
607  }
608  }
609 }
610 #elif GEOMETRY == SPHERICAL
611 
612 /* ***********************************************************
613 
614  Compute right-hand side of the RMHD/RHD
615  equations in SPHERICAL geometry.
616 
617  *********************************************************** */
618 {
619  double r, th, phi;
620  double r2, r3, r_1;
621  double s, s2, ct, s_1;
622  double vB, vel2, lor2, wt, B2;
623  double mth = 0.0, mphi = 0.0;
624 
625  if (g_dir == IDIR) {
626  double Sm;
627 
628  /* ****************************************************
629  Spherical radial direction:
630  multiply fluxes by interface area
631  **************************************************** */
632 
633  th = x2[j]; s = sin(th);
634  phi = x3[k];
635  for (i = beg - 1; i <= end; i++){
636  r = x1p[i];
637  r2 = r*r;
638  r3 = r2*r;
639 
640  fA[i][RHO] = flux[i][RHO]*r2;
641  EXPAND(fA[i][iMR] = flux[i][iMR]*r2; ,
642  fA[i][iMTH] = flux[i][iMTH]*r2; ,
643  fA[i][iMPHI] = flux[i][iMPHI]*r3;)
644  #if PHYSICS == RMHD
645  EXPAND(fA[i][iBR] = flux[i][iBR]*r2; ,
646  fA[i][iBTH] = flux[i][iBTH]*r; ,
647  fA[i][iBPHI] = flux[i][iBPHI]*r;)
648  #endif
649  #if HAVE_ENERGY
650  fA[i][ENG] = flux[i][ENG]*r2;
651  #endif
652  #ifdef GLM_MHD
653  fA[i][PSI_GLM] = flux[i][PSI_GLM]*r2;
654  #endif
655  NSCL_LOOP(nv) fA[i][nv] = flux[i][nv]*r2;
656  }
657 
658  /* ****************************************************
659  Spherical radial direction:
660 
661  - initialize rhs with flux differences (I1)
662  - add source terms (I2)
663  - enforce conservation of total angular
664  momentum and/or energy (I3)
665  - add gravity (I4)
666  **************************************************** */
667 
668  #if COMPONENTS >= 2 /* -- need enthalpy for source term computation -- */
669  Enthalpy (state->v, h, beg, end);
670  #endif
671  for (i = beg; i <= end; i++) {
672  r = x1[i];
673  dtdV = dt/dV1[i];
674  dtdx = dt/dx1[i];
675  r_1 = grid[IDIR].r_1[i];
676 
677  /* -----------------------------------------------
678  I1. initialize rhs with flux difference
679  ----------------------------------------------- */
680 
681  rhs[i][RHO] = -dtdV*(fA[i][RHO] - fA[i-1][RHO]);
682  EXPAND(
683  rhs[i][iMR] = - dtdV*(fA[i][iMR] - fA[i-1][iMR])
684  - dtdx*(p[i] - p[i-1]); ,
685  rhs[i][iMTH] = -dtdV*(fA[i][iMTH] - fA[i-1][iMTH]); ,
686  rhs[i][iMPHI] = -dtdV*(fA[i][iMPHI] - fA[i-1][iMPHI])*r_1;
687  )
688  #if PHYSICS == RMHD
689  EXPAND(
690  rhs[i][iBR] = -dtdV*(fA[i][iBR] - fA[i-1][iBR]); ,
691  rhs[i][iBTH] = -dtdx*(fA[i][iBTH] - fA[i-1][iBTH])*r_1; ,
692  rhs[i][iBPHI] = -dtdx*(fA[i][iBPHI] - fA[i-1][iBPHI])*r_1;
693  )
694  #ifdef GLM_MHD
695  rhs[i][iBR] = -dtdx*(flux[i][iBR] - flux[i-1][iBR]);
696  rhs[i][PSI_GLM] = -dtdV*(fA[i][PSI_GLM] - fA[i-1][PSI_GLM]);
697  #endif
698  #endif
699  #if HAVE_ENERGY
700  rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);
701  #endif
702 
703 
704  NSCL_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
705 
706 
707  /* ----------------------------------------------------
708  I2. Add source terms
709  ---------------------------------------------------- */
710 
711  v = state->vh[i];
712 
713  vel2 = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
714  lor2 = 1.0/(1.0 - vel2);
715  #if PHYSICS == RMHD
716  vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
717  B2 = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]);
718  wt = v[RHO]*h[i]*lor2 + B2;
719  EXPAND( ,
720  mth = wt*v[iVTH] - vB*v[iBTH]; ,
721  mphi = wt*v[iVPHI] - vB*v[iBPHI];)
722  #elif PHYSICS == RHD
723  wt = v[RHO]*h[i]*lor2;
724  EXPAND( ; ,
725  mth = wt*v[iVTH]; ,
726  mphi = wt*v[iVPHI];)
727  #endif
728 
729  Sm = EXPAND( 0.0, + mth*v[iVTH], + mphi*v[iVPHI]);
730  #if PHYSICS == RMHD
731  Sm += EXPAND( 0.0, - (v[iBTH]/lor2 + vB*v[iVTH])*v[iBTH],
732  - (v[iBPHI]/lor2 + vB*v[iVPHI])*v[iBPHI]);
733  #endif
734  rhs[i][iMR] += dt*Sm*r_1;
735 
736  /* ----------------------------------------------------
737  I4. Include gravity
738  ---------------------------------------------------- */
739 
740  #if (BODY_FORCE & VECTOR)
741  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
742  rhs[i][iMR] += dt*v[RHO]*g[IDIR];
743  #if HAVE_ENERGY
744  rhs[i][ENG] += dt*0.5*(flux[i][RHO] + flux[i-1][RHO])*g[IDIR];
745  #endif
746  #endif
747  }
748 
749  } else if (g_dir == JDIR) {
750  double Sm;
751 
752  /* ****************************************************
753  Spherical meridional direction:
754  multiply fluxes by zone-interface area
755  **************************************************** */
756 
757  r = x1[i];
758  phi = x3[k];
759  for (j = beg - 1; j <= end; j++){
760  s = grid[JDIR].A[j];
761  s2 = s*s;
762 
763  fA[j][RHO] = flux[j][RHO]*s;
764  EXPAND(fA[j][iMR] = flux[j][iMR]*s; ,
765  fA[j][iMTH] = flux[j][iMTH]*s; ,
766  fA[j][iMPHI] = flux[j][iMPHI]*s2;)
767  #if PHYSICS == RMHD
768  EXPAND(fA[j][iBR] = flux[j][iBR]*s; ,
769  fA[j][iBTH] = flux[j][iBTH]*s; ,
770  fA[j][iBPHI] = flux[j][iBPHI];)
771  #endif
772  #if HAVE_ENERGY
773  fA[j][ENG] = flux[j][ENG]*s;
774  #endif
775  #ifdef GLM_MHD
776  fA[j][PSI_GLM] = flux[j][PSI_GLM]*s;
777  #endif
778  NSCL_LOOP(nv) fA[j][nv] = flux[j][nv]*s;
779  }
780 
781  /* ****************************************************
782  Spherical meridional direction:
783 
784  - initialize rhs with flux differences (J1)
785  - add source terms (J2)
786  - enforce conservation of total angular
787  momentum and/or energy (J3)
788  - add gravity (J4)
789  **************************************************** */
790 
791  r_1 = 0.5*(x1p[i]*x1p[i] - x1p[i-1]*x1p[i-1])/dV1[i];
792 
793  #if COMPONENTS >= 2 /* -- need enthalpy for source term computation -- */
794  Enthalpy (state->v, h, beg, end);
795  #endif
796  for (j = beg; j <= end; j++){
797  th = x2[j];
798  dtdV = dt/dV2[j]*r_1;
799  dtdx = dt/dx2[j]*r_1;
800  s = sin(th);
801  s_1 = 1.0/s;
802  ct = grid[JDIR].ct[j]; /* = cot(theta) */
803 
804  /* -----------------------------------------------
805  J1. initialize rhs with flux difference
806  ----------------------------------------------- */
807 
808  rhs[j][RHO] = -dtdV*(fA[j][RHO] - fA[j-1][RHO]);
809  EXPAND(
810  rhs[j][iMR] = - dtdV*(fA[j][iMR] - fA[j-1][iMR]); ,
811  rhs[j][iMTH] = - dtdV*(fA[j][iMTH] - fA[j-1][iMTH])
812  - dtdx*(p[j] - p[j-1]); ,
813  rhs[j][iMPHI] = - dtdV*(fA[j][iMPHI] - fA[j-1][iMPHI])*fabs(s_1);
814  )
815  #if PHYSICS == RMHD
816  EXPAND(
817  rhs[j][iBR] = -dtdV*(fA[j][iBR] - fA[j-1][iBR]); ,
818  rhs[j][iBTH] = -dtdV*(fA[j][iBTH] - fA[j-1][iBTH]); ,
819  rhs[j][iBPHI] = -dtdx*(fA[j][iBPHI] - fA[j-1][iBPHI]);
820  )
821  #ifdef GLM_MHD
822  rhs[j][iBTH] = -dtdx*(flux[j][iBTH] - flux[j-1][iBTH]);
823  rhs[j][PSI_GLM] = -dtdV*(fA[j][PSI_GLM] - fA[j-1][PSI_GLM]);
824  #endif
825  #endif
826  #if HAVE_ENERGY
827  rhs[j][ENG] = -dtdV*(fA[j][ENG] - fA[j-1][ENG]);
828  #endif
829 
830  NSCL_LOOP(nv) rhs[j][nv] = -dtdV*(fA[j][nv] - fA[j-1][nv]);
831 
832 
833  /* ----------------------------------------------------
834  J2. Add source terms
835  ---------------------------------------------------- */
836 
837  v = state->vh[j];
838  vel2 = EXPAND(v[VX1]*v[VX1], + v[VX2]*v[VX2], + v[VX3]*v[VX3]);
839  lor2 = 1.0/(1.0 - vel2);
840  #if PHYSICS == RMHD
841  vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
842  B2 = EXPAND(v[BX1]*v[BX1], + v[BX2]*v[BX2], + v[BX3]*v[BX3]);
843  wt = v[RHO]*h[j]*lor2 + B2;
844  EXPAND( ,
845  mth = wt*v[iVTH] - vB*v[iBTH]; ,
846  mphi = wt*v[iVPHI] - vB*v[iBPHI];)
847  #elif PHYSICS == RHD
848  wt = v[RHO]*h[j]*lor2;
849  EXPAND( ; ,
850  mth = wt*v[iVTH]; ,
851  mphi = wt*v[iVPHI];)
852  #endif
853 
854  Sm = EXPAND( 0.0, - mth*v[iVR], + ct*mphi*v[iVPHI]);
855  #if PHYSICS == RMHD
856  Sm += EXPAND( 0.0, + (v[iBTH]/lor2 + vB*v[iVTH])*v[iBR],
857  - ct*(v[iBPHI]/lor2 + vB*v[iVPHI])*v[iBPHI]);
858  #endif
859  rhs[j][iMTH] += dt*Sm*r_1;
860 
861  /* ----------------------------------------------------
862  J4. Include gravity
863  ---------------------------------------------------- */
864 
865  #if (BODY_FORCE & VECTOR)
866  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
867  rhs[j][iMTH] += dt*v[RHO]*g[JDIR];
868  #if HAVE_ENERGY
869  rhs[j][ENG] += dt*0.5*(flux[j][RHO] + flux[j-1][RHO])*g[JDIR];
870  #endif
871  #endif
872  }
873 
874  } else if (g_dir == KDIR) {
875 
876  /* ****************************************************
877  Spherical azimuthal direction:
878 
879  - initialize rhs with flux differences (K1)
880  - add gravity (K4)
881  **************************************************** */
882 
883  r = x1[i];
884  th = x2[j];
885  r_1 = 0.5*(x1p[i]*x1p[i] - x1p[i-1]*x1p[i-1])/dV1[i];
886  scrh = dt*r_1*dx2[j]/dV2[j];
887 
888  for (k = beg; k <= end; k++) {
889  phi = x3[k];
890  dtdx = scrh/dx3[k];
891 
892  /* ------------------------------------------------
893  K1. initialize rhs with flux difference
894  ------------------------------------------------ */
895 
896  NVAR_LOOP(nv) rhs[k][nv] = -dtdx*(flux[k][nv] - flux[k-1][nv]);
897  rhs[k][iMPHI] -= dtdx*(p[k] - p[k-1]);
898 
899  /* -------------------------------------------------------
900  K4. Include gravity
901  ------------------------------------------------------- */
902 
903  v = state->vh[k];
904  #if (BODY_FORCE & VECTOR)
905  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
906  rhs[k][iMPHI] += dt*v[RHO]*g[KDIR];
907  #if HAVE_ENERGY
908  rhs[k][ENG] += dt*0.5*(flux[k][RHO] + flux[k-1][RHO])*g[KDIR];
909  #endif
910  #endif
911  }
912  }
913 }
914 #endif /* GEOMETRY == SPHERICAL */
915 
916 /* --------------------------------------------------
917  Powell's source terms
918  -------------------------------------------------- */
919 
920  #if PHYSICS == RMHD
921  #if DIVB_CONTROL == EIGHT_WAVES
922  for (i = beg; i <= end; i++) {
923  EXPAND(rhs[i][MX1] += dt*state->src[i][MX1]; ,
924  rhs[i][MX2] += dt*state->src[i][MX2]; ,
925  rhs[i][MX3] += dt*state->src[i][MX3];)
926 
927  EXPAND(rhs[i][BX1] += dt*state->src[i][BX1]; ,
928  rhs[i][BX2] += dt*state->src[i][BX2]; ,
929  rhs[i][BX3] += dt*state->src[i][BX3];)
930  #if HAVE_ENERGY
931  rhs[i][ENG] += dt*state->src[i][ENG];
932  #endif
933  }
934  #endif
935  #endif
936 
937 /* -------------------------------------------------
938  Extended GLM source terms
939  ------------------------------------------------- */
940 
941  #if (defined GLM_MHD) && (GLM_EXTENDED == YES)
942  print1 ("! RightHandSide(): Extended GLM source terms not defined for RMHD\n");
943  QUIT_PLUTO(1);
944  #endif
945 
946 /* --------------------------------------------------
947  Reset right hand side in internal boundary zones
948  -------------------------------------------------- */
949 
950  #if INTERNAL_BOUNDARY == YES
951  InternalBoundaryReset(state, Dts, beg, end, grid);
952  #endif
953 
954 /* --------------------------------------------------
955  Time step determination
956  -------------------------------------------------- */
957 
958 #if !GET_MAX_DT
959 return;
960 #endif
961 
962  cl = 0.0;
963  for (i = beg-1; i <= end; i++) {
964  scrh = Dts->cmax[i]*grid[g_dir].inv_dxi[i];
965  cl = MAX(cl, scrh);
966  }
967  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
968  if (g_dir == JDIR) {
969  cl /= fabs(grid[IDIR].xgc[g_i]);
970  }
971  #if GEOMETRY == SPHERICAL
972  if (g_dir == KDIR){
973  cl /= fabs(grid[IDIR].xgc[g_i])*sin(grid[JDIR].xgc[g_j]);
974  }
975  #endif
976  #endif
977  Dts->inv_dta = MAX(cl, Dts->inv_dta);
978 }
MX3
#define MX3
Definition: mod_defs.h:22
MAX
#define MAX(a, b)
Definition: macros.h:101
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
MX1
#define MX1
Definition: mod_defs.h:20
STATE_1D::vh
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
VX2
#define VX2
Definition: mod_defs.h:29
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
iMTH
#define iMTH
Definition: mod_defs.h:70
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
STATE_1D::rhs
double ** rhs
Conservative right hand side.
Definition: structs.h:163
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
YES
#define YES
Definition: pluto.h:25
Enthalpy
void Enthalpy(double **v, real *h, int beg, int end)
Definition: eos.c:48
iVPHI
#define iVPHI
Definition: mod_defs.h:67
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
iMPHI
#define iMPHI
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
B2
#define B2
Definition: cooling_ode_solver.c:362
GRID::dx
double * dx
Definition: structs.h:83
TIME_STEP::cmax
double * cmax
Maximum signal velocity for hyperbolic eqns.
Definition: structs.h:214
iBR
#define iBR
Definition: mod_defs.h:150
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
iMZ
#define iMZ
Definition: mod_defs.h:59
iMR
#define iMR
Definition: mod_defs.h:69
TIME_STEP::inv_dta
double inv_dta
Inverse of advection (hyperbolic) time step, .
Definition: structs.h:215
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
STATE_1D::src
double ** src
Definition: structs.h:160
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
NVAR_LOOP
#define NVAR_LOOP(n)
Definition: pluto.h:618
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
PHYSICS
#define PHYSICS
Definition: definitions_01.h:1
j
int j
Definition: analysis.c:2
iVTH
#define iVTH
Definition: mod_defs.h:66
TIME_STEP
Definition: structs.h:213
MX2
#define MX2
Definition: mod_defs.h:21
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
iBPHI
#define iBPHI
Definition: mod_defs.h:152
k
int k
Definition: analysis.c:2
MXn
int MXn
Definition: globals.h:74
GRID::x
double * x
Definition: structs.h:80
RightHandSide
void RightHandSide(const State_1D *state, Time_Step *Dts, int beg, int end, double dt, Grid *grid)
Definition: rhs.c:88
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
iBTH
#define iBTH
Definition: mod_defs.h:151
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
USE_PR_GRADIENT
#define USE_PR_GRADIENT
Definition: rhs.c:9
pluto.h
PLUTO main header file.
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
STATE_1D
Definition: structs.h:133
i
int i
Definition: analysis.c:2
InternalBoundaryReset
void InternalBoundaryReset(const State_1D *, Time_Step *, int, int, Grid *)
BodyForceVector
void BodyForceVector(double *, double *, double, double, double)
Definition: init.c:441
BX1
#define BX1
Definition: mod_defs.h:25
GRID::ct
double * ct
Geometrical factor cot(theta).
Definition: structs.h:89
VX3
#define VX3
Definition: mod_defs.h:30
RHD
#define RHD
Definition: pluto.h:110
STATE_1D::press
double * press
Upwind pressure term computed with the Riemann solver.
Definition: structs.h:164
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
GRID::inv_dxi
double * inv_dxi
inverse of the distance between the center of two cells, inv_dxi = .
Definition: structs.h:91
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
RMHD
#define RMHD
Definition: pluto.h:112
AdvectFlux
void AdvectFlux(const State_1D *state, int beg, int end, Grid *grid)
Definition: adv_flux.c:47
GLM_MHD
#define GLM_MHD
Definition: glm.h:25
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395