PLUTO
Macros | Functions
rhs.c File Reference

Compute the right hand side of the conservative HD/MHD equations. More...

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

Go to the source code of this file.

Macros

#define USE_PR_GRADIENT   YES /* -- default, do not change!! -- */
 
#define iMPHI   MX2 /* -- for Cartesian coordinates -- */
 

Functions

static void TotalFlux (const State_1D *, double *, int, int, Grid *)
 
void RightHandSide (const State_1D *state, Time_Step *Dts, int beg, int end, double dt, Grid *grid)
 

Detailed Description

Compute the right hand side of the conservative HD/MHD equations.

This function constructs the one-dimensional right hand side of the conservative MHD or HD equations in the direction given by g_dir in different geometries. The right hand side in the d direction is computed as a two-point flux difference term plus a source term:

\[ {\cal R}_{\ivec}^{(d)} = -\frac{\Delta t}{\Delta{\cal V}_{\ivec}}\ \Big[ (A{\cal F})_{\ivec+\HALF\hvec{e}_d} -(A{\cal F})_{\ivec-\HALF\hvec{e}_d} \Big] + \Delta t{\cal S}_{\ivec} \]

where $\ivec = (i,j,k)$ while

  • $ A $: interface areas
  • $ {\cal V} $: cell volume
  • $ \F $: interface fluxes
  • $ \Delta t $: time step
  • $ {\cal S} $: source term including geometrical terms and body forces.

See also the RHS_page. The right hand side is assembled through the following steps:

  • If either one of FARGO, ROTATION or gravitational potential is used, fluxes are combined to enforce conservation of total angular momentum and/or energy, see TotalFlux();
  • initialize rhs with flux differences;
  • add dissipative effects (viscosity and thermal conduction) to entropy equation. Ohmic dissipation is included in a separate step.

Source terms are added later in RightHandSideSource().

For the entropy equation, the dissipative contributions can be recovered from the internal energy equation, for which (Boyd, Eq. [3.27]):

\[ \pd{p}{t} + \vec{v}\cdot\nabla p + \Gamma p \nabla\cdot\vec{v} = (\Gamma-1)\left[\nabla\cdot\left(\kappa\nabla T\right) + \eta\vec{J}\cdot\vec{J} + \mu\left(\pd{v_i}{x_j}+\pd{v_j}{x_i} - \frac{2}{3}\delta_{ij}\nabla\cdot\vec{v}\right)\pd{v_i}{x_j}\right] \]

To obtain the corresponding contribution to the (conservative form of) entropy equation, just divide the previous one by $\rho^{\Gamma-1}$:

\[ \pd{(\rho s)}{t} + \nabla\cdot(\rho s\vec{v}) = (\Gamma-1)\rho^{1-\Gamma}\left[\nabla\cdot\left(\kappa\nabla T\right) + \eta\vec{J}^2 + \mu \Pi_{ij}\pd{v_i}{x_j}\right] \]

See the books by

  • Goedbloed & Poedts, "Principle of MHD", page 165-166
  • Boyd, "The Physics of Plasmas", page 57, Eq. 3.27

References

  • "PLUTO: A Numerical Code for Computational Astrophysics." Mignone et al, ApJS (2007) 170, 228
  • "The PLUTO Code for Adaptive Mesh Computations in Astrophysical Fluid Dynamics" Mignone et al, ApJS (2012) 198, 7M
Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Aug 20, 2015
Todo:
  • merge GetAreaFlux and TotalFlux ?

Definition in file rhs.c.

Macro Definition Documentation

#define iMPHI   MX2 /* -- for Cartesian coordinates -- */
#define USE_PR_GRADIENT   YES /* -- default, do not change!! -- */

Definition at line 83 of file rhs.c.

Function Documentation

void RightHandSide ( const State_1D *  state,
Time_Step *  Dts,
int  beg,
int  end,
double  dt,
Grid *  grid 
)
Parameters
[in,out]statepointer to State_1D structure
[in]Dtspointer to time step structure
[in]beginitial index of computation
[in]endfinal index of computation
[in]dttime increment
[in]gridpointer to Grid structure
Returns
This function has no return value.
Note
–
Todo:
–

Definition at line 88 of file rhs.c.

103 {
104  int i, j, k, nv;
105  double dtdx, dtdV, scrh, rhog;
106  double *x1 = grid[IDIR].x, *x2 = grid[JDIR].x, *x3 = grid[KDIR].x;
107  double *x1p = grid[IDIR].xr, *x2p = grid[JDIR].xr, *x3p = grid[KDIR].xr;
108  double *dx1 = grid[IDIR].dx, *dx2 = grid[JDIR].dx, *dx3 = grid[KDIR].dx;
109  double *dV1 = grid[IDIR].dV, *dV2 = grid[JDIR].dV, *dV3 = grid[KDIR].dV;
110  double **rhs = state->rhs;
111  double **flux = state->flux, *p = state->press;
112  double **vh = state->vh, **vp = state->vp, **vm = state->vm;
113  double *A = grid[g_dir].A, *dV = grid[g_dir].dV;
114 
115  double cl;
116  double w, wp, vphi, phi_c;
117  double g[3];
118  static double **fA, *phi_p;
119  static double **fvA;
120 #if ENTROPY_SWITCH
121  double rhs_entr;
122  double **visc_flux = state->visc_flux;
123  double **visc_src = state->visc_src;
124  double **tc_flux = state->tc_flux;
125  double **res_flux = state->res_flux;
126  if (fvA == NULL) fvA = ARRAY_2D(NMAX_POINT, NVAR, double);
127 #endif
128 
129  #if GEOMETRY != CARTESIAN
130  if (fA == NULL) fA = ARRAY_2D(NMAX_POINT, NVAR, double);
131  #endif
132  if (phi_p == NULL) phi_p = ARRAY_1D(NMAX_POINT, double);
133 
134 /* --------------------------------------------------
135  1. Compute fluxes for passive scalars and dust
136  -------------------------------------------------- */
137 
138 #if NSCL > 0
139  AdvectFlux (state, beg - 1, end, grid);
140 #endif
141 
142 #if DUST == YES
143  Dust_Solver(state, beg - 1, end, Dts->cmax, grid);
144 #endif
145 
146  i = g_i; /* will be redefined during x1-sweep */
147  j = g_j; /* will be redefined during x2-sweep */
148  k = g_k; /* will be redefined during x3-sweep */
149 
150 /* ------------------------------------------------
151  Add pressure to normal component of
152  momentum flux if necessary.
153  ------------------------------------------------ */
154 
155 #if USE_PR_GRADIENT == NO
156  for (i = beg - 1; i <= end; i++) flux[i][MXn] += p[i];
157 #endif
158 
159 /* -----------------------------------------------------
160  Compute gravitational potential at cell interfaces
161  ----------------------------------------------------- */
162 
163 #if (BODY_FORCE & POTENTIAL)
164  if (g_dir == IDIR) {
165  for (i = beg-1; i <= end; i++) {
166  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
167  }
168  }else if (g_dir == JDIR){
169  for (j = beg-1; j <= end; j++) {
170  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
171  }
172  }else if (g_dir == KDIR){
173  for (k = beg-1; k <= end; k++) {
174  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
175  }
176  }
177 #endif
178 
179 /* -----------------------------------------------------
180  Compute total flux
181  ----------------------------------------------------- */
182 
183  #if (defined FARGO && !defined SHEARINGBOX) ||\
184  (ROTATING_FRAME == YES) || (BODY_FORCE & POTENTIAL)
185  TotalFlux(state, phi_p, beg-1, end, grid);
186  #endif
187 
188 #if GEOMETRY == CARTESIAN
189 
190  if (g_dir == IDIR){
191 
192  for (i = beg; i <= end; i++) {
193  dtdx = dt/dx1[i];
194 
195  /* -- I1. initialize rhs with flux difference -- */
196 
197  NVAR_LOOP(nv) rhs[i][nv] = -dtdx*(flux[i][nv] - flux[i-1][nv]);
198  #if USE_PR_GRADIENT == YES
199  rhs[i][MX1] -= dtdx*(p[i] - p[i-1]);
200  #endif
201 
202  /* -- I5. Add dissipative terms to entropy equation -- */
203 
204  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
205  rhog = vh[i][RHO];
206  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
207 
208  rhs_entr = 0.0;
209  #if VISCOSITY == EXPLICIT
210  rhs_entr += (visc_flux[i][ENG] - visc_flux[i-1][ENG]);
211  rhs_entr -= EXPAND( vh[i][VX1]*(visc_flux[i][MX1] - visc_flux[i-1][MX1]) ,
212  + vh[i][VX2]*(visc_flux[i][MX2] - visc_flux[i-1][MX2]) ,
213  + vh[i][VX3]*(visc_flux[i][MX3] - visc_flux[i-1][MX3]));
214  #endif
215 
216  #if THERMAL_CONDUCTION == EXPLICIT
217  rhs_entr += (tc_flux[i][ENG] - tc_flux[i-1][ENG]);
218  #endif
219  rhs[i][ENTR] += rhs_entr*dtdx*rhog;
220  #endif
221 
222  }
223  } else if (g_dir == JDIR){
224 
225  for (j = beg; j <= end; j++) {
226  dtdx = dt/dx2[j];
227 
228  /* -- J1. initialize rhs with flux difference -- */
229 
230  NVAR_LOOP(nv) rhs[j][nv] = -dtdx*(flux[j][nv] - flux[j-1][nv]);
231  #if USE_PR_GRADIENT == YES
232  rhs[j][MX2] -= dtdx*(p[j] - p[j-1]);
233  #endif
234 
235  /* -- J5. Add dissipative terms to entropy equation -- */
236 
237  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
238  rhog = vh[j][RHO];
239  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
240 
241  rhs_entr = 0.0;
242  #if VISCOSITY == EXPLICIT
243  rhs_entr += (visc_flux[j][ENG] - visc_flux[j-1][ENG]);
244  rhs_entr -= EXPAND( vh[j][VX1]*(visc_flux[j][MX1] - visc_flux[j-1][MX1]) ,
245  + vh[j][VX2]*(visc_flux[j][MX2] - visc_flux[j-1][MX2]) ,
246  + vh[j][VX3]*(visc_flux[j][MX3] - visc_flux[j-1][MX3]));
247  #endif
248 
249  #if THERMAL_CONDUCTION == EXPLICIT
250  rhs_entr += (tc_flux[j][ENG] - tc_flux[j-1][ENG]);
251  #endif
252  rhs[j][ENTR] += rhs_entr*dtdx*rhog;
253  #endif
254 
255  }
256 
257  }else if (g_dir == KDIR){
258 
259  for (k = beg; k <= end; k++) {
260  dtdx = dt/dx3[k];
261 
262  /* -- K1. initialize rhs with flux difference -- */
263 
264  NVAR_LOOP(nv) rhs[k][nv] = -dtdx*(flux[k][nv] - flux[k-1][nv]);
265  #if USE_PR_GRADIENT == YES
266  rhs[k][MX3] -= dtdx*(p[k] - p[k-1]);
267  #endif
268 
269  /* -- K5. Add dissipative terms to entropy equation -- */
270 
271  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
272  rhog = vh[k][RHO];
273  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
274  rhs_entr = 0.0;
275  #if VISCOSITY == EXPLICIT
276  rhs_entr += (visc_flux[k][ENG] - visc_flux[k-1][ENG]);
277  rhs_entr -= EXPAND( vh[k][VX1]*(visc_flux[k][MX1] - visc_flux[k-1][MX1]) ,
278  + vh[k][VX2]*(visc_flux[k][MX2] - visc_flux[k-1][MX2]) ,
279  + vh[k][VX3]*(visc_flux[k][MX3] - visc_flux[k-1][MX3]));
280  #endif
281 
282  #if THERMAL_CONDUCTION == EXPLICIT
283  rhs_entr += (tc_flux[k][ENG] - tc_flux[k-1][ENG]);
284  #endif
285  rhs[k][ENTR] += rhs_entr*dtdx*rhog;
286  #endif
287  }
288  }
289 
290 #elif GEOMETRY == CYLINDRICAL
291 
292 {
293  double R, z, phi, R_1;
294 
295 #if DUST == YES
296  #error "DUST not implemented in CYLINDRICAL coordinates"
297 #endif
298 
299  if (g_dir == IDIR) {
300  double vc[NVAR];
301 
302  GetAreaFlux (state, fA, fvA, beg, end, grid);
303  for (i = beg; i <= end; i++){
304  R = x1[i];
305  dtdV = dt/dV1[i];
306  dtdx = dt/dx1[i];
307  R_1 = 1.0/R;
308 
309  /* -- I1. initialize rhs with flux difference -- */
310 
311  rhs[i][RHO] = -dtdV*(fA[i][RHO] - fA[i-1][RHO]);
312  EXPAND(rhs[i][iMR] = - dtdV*(fA[i][iMR] - fA[i-1][iMR]); ,
313  rhs[i][iMZ] = - dtdV*(fA[i][iMZ] - fA[i-1][iMZ]); ,
314  rhs[i][iMPHI] = - dtdV*(fA[i][iMPHI] - fA[i-1][iMPHI])*fabs(R_1);)
315  #if USE_PR_GRADIENT == YES
316  rhs[i][iMR] -= dtdx*(p[i] - p[i-1]);
317  #endif
318  #if PHYSICS == MHD
319  EXPAND(rhs[i][iBR] = - dtdV*(fA[i][iBR] - fA[i-1][iBR]); ,
320  rhs[i][iBZ] = - dtdV*(fA[i][iBZ] - fA[i-1][iBZ]); ,
321  rhs[i][iBPHI] = - dtdx*(flux[i][iBPHI] - flux[i-1][iBPHI]);)
322  #ifdef GLM_MHD
323  rhs[i][iBR] = - dtdx*(flux[i][iBR] - flux[i-1][iBR]);
324  rhs[i][PSI_GLM] = - dtdV*(fA[i][PSI_GLM] - fA[i-1][PSI_GLM]);
325  #endif
326  #endif
327  IF_ENERGY(rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);)
328  NSCL_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
329 
330  /* -- I5. Add dissipative terms to entropy equation -- */
331 
332  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
333  NVAR_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
334  rhog = vc[RHO];
335  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
336 
337  rhs_entr = 0.0;
338  #if VISCOSITY == EXPLICIT
339  rhs_entr += (fvA[i][ENG] - fvA[i-1][ENG]);
340  rhs_entr -= EXPAND( vc[VX1]*(fvA[i][MX1] - fvA[i-1][MX1]) ,
341  + vc[VX2]*(fvA[i][MX2] - fvA[i-1][MX2]) ,
342  + vc[VX3]*(fvA[i][MX3] - fvA[i-1][MX3])*fabs(R_1));
343 
344  rhs_entr -= EXPAND( vc[VX1]*visc_src[i][MX1],
345  + vc[VX2]*visc_src[i][MX2],
346  + vc[VX3]*visc_src[i][MX3]);
347  #endif
348 
349  #if THERMAL_CONDUCTION == EXPLICIT
350  rhs_entr += (A[i]*tc_flux[i][ENG] - A[i-1]*tc_flux[i-1][ENG]);
351  #endif
352  rhs[i][ENTR] += rhs_entr*dtdV*rhog;
353  #endif
354  }
355 
356  } else if (g_dir == JDIR) {
357 
358  for (j = beg; j <= end; j++){
359  dtdx = dt/dx2[j];
360 
361  /* -- J1. initialize rhs with flux difference -- */
362 
363  NVAR_LOOP(nv) rhs[j][nv] = -dtdx*(flux[j][nv] - flux[j-1][nv]);
364  #if USE_PR_GRADIENT == YES
365  rhs[j][iMZ] += - dtdx*(p[j] - p[j-1]);
366  #endif
367 
368  /* -- J5. Add dissipative terms to entropy equation -- */
369 
370  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
371  rhog = vh[j][RHO];
372  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
373 
374  rhs_entr = 0.0;
375  #if VISCOSITY == EXPLICIT
376  rhs_entr += (visc_flux[j][ENG] - visc_flux[j-1][ENG]);
377  rhs_entr -= EXPAND( vc[VX1]*(visc_flux[j][MX1] - visc_flux[j-1][MX1]) ,
378  + vc[VX2]*(visc_flux[j][MX2] - visc_flux[j-1][MX2]) ,
379  + vc[VX3]*(visc_flux[j][MX3] - visc_flux[j-1][MX3]));
380  #endif
381  #if THERMAL_CONDUCTION == EXPLICIT
382  rhs_entr += (tc_flux[j][ENG] - tc_flux[j-1][ENG]);
383  #endif
384  rhs[j][ENTR] += rhs_entr*dtdx*rhog;
385  #endif
386  }
387  }
388 }
389 
390 #elif GEOMETRY == POLAR
391 {
392  double R, phi, z;
393  double R_1;
394 
395  if (g_dir == IDIR) {
396  double vc[NVAR];
397 
398  GetAreaFlux (state, fA, fvA, beg, end, grid);
399  for (i = beg; i <= end; i++) {
400  R = x1[i];
401  dtdV = dt/dV1[i];
402  dtdx = dt/dx1[i];
403  R_1 = grid[IDIR].r_1[i];
404 
405  /* -- I1. initialize rhs with flux difference -- */
406 
407  rhs[i][RHO] = -dtdV*(fA[i][RHO] - fA[i-1][RHO]);
408  EXPAND(rhs[i][iMR] = - dtdV*(fA[i][iMR] - fA[i-1][iMR])
409  - dtdx*(p[i] - p[i-1]); ,
410  rhs[i][iMPHI] = - dtdV*(fA[i][iMPHI] - fA[i-1][iMPHI])*R_1; ,
411  rhs[i][iMZ] = - dtdV*(fA[i][iMZ] - fA[i-1][iMZ]);)
412  #if PHYSICS == MHD
413  EXPAND(rhs[i][iBR] = -dtdV*(fA[i][iBR] - fA[i-1][iBR]); ,
414  rhs[i][iBPHI] = -dtdx*(flux[i][iBPHI] - flux[i-1][iBPHI]); ,
415  rhs[i][iBZ] = -dtdV*(fA[i][iBZ] - fA[i-1][iBZ]);)
416  #ifdef GLM_MHD
417  rhs[i][iBR] = -dtdx*(flux[i][iBR] - flux[i-1][iBR]);
418  rhs[i][PSI_GLM] = -dtdV*(fA[i][PSI_GLM] - fA[i-1][PSI_GLM]);
419  #endif
420  #endif
421  IF_ENERGY(rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);)
422 
423  NSCL_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
424 
425  IF_DUST(NDUST_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
426  rhs[i][MX2_D] *= R_1;)
427 
428  /* -- I5. Add dissipative terms to entropy equation -- */
429 
430  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
431  for (nv = NVAR; nv--; ) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
432  rhog = vc[RHO];
433  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
434 
435  rhs_entr = 0.0;
436  #if VISCOSITY == EXPLICIT
437  rhs_entr += (fvA[i][ENG] - fvA[i-1][ENG]);
438  rhs_entr -= EXPAND( vc[VX1]*(fvA[i][iMR] - fvA[i-1][iMR]) ,
439  + vc[VX2]*(fvA[i][iMPHI] - fvA[i-1][iMPHI])*R_1 ,
440  + vc[VX3]*(fvA[i][iMZ] - fvA[i-1][iMZ]));
441 
442  rhs_entr -= EXPAND( vc[VX1]*visc_src[i][MX1],
443  + vc[VX2]*visc_src[i][MX2],
444  + vc[VX3]*visc_src[i][MX3]);
445  #endif
446 
447  #if THERMAL_CONDUCTION == EXPLICIT
448  rhs_entr += (A[i]*tc_flux[i][ENG] - A[i-1]*tc_flux[i-1][ENG]);
449  #endif
450  rhs[i][ENTR] += rhs_entr*dtdV*rhog;
451  #endif
452 
453  }
454 
455  } else if (g_dir == JDIR) {
456 
457  scrh = dt/x1[i];
458  for (j = beg; j <= end; j++){
459  dtdx = scrh/dx2[j];
460 
461  /* -- J1. Compute equations rhs for phi-contributions -- */
462 
463  NVAR_LOOP(nv) rhs[j][nv] = -dtdx*(flux[j][nv] - flux[j-1][nv]);
464  rhs[j][iMPHI] -= dtdx*(p[j] - p[j-1]);
465 
466  /* -- J5. Add dissipative terms to entropy equation -- */
467 
468  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
469  rhog = vh[j][RHO];
470  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
471 
472  rhs_entr = 0.0;
473  #if VISCOSITY == EXPLICIT
474  rhs_entr += (visc_flux[j][ENG] - visc_flux[j-1][ENG]);
475  rhs_entr -= EXPAND( vh[j][VX1]*(visc_flux[j][MX1] - visc_flux[j-1][MX1]) ,
476  + vh[j][VX2]*(visc_flux[j][MX2] - visc_flux[j-1][MX2]) ,
477  + vh[j][VX3]*(visc_flux[j][MX3] - visc_flux[j-1][MX3]));
478 
479  rhs_entr -= EXPAND( vh[j][VX1]*visc_src[j][MX1],
480  + vh[j][VX2]*visc_src[j][MX2],
481  + vh[j][VX3]*visc_src[j][MX3]);
482  #endif
483  #if THERMAL_CONDUCTION == EXPLICIT
484  rhs_entr += (tc_flux[j][ENG] - tc_flux[j-1][ENG]);
485  #endif
486  rhs[j][ENTR] += rhs_entr*dtdx*rhog;
487  #endif
488 
489  }
490 
491  } else if (g_dir == KDIR) {
492 
493  for (k = beg; k <= end; k++){
494  dtdx = dt/dx3[k];
495 
496  /* -- K1. initialize rhs with flux difference -- */
497 
498  VAR_LOOP(nv) rhs[k][nv] = -dtdx*(flux[k][nv] - flux[k-1][nv]);
499  rhs[k][MX3] -= dtdx*(p[k] - p[k-1]);
500 
501  /* -- K5. Add dissipative terms to entropy equation -- */
502 
503  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
504  rhog = vh[k][RHO];
505  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
506 
507  rhs_entr = 0.0;
508  #if VISCOSITY == EXPLICIT
509  rhs_entr += (visc_flux[k][ENG] - visc_flux[k-1][ENG]);
510  rhs_entr -= EXPAND( vh[j][VX1]*(visc_flux[k][MX1] - visc_flux[k-1][MX1]) ,
511  + vh[j][VX2]*(visc_flux[k][MX2] - visc_flux[k-1][MX2]) ,
512  + vh[j][VX3]*(visc_flux[k][MX3] - visc_flux[k-1][MX3]));
513  #endif
514 
515  #if THERMAL_CONDUCTION == EXPLICIT
516  rhs_entr += (tc_flux[k][ENG] - tc_flux[k-1][ENG]);
517  #endif
518  rhs[k][ENTR] += rhs_entr*dtdx*rhog;
519  #endif
520 
521  }
522  }
523 }
524 #elif GEOMETRY == SPHERICAL
525 {
526  double r_1, th, s, s_1;
527 
528  if (g_dir == IDIR) {
529  double vc[NVAR], dVdx;
530 
531  GetAreaFlux (state, fA, fvA, beg, end, grid);
532  for (i = beg; i <= end; i++) {
533  dtdV = dt/dV1[i];
534  dtdx = dt/dx1[i];
535  r_1 = grid[IDIR].r_1[i];
536 
537  /* -- I1. initialize rhs with flux difference -- */
538 
539 /* Alternative sequence
540 dVdx = dV1[i]/dx1[i];
541 NVAR_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
542 rhs[i][MX1] -= dtdx*(p[i] - p[i-1]);
543 rhs[i][MX3] *= r_1;
544 #if PHYSICS == MHD
545  EXPAND(
546  ,
547  rhs[i][iBTH] *= dVrdx; ,
548  rhs[i][iBPHI] *= dVrdx;
549  )
550  #ifdef GLM_MHD
551  rhs[i][iBR] *= dVdx;
552  #endif
553 #endif
554 IF_DUST(rhs[i][MX3_D] *= r_1;)
555 */
556 
557  rhs[i][RHO] = -dtdV*(fA[i][RHO] - fA[i-1][RHO]);
558  EXPAND(
559  rhs[i][iMR] = - dtdV*(fA[i][iMR] - fA[i-1][iMR])
560  - dtdx*(p[i] - p[i-1]); ,
561  rhs[i][iMTH] = -dtdV*(fA[i][iMTH] - fA[i-1][iMTH]); ,
562  rhs[i][iMPHI] = -dtdV*(fA[i][iMPHI] - fA[i-1][iMPHI])*r_1;
563  )
564  #if PHYSICS == MHD
565  EXPAND(
566  rhs[i][iBR] = -dtdV*(fA[i][iBR] - fA[i-1][iBR]); ,
567  rhs[i][iBTH] = -dtdx*(fA[i][iBTH] - fA[i-1][iBTH])*r_1; ,
568  rhs[i][iBPHI] = -dtdx*(fA[i][iBPHI] - fA[i-1][iBPHI])*r_1;
569  )
570  #ifdef GLM_MHD
571  rhs[i][iBR] = -dtdx*(flux[i][iBR] - flux[i-1][iBR]);
572  rhs[i][PSI_GLM] = -dtdV*(fA[i][PSI_GLM] - fA[i-1][PSI_GLM]);
573  #endif
574  #endif
575  IF_ENERGY(rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);)
576 
577  NSCL_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
578  IF_DUST(NDUST_LOOP(nv) rhs[i][nv] = -dtdV*(fA[i][nv] - fA[i-1][nv]);
579  rhs[i][MX3_D] *= r_1;)
580 
581  /* -- I5. Add dissipative terms to entropy equation -- */
582 
583  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
584  NVAR_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
585 
586  rhog = vc[RHO];
587  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
588 
589  rhs_entr = 0.0;
590  #if VISCOSITY == EXPLICIT
591  rhs_entr += (fvA[i][ENG] - fvA[i-1][ENG]);
592  rhs_entr -= EXPAND( vc[VX1]*(fvA[i][iMR] - fvA[i-1][iMR]) ,
593  + vc[VX2]*(fvA[i][iMTH] - fvA[i-1][iMTH]) ,
594  + vc[VX3]*(fvA[i][iMPHI] - fvA[i-1][iMPHI])*r_1);
595 
596  rhs_entr -= EXPAND( vc[VX1]*visc_src[i][MX1],
597  + vc[VX2]*visc_src[i][MX2],
598  + vc[VX3]*visc_src[i][MX3]);
599  #endif
600  #if THERMAL_CONDUCTION == EXPLICIT
601  rhs_entr += (A[i]*tc_flux[i][ENG] - A[i-1]*tc_flux[i-1][ENG]);
602  #endif
603  rhs[i][ENTR] += rhs_entr*dtdV*rhog;
604  #endif
605  }
606 
607  } else if (g_dir == JDIR) {
608 
609  GetAreaFlux (state, fA, fvA, beg, end, grid);
610  r_1 = 0.5*(x1p[i]*x1p[i] - x1p[i-1]*x1p[i-1])/dV1[i];
611  scrh = dt*r_1;
612  for (j = beg; j <= end; j++){
613  th = x2[j];
614  dtdV = scrh/dV2[j];
615  dtdx = scrh/dx2[j];
616  s = sin(th);
617  s_1 = 1.0/s;
618 
619  /* -- J1. initialize rhs with flux difference -- */
620 
621  rhs[j][RHO] = -dtdV*(fA[j][RHO] - fA[j-1][RHO]);
622  EXPAND(
623  rhs[j][iMR] = - dtdV*(fA[j][iMR] - fA[j-1][iMR]); ,
624  rhs[j][iMTH] = - dtdV*(fA[j][iMTH] - fA[j-1][iMTH])
625  - dtdx*(p[j] - p[j-1]); ,
626  rhs[j][iMPHI] = - dtdV*(fA[j][iMPHI] - fA[j-1][iMPHI])*fabs(s_1);
627  )
628  #if PHYSICS == MHD
629  EXPAND(
630  rhs[j][iBR] = -dtdV*(fA[j][iBR] - fA[j-1][iBR]); ,
631  rhs[j][iBTH] = -dtdV*(fA[j][iBTH] - fA[j-1][iBTH]); ,
632  rhs[j][iBPHI] = -dtdx*(flux[j][iBPHI] - flux[j-1][iBPHI]);
633  )
634  #ifdef GLM_MHD
635  rhs[j][iBTH] = -dtdx*(flux[j][iBTH] - flux[j-1][iBTH]);
636  rhs[j][PSI_GLM] = -dtdV*(fA[j][PSI_GLM] - fA[j-1][PSI_GLM]);
637  #endif
638  #endif
639  IF_ENERGY(rhs[j][ENG] = -dtdV*(fA[j][ENG] - fA[j-1][ENG]);)
640 
641  NSCL_LOOP(nv) rhs[j][nv] = -dtdV*(fA[j][nv] - fA[j-1][nv]);
642 
643  IF_DUST(NDUST_LOOP(nv) rhs[j][nv] = -dtdV*(fA[j][nv] - fA[j-1][nv]);
644  rhs[j][MX3_D] *= fabs(s_1);)
645 
646  /* -- J5. Add TC dissipative term to entropy equation -- */
647 
648  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
649  rhog = vh[j][RHO];
650  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
651 
652  rhs_entr = 0.0;
653  #if VISCOSITY == EXPLICIT
654  rhs_entr += (fvA[j][ENG] - fvA[j-1][ENG]);
655  rhs_entr -= EXPAND( vc[VX1]*(fvA[j][iMR] - fvA[j-1][iMR]) ,
656  + vc[VX2]*(fvA[j][iMTH] - fvA[j-1][iMTH]) ,
657  + vc[VX3]*(fvA[j][iMPHI] - fvA[j-1][iMPHI])*fabs(s_1));
658 
659  rhs_entr -= EXPAND( vc[VX1]*visc_src[j][MX1],
660  + vc[VX2]*visc_src[j][MX2],
661  + vc[VX3]*visc_src[j][MX3]);
662  #endif
663  #if THERMAL_CONDUCTION == EXPLICIT
664  rhs_entr += (A[j]*tc_flux[j][ENG] - A[j-1]*tc_flux[j-1][ENG]);
665  #endif
666  rhs[j][ENTR] += rhs_entr*dtdV*rhog;
667  #endif
668 
669  }
670 
671  } else if (g_dir == KDIR) {
672 
673  r_1 = 0.5*(x1p[i]*x1p[i] - x1p[i-1]*x1p[i-1])/dV1[i];
674  scrh = dt*r_1*dx2[j]/dV2[j];
675  for (k = beg; k <= end; k++) {
676  dtdx = scrh/dx3[k];
677 
678  /* -- K1. initialize rhs with flux difference -- */
679 
680  VAR_LOOP(nv) rhs[k][nv] = -dtdx*(flux[k][nv] - flux[k-1][nv]);
681  rhs[k][iMPHI] -= dtdx*(p[k] - p[k-1]);
682 
683  /* -- K5. Add dissipative terms to entropy equation -- */
684 
685  #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
686  rhog = vh[k][RHO];
687  rhog = (g_gamma - 1.0)*pow(rhog,1.0-g_gamma);
688 
689  rhs_entr = 0.0;
690  #if VISCOSITY == EXPLICIT
691  rhs_entr += (visc_flux[k][ENG] - visc_flux[k-1][ENG]);
692  rhs_entr -= EXPAND( vc[VX1]*(visc_flux[k][MX1] - visc_flux[k-1][MX1]) ,
693  + vc[VX2]*(visc_flux[k][MX2] - visc_flux[k-1][MX2]) ,
694  + vc[VX3]*(visc_flux[k][MX3] - visc_flux[k-1][MX3]));
695  #endif
696  #if THERMAL_CONDUCTION == EXPLICIT
697  rhs_entr += (tc_flux[k][ENG] - tc_flux[k-1][ENG]);
698  #endif
699  rhs[k][ENTR] += rhs_entr*dtdx*rhog;
700  #endif
701  }
702  }
703 }
704 #endif /* GEOMETRY == SPHERICAL */
705 
706 /* --------------------------------------------------------------
707  Add source terms
708  -------------------------------------------------------------- */
709 
710  RightHandSideSource (state, Dts, beg, end, dt, phi_p, grid);
711 
712 /* --------------------------------------------------
713  Reset right hand side in internal boundary zones
714  -------------------------------------------------- */
715 
716  #if INTERNAL_BOUNDARY == YES
717  InternalBoundaryReset(state, Dts, beg, end, grid);
718  #endif
719 
720 /* --------------------------------------------------
721  Time step determination
722  -------------------------------------------------- */
723 
724 #if !GET_MAX_DT
725 return;
726 #endif
727 
728  cl = 0.0;
729  for (i = beg-1; i <= end; i++) {
730  scrh = Dts->cmax[i]*grid[g_dir].inv_dxi[i];
731  cl = MAX(cl, scrh);
732  }
733  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
734  if (g_dir == JDIR) {
735  cl /= fabs(grid[IDIR].xgc[g_i]);
736  }
737  #if GEOMETRY == SPHERICAL
738  if (g_dir == KDIR){
739  cl /= fabs(grid[IDIR].xgc[g_i])*sin(grid[JDIR].xgc[g_j]);
740  }
741  #endif
742  #endif
743  Dts->inv_dta = MAX(cl, Dts->inv_dta);
744 }
MX3
#define MX3
Definition: mod_defs.h:22
MAX
#define MAX(a, b)
Definition: macros.h:101
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
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
IF_ENERGY
#define IF_ENERGY(a)
Definition: pluto.h:631
STATE_1D::visc_src
double ** visc_src
Viscosity source term.
Definition: structs.h:151
iMTH
#define iMTH
Definition: mod_defs.h:70
RightHandSideSource
void RightHandSideSource(const State_1D *state, Time_Step *Dts, int beg, int end, double dt, double *phi_p, Grid *grid)
Definition: rhs_source.c:68
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
TotalFlux
static void TotalFlux(const State_1D *, double *, int, int, Grid *)
Definition: rhs.c:747
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
STATE_1D::res_flux
double ** res_flux
Resistive flux (current)
Definition: structs.h:153
iMPHI
#define iMPHI
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
GetAreaFlux
void GetAreaFlux(const State_1D *state, double **fA, double **fvA, int beg, int end, Grid *grid)
Definition: get_area_flux.c:4
NDUST_LOOP
#define NDUST_LOOP(n)
Definition: pluto.h:617
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
STATE_1D::visc_flux
double ** visc_flux
Viscosity flux.
Definition: structs.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
STATE_1D::vp
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
Definition: structs.h:142
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
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
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
PHYSICS
#define PHYSICS
Definition: definitions_01.h:1
j
int j
Definition: analysis.c:2
MX2
#define MX2
Definition: mod_defs.h:21
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
STATE_1D::tc_flux
double ** tc_flux
Thermal conduction flux.
Definition: structs.h:152
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
GRID::r_1
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
USE_PR_GRADIENT
#define USE_PR_GRADIENT
Definition: rhs.c:83
STATE_1D::vm
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
iBTH
#define iBTH
Definition: mod_defs.h:151
MHD
#define MHD
Definition: pluto.h:111
s
#define s
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
i
int i
Definition: analysis.c:2
InternalBoundaryReset
void InternalBoundaryReset(const State_1D *, Time_Step *, int, int, Grid *)
VX3
#define VX3
Definition: mod_defs.h:30
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
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
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
IF_DUST
#define IF_DUST(a)
Definition: pluto.h:625
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the caller graph for this function:

void TotalFlux ( const State_1D *  state,
double *  phi_p,
int  beg,
int  end,
Grid *  grid 
)
static

Compute the total flux in order to enforce conservation of angular momentum and energy in presence of FARGO source terms, rotation or gravitational potential.

In Cartesian coordinates, for instance, we modify the energy and the y-momentum flux as:

\begin{eqnarray*} \F^{E}_{i+\HALF} &\quad \to\quad& \F^{E}_{i+\HALF} + w_{i+\HALF}\left(\HALF w_{i+\HALF}\F^{\rho}_{i+\HALF} + \F^{m_y}_{i+\HALF}\right) + \Phi_{i+\HALF}\F^{\rho}_{i+\HALF} \\ \noalign{\medskip} \F^{m_y}_{i+\HALF} &\quad \to\quad& \F^{m_y}_{i+\HALF} + w_{i+\HALF}\F^{\rho}_{i+\HALF} \end{eqnarray*}

where $w$ is the average orbital velocity computed during the FARGO algorithm and $\Phi$ is the gravitational potential. This is described in Mignone et al. (A&A 2012), see appendix there.

Note
This function does not change the fluxes when the Shearing-box and FARGO modules are both enabled. Instead, we incorporate the source terms elsewhere.
Parameters
[in]statepointer to State_1D structure;
[in,out]phi_p1D array defining the gravitational potential at grid interfaces;
[in]beginitial index of computation;
[in]endfinal index of computation;
[in]gridpointer to Grid structure;

Reference

  • "A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code" Mignone et al., A&A (2012) 545A, 152M

Definition at line 747 of file rhs.c.

789 {
790  int i,j,k;
791  double wp, R;
792  double **flux, *vp;
793  double *x1, *x2, *x3;
794  double *x1p, *x2p, *x3p;
795 #ifdef FARGO
796  double **wA;
797  wA = FARGO_GetVelocity();
798 #endif
799 
800  flux = state->flux;
801  x1 = grid[IDIR].x; x1p = grid[IDIR].xr;
802  x2 = grid[JDIR].x; x2p = grid[JDIR].xr;
803  x3 = grid[KDIR].x; x3p = grid[KDIR].xr;
804 
805  i = g_i; /* will be redefined during x1-sweep */
806  j = g_j; /* will be redefined during x2-sweep */
807  k = g_k; /* will be redefined during x3-sweep */
808 
809  if (g_dir == IDIR){
810  for (i = beg; i <= end; i++){
811 
812  /* ----------------------------------------------------
813  include flux contributions from FARGO or Rotation
814  Note: ShearingBox terms are not included here but
815  only within the BodyForce function.
816  ---------------------------------------------------- */
817 
818  #if (defined FARGO && !defined SHEARINGBOX) || (ROTATING_FRAME == YES)
819  wp = 0.0;
820  #if GEOMETRY == SPHERICAL
821  IF_FARGO(wp = 0.5*(wA[j][i] + wA[j][i+1]);)
822  R = x1p[i]*sin(x2[j]); /* -- cylindrical radius -- */
823  #else
824  IF_FARGO(wp = 0.5*(wA[k][i] + wA[k][i+1]);)
825  R = x1p[i]; /* -- cylindrical radius -- */
826  #endif
827  IF_ROTATING_FRAME(wp += g_OmegaZ*R;)
828  IF_ENERGY (flux[i][ENG] += wp*(0.5*wp*flux[i][RHO] + flux[i][iMPHI]);)
829  flux[i][iMPHI] += wp*flux[i][RHO];
830  #endif
831 
832  /* -- gravitational potential -- */
833 
834  #if (BODY_FORCE & POTENTIAL) && HAVE_ENERGY
835  flux[i][ENG] += flux[i][RHO]*phi_p[i];
836  #endif
837  }
838  }else if (g_dir == JDIR){
839  for (j = beg; j <= end; j++){
840 
841  /* ----------------------------------------------------
842  include flux contributions from FARGO and Rotation
843  ---------------------------------------------------- */
844 
845  #if GEOMETRY == SPHERICAL
846  #if (defined FARGO) || (ROTATING_FRAME == YES)
847  wp = 0.0;
848  R = x1[i]*sin(x2p[j]);
849  IF_FARGO (wp += 0.5*(wA[j][i] + wA[j+1][i]);)
850  IF_ROTATING_FRAME(wp += g_OmegaZ*R;)
851  IF_ENERGY (flux[j][ENG] += wp*(0.5*wp*flux[j][RHO] + flux[j][iMPHI]);)
852  flux[j][iMPHI] += wp*flux[j][RHO];
853  #endif
854  #endif
855 
856  #if (BODY_FORCE & POTENTIAL) && HAVE_ENERGY
857  flux[j][ENG] += flux[j][RHO]*phi_p[j];
858  #endif
859 
860  }
861  }else if (g_dir == KDIR){
862  for (k = beg; k <= end; k++) {
863 
864  /* ----------------------------------------------------
865  include flux contributions from FARGO
866  (polar/cartesian geometries only)
867  ---------------------------------------------------- */
868 
869  #if (GEOMETRY != SPHERICAL) && (defined FARGO && !defined SHEARINGBOX)
870  wp = 0.5*(wA[k][i] + wA[k+1][i]);
871  IF_ENERGY(flux[k][ENG] += wp*(0.5*wp*flux[k][RHO] + flux[k][iMPHI]);)
872  flux[k][iMPHI] += wp*flux[k][RHO];
873  #endif
874 
875  #if (BODY_FORCE & POTENTIAL) && HAVE_ENERGY
876  flux[k][ENG] += flux[k][RHO]*phi_p[k];
877  #endif
878  }
879  }
880 }
IF_ROTATING_FRAME
#define IF_ROTATING_FRAME(a)
Definition: pluto.h:643
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
IF_ENERGY
#define IF_ENERGY(a)
Definition: pluto.h:631
GRID::xr
double * xr
Definition: structs.h:81
RHO
#define RHO
Definition: mod_defs.h:19
wA
static double ** wA
Definition: fargo_velocity.c:24
g_OmegaZ
#define g_OmegaZ
Definition: init.c:64
iMPHI
#define iMPHI
KDIR
#define KDIR
Definition: pluto.h:195
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
IF_FARGO
#define IF_FARGO(a)
Definition: pluto.h:637
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
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
i
int i
Definition: analysis.c:2
FARGO_GetVelocity
double ** FARGO_GetVelocity(void)
Definition: fargo_velocity.c:270
JDIR
#define JDIR
Definition: pluto.h:194

Here is the call graph for this function:

Here is the caller graph for this function: