PLUTO
rhs_source.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Add source terms to the right hand side of HD/MHD eqns.
5 
6  Add source terms to the right hand side of hydro/MHD equations.
7  These include
8 
9  -# Shearing box source terms;
10  -# Geometrical source term in curvilinear coordinates;
11  -# FARGO or ROTATION source terms (included to enforce conservation of
12  total energy or angular mometum);
13  -# Body forces;
14  -# Powell's 8-waves source terms;
15  -# Extended GLM source terms.
16 
17  \note
18  <B> Shearing-Box source terms: </B>
19  With the shearing-box module, the user is required to provide only
20  gravity while Coriolis source terms are automatically added
21  here.
22  Without FARGO, the source terms of the momentum and energy equations are
23  \f{eqnarray*}{
24  \vec{S}_{\vec{m}} &=&
25  \Big[ 2\Omega_0^2 q x \rho + 2\Omega_0\rho v_y\Big]\hvec{i}
26  + \Big[-2\Omega_0\rho v_x\Big]\hvec{j}
27  + \Big[-\rho\Omega_0^2z\Big]\hvec{k} \\
28  S_E &=& \rho \vec{v} \cdot\vec{g}
29  = \rho v_x (2\Omega_0^2 qx) + \rho v_z (-\Omega_0^2 z)
30  \f}
31  When the FARGO module is used, they modify to
32  \f{eqnarray*}{
33  \vec{S}_{\vec{m}} &=&
34  \Big[ 2\Omega_0\rho v'_y\Big]\hvec{i}
35  + \Big[(q-2)\Omega_0\rho v_x\Big]\hvec{j}
36  + \Big[-\rho\Omega_0^2z\Big]\hvec{k} \\ \noalign{\medskip}
37  S_E &=& (\rho v'_yv_x - B_yB_x)q\Omega_0 + \rho v_z(-\Omega_0^2 z)
38  \f}
39  The energy source term is included during the Fargo transport step,
40  see FARGO_Source().
41 
42  \b Reference
43  - "A conservative orbital advection scheme for simulations
44  of magnetized shear flows with the PLUTO code"\n
45  Mignone et al, A&A (2012) 545, A152 (--> Appendix A.1.1)
46 
47  \author A. Mignone (mignone@ph.unito.it)
48  \date Aug 26, 2015
49 */
50 /* ///////////////////////////////////////////////////////////////////// */
51 #include "pluto.h"
52 
53 #if PHYSICS == MHD
54  #if BACKGROUND_FIELD == YES
55  #define TotBB(v, b0, a, b) (v[a]*(v[b] + b0[b]) + b0[a]*v[b])
56  #else
57  #define TotBB(v, b0, a, b) (v[a]*v[b])
58  #endif
59 #else
60  #define TotBB(v, b0, a, b) 0.0
61 #endif
62 
63 #ifndef iMPHI
64  #define iMPHI MX2 /* -- for Cartesian coordinates -- */
65 #endif
66 
67 /* *********************************************************************** */
68 void RightHandSideSource (const State_1D *state, Time_Step *Dts,
69  int beg, int end, double dt, double *phi_p, Grid *grid)
70 /*!
71  *
72  * \param [in,out] state pointer to State_1D structure
73  * \param [in] Dts pointer to time step structure
74  * \param [in] beg initial index of computation
75  * \param [in] end final index of computation
76  * \param [in] dt time increment
77  * \param [in] phi_p force potential at interfaces
78  * \param [in] grid pointer to Grid structure
79  *
80  * \return This function has no return value.
81  ************************************************************************* */
82 {
83  int i, j, k, nv;
84  double r_1, g[3];
85 
86  double dtdx, scrh, th, ct, s;
87  double Sm;
88  double *x1 = grid[IDIR].x, *x2 = grid[JDIR].x, *x3 = grid[KDIR].x;
89  double *x1p = grid[IDIR].xr, *x2p = grid[JDIR].xr, *x3p = grid[KDIR].xr;
90  double *dx1 = grid[IDIR].dx, *dx2 = grid[JDIR].dx, *dx3 = grid[KDIR].dx;
91  double *dV1 = grid[IDIR].dV, *dV2 = grid[JDIR].dV, *dV3 = grid[KDIR].dV;
92 
93  double **rhs = state->rhs;
94  double **flux = state->flux;
95  double **vh = state->vh;
96  double **vp = state->vp;
97  double **vm = state->vm;
98  double *p;
99  double *A, *dV;
100  double cl;
101  double **Bg0, **wA, w, wp, vphi, phi_c;
102  double vphi_d, Sm_d;
103  double vc[NVAR], *vg;
104 
105 #ifdef FARGO
106  wA = FARGO_GetVelocity();
107 #endif
108 
109 #if (PHYSICS == MHD) && (BACKGROUND_FIELD == YES)
110  Bg0 = GetBackgroundField (beg, end, CELL_CENTER, grid);
111 #endif
112 
113 /* --------------------------
114  pointer shortcuts
115  -------------------------- */
116 
117  p = state->press;
118  A = grid[g_dir].A;
119  dV = grid[g_dir].dV;
120 
121  i = g_i; /* will be redefined during x1-sweep */
122  j = g_j; /* will be redefined during x2-sweep */
123  k = g_k; /* will be redefined during x3-sweep */
124 
125  if (g_dir == IDIR){
126 
127 #if GEOMETRY == SPHERICAL
128  th = x2[j];
129  s = sin(th);
130 #endif
131 
132  for (i = beg; i <= end; i++) {
133  dtdx = dt/dx1[i];
134 
135  /* --------------------------------------------
136  I1. Add geometrical source term
137  -------------------------------------------- */
138 
139 #if GEOMETRY == CARTESIAN
140 
141  vg = state->vh[i];
142  #ifdef SHEARINGBOX /* Include Coriolis term for x-momentum */
143  rhs[i][MX1] += dt*2.0*vh[i][RHO]*vh[i][VX2]*SB_OMEGA;
144  #endif
145 
146  #if (defined FARGO && !defined SHEARINGBOX)
147  w = wA[k][i];
148  #endif
149 
150 #elif GEOMETRY == CYLINDRICAL
151 
152  r_1 = 1.0/x1[i];
153  vg = vc;
154  NFLX_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
155  #if COMPONENTS == 3
156  vphi = vc[iVPHI];
157  IF_ROTATING_FRAME(w = g_OmegaZ*x1[i];
158  vphi += w;)
159  rhs[i][MX1] += dt*(vc[RHO]*vphi*vphi - TotBB(vc, Bg0[i], iBPHI, iBPHI))*r_1;
160  #endif
161 
162 #elif GEOMETRY == POLAR
163 
164  r_1 = grid[IDIR].r_1[i];
165  vg = vc;
166  NFLX_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
167  vphi = vc[iVPHI];
168  #if (defined FARGO) || (ROTATING_FRAME == YES)
169  w = 0.0;
170  IF_FARGO (w += wA[k][i];)
171  IF_ROTATING_FRAME(w += g_OmegaZ*x1[i];)
172  vphi += w;
173  #endif
174  rhs[i][MX1] += dt*( vc[RHO]*vphi*vphi
175  - TotBB(vc, Bg0[i], iBPHI, iBPHI))*r_1;
176 
177  IF_DUST(vphi_d = SELECT(0.0, vc[VX3_D], 0.0);
178  #if (defined FARGO) || (ROTATING_FRAME == YES)
179  vphi_d += w;
180  #endif
181  rhs[i][MX1_D] += dt*vc[RHO_D]*vphi_d*vphi_d*r_1;)
182 
183 #elif GEOMETRY == SPHERICAL
184 
185  r_1 = 1.0/x1[i];
186  vg = vc;
187  NFLX_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
188  vphi = SELECT(0.0, 0.0, vc[iVPHI]);
189  #if (defined FARGO) || (ROTATING_FRAME == YES)
190  w = 0.0;
191  IF_FARGO (w += wA[j][i];)
192  IF_ROTATING_FRAME(w += g_OmegaZ*x1[i]*s;)
193  vphi += w;
194  #endif
195  Sm = vc[RHO]*(EXPAND( 0.0, + vc[VX2]*vc[VX2], + vphi*vphi));
196  Sm += EXPAND( 0.0, - TotBB(vc, Bg0[i], iBTH, iBTH),
197  - TotBB(vc, Bg0[i], iBPHI,iBPHI));
198  rhs[i][MX1] += dt*Sm*r_1;
199 
200  IF_DUST(vphi_d = SELECT(0.0, 0.0, vc[VX3_D]);
201  #if (defined FARGO) || (ROTATING_FRAME == YES)
202  vphi_d += w;
203  #endif
204  Sm_d = vc[RHO_D]*(EXPAND( 0.0, + vc[VX2_D]*vc[VX2_D], + vphi_d*vphi_d));
205  rhs[i][MX1_D] += dt*Sm_d*r_1;)
206 #endif
207 
208  /* ----------------------------------------------------
209  I2. Modify rhs to enforce conservation
210  ---------------------------------------------------- */
211 
212 #if ((defined FARGO && !defined SHEARINGBOX) || (ROTATING_FRAME == YES))
213  rhs[i][iMPHI] -= w*rhs[i][RHO];
214  IF_DUST (rhs[i][iMPHI_D] -= w*rhs[i][RHO_D];)
215  IF_ENERGY(rhs[i][ENG] -= w*(rhs[i][iMPHI] + 0.5*w*rhs[i][RHO]);)
216 #endif
217 
218  /* ----------------------------------------------------
219  I3. Include body forces
220  ---------------------------------------------------- */
221 
222 #if (BODY_FORCE & VECTOR)
223  BodyForceVector(vg, g, x1[i], x2[j], x3[k]);
224  rhs[i][MX1] += dt*vg[RHO]*g[IDIR];
225  IF_DUST (rhs[i][MX1_D] += dt*vg[RHO_D]*g[IDIR];)
226  IF_ENERGY(rhs[i][ENG] += dt*0.5*(flux[i][RHO] + flux[i-1][RHO])*g[IDIR];)
227 #endif
228 
229 #if (BODY_FORCE & POTENTIAL)
230  rhs[i][MX1] -= dtdx*vg[RHO]*(phi_p[i] - phi_p[i-1]);
231  IF_DUST (rhs[i][MX1_D] -= dtdx*vg[RHO_D]*(phi_p[i] - phi_p[i-1]);)
232  IF_ENERGY(phi_c = BodyForcePotential(x1[i], x2[j], x3[k]);
233  rhs[i][ENG] -= phi_c*rhs[i][RHO];)
234 #endif
235 
236 
237  }
238 
239  } else if (g_dir == JDIR){
240 
241  scrh = 1.0;
242 #if GEOMETRY == POLAR
243  scrh = 1.0/x1[i];
244 #elif GEOMETRY == SPHERICAL
245  scrh = r_1 = 0.5*(x1p[i]*x1p[i] - x1p[i-1]*x1p[i-1])/dV1[i];
246 #endif
247  for (j = beg; j <= end; j++) {
248  dtdx = dt/dx2[j]*scrh;
249 
250  /* --------------------------------------------
251  J1. Add geometrical source terms
252  -------------------------------------------- */
253 
254 #if GEOMETRY != SPHERICAL
255 
256  vg = state->vh[j];
257 
258 #ifdef SHEARINGBOX /* Include Coriolis term for y-momentum */
259  #ifdef FARGO
260  rhs[j][MX2] += dt*(SB_Q-2.0)*vh[j][RHO]*vh[j][VX1]*SB_OMEGA;
261  #else
262  rhs[j][MX2] -= dt*2.0*vh[j][RHO]*vh[j][VX1]*SB_OMEGA;
263  #endif
264 #endif
265 
266 
267 #elif GEOMETRY == SPHERICAL
268 
269  th = x2[j];
270  s = sin(th);
271  ct = grid[JDIR].ct[j];
272 
273  vg = vc;
274  NFLX_LOOP(nv) vc[nv] = state->vh[j][nv];
275  vphi = SELECT(0.0, 0.0, vc[iVPHI]);
276  #if (defined FARGO) || (ROTATING_FRAME == YES)
277  w = 0.0;
278  IF_FARGO (w += wA[j][i];)
279  IF_ROTATING_FRAME(w += g_OmegaZ*x1[i]*s;)
280  vphi += w;
281  #endif
282  Sm = vc[RHO]*(EXPAND( 0.0, - vc[iVTH]*vc[iVR], + ct*vphi*vphi));
283  Sm += EXPAND(0.0, + TotBB(vc, Bg0[j], iBTH, iBR),
284  - ct*TotBB(vc, Bg0[j], iBPHI, iBPHI));
285  rhs[j][MX2] += dt*Sm*r_1;
286 
287  IF_DUST(vphi_d = SELECT(0.0, 0.0, vc[VX3_D]);
288  #if (defined FARGO) || (ROTATING_FRAME == YES)
289  vphi_d += w;
290  #endif
291  Sm_d = vc[RHO_D]*(EXPAND( 0.0, - vc[VX2_D]*vc[VX1_D], + ct*vphi_d*vphi_d));
292  rhs[j][MX2_D] += dt*Sm_d*r_1;)
293 
294  /* ----------------------------------------------------
295  J2. Modify rhs to enforce conservation
296  ---------------------------------------------------- */
297 
298  #if (defined FARGO) || (ROTATING_FRAME == YES)
299  rhs[j][MX3] -= w*rhs[j][RHO];
300  IF_DUST (rhs[j][MX3_D] -= w*rhs[j][RHO_D];)
301  IF_ENERGY(rhs[j][ENG] -= w*(rhs[j][MX3] + 0.5*w*rhs[j][RHO]);)
302  #endif
303 
304 #endif /* GEOMETRY == SPHERICAL */
305 
306  /* ----------------------------------------------------
307  J3. Include Body force
308  ---------------------------------------------------- */
309 
310 #if (BODY_FORCE & VECTOR)
311  BodyForceVector(vg, g, x1[i], x2[j], x3[k]);
312  rhs[j][MX2] += dt*vg[RHO]*g[JDIR];
313  IF_DUST (rhs[j][MX2_D] += dt*vg[RHO_D]*g[JDIR];)
314  IF_ENERGY(rhs[j][ENG] += dt*0.5*(flux[j][RHO] + flux[j-1][RHO])*g[JDIR];)
315 #endif
316 
317 #if (BODY_FORCE & POTENTIAL)
318  rhs[j][MX2] -= dtdx*vg[RHO]*(phi_p[j] - phi_p[j-1]);
319  IF_DUST (rhs[j][MX2_D] -= dtdx*vg[RHO_D]*(phi_p[j] - phi_p[j-1]);)
320  IF_ENERGY(phi_c = BodyForcePotential(x1[i], x2[j], x3[k]);
321  rhs[j][ENG] -= phi_c*rhs[j][RHO];)
322 #endif
323  }
324 
325  }else if (g_dir == KDIR){
326 
327  scrh = dt;
328 #if GEOMETRY == SPHERICAL
329  r_1 = 0.5*(x1p[i]*x1p[i] - x1p[i-1]*x1p[i-1])/dV1[i];
330  scrh = dt*r_1*dx2[j]/dV2[j];
331 #endif
332 
333  for (k = beg; k <= end; k++) {
334  dtdx = scrh/dx3[k];
335  vg = state->vh[k];
336 
337 #if (GEOMETRY == CARTESIAN) || (GEOMETRY == POLAR)
338  /* ------------------------------------------------------
339  K2. modify rhs to enforce conservation (FARGO only)
340  (solid body rotations are not included the
341  velocity depends on the cylindrical radius only)
342  ------------------------------------------------------ */
343 
344  #if (defined FARGO && !defined SHEARINGBOX)
345  w = wA[k][i];
346  rhs[k][MX2] -= w*rhs[k][RHO];
347  IF_DUST (rhs[k][MX2_D] -= w*rhs[k][RHO_D];)
348  IF_ENERGY(rhs[k][ENG] -= w*(rhs[k][MX2] + 0.5*w*rhs[k][RHO]);)
349  #endif
350 #endif
351 
352  /* ----------------------------------------------------
353  K3. Include body forces
354  ---------------------------------------------------- */
355 
356 #if (BODY_FORCE & VECTOR)
357  BodyForceVector(vg, g, x1[i], x2[j], x3[k]);
358  rhs[k][MX3] += dt*vg[RHO]*g[KDIR];
359  IF_DUST(rhs[k][MX3_D] += dt*vg[RHO_D]*g[KDIR];)
360  IF_ENERGY(rhs[k][ENG] += dt*0.5*(flux[k][RHO] + flux[k-1][RHO])*g[KDIR];)
361 #endif
362 
363 #if (BODY_FORCE & POTENTIAL)
364  rhs[k][MX3] -= dtdx*vg[RHO]*(phi_p[k] - phi_p[k-1]);
365  IF_DUST (rhs[k][MX3_D] -= dtdx*vg[RHO_D]*(phi_p[k] - phi_p[k-1]);)
366  IF_ENERGY(phi_c = BodyForcePotential(x1[i], x2[j], x3[k]);
367  rhs[k][ENG] -= phi_c*rhs[k][RHO];)
368 #endif
369  }
370  }
371 
372 /* ---------------------------------------------------------------
373  Source terms coming from tensor discretazion of parabolic
374  terms in curvilinear coordinates (only for viscosity)
375  ---------------------------------------------------------------- */
376 
377 #if (VISCOSITY == EXPLICIT) && (GEOMETRY != CARTESIAN)
378  for (i = beg; i <= end; i++) {
379  EXPAND(rhs[i][MX1] += dt*state->visc_src[i][MX1]; ,
380  rhs[i][MX2] += dt*state->visc_src[i][MX2]; ,
381  rhs[i][MX3] += dt*state->visc_src[i][MX3];)
382  }
383 #endif
384 
385 /* --------------------------------------------------
386  Dust drag force term
387  -------------------------------------------------- */
388 
389 #if DUST == YES
390  Dust_DragForce(state, beg, end, dt, grid);
391 #endif
392 
393 /* --------------------------------------------------
394  Powell's source terms
395  -------------------------------------------------- */
396 
397 #if (PHYSICS == MHD) && (DIVB_CONTROL == EIGHT_WAVES)
398  for (i = beg; i <= end; i++) {
399  EXPAND(rhs[i][MX1] += dt*state->src[i][MX1]; ,
400  rhs[i][MX2] += dt*state->src[i][MX2]; ,
401  rhs[i][MX3] += dt*state->src[i][MX3];)
402 
403  EXPAND(rhs[i][BX1] += dt*state->src[i][BX1]; ,
404  rhs[i][BX2] += dt*state->src[i][BX2]; ,
405  rhs[i][BX3] += dt*state->src[i][BX3];)
406  #if HAVE_ENERGY
407  rhs[i][ENG] += dt*state->src[i][ENG];
408  #endif
409  }
410 #endif
411 
412 /* -------------------------------------------------
413  Extended GLM source terms
414  ------------------------------------------------- */
415 
416 #if (defined GLM_MHD) && (GLM_EXTENDED == YES)
417  GLM_ExtendedSource (state, dt, beg, end, grid);
418 #endif
419 
420 }
421 #undef TotBB
#define MX3
Definition: mod_defs.h:22
#define IF_ROTATING_FRAME(a)
Definition: pluto.h:643
#define iMPHI
Definition: rhs_source.c:64
#define MX1
Definition: mod_defs.h:20
double ** vh
Primitive state at n+1/2 (only for one step method)
Definition: structs.h:162
#define VX2
Definition: mod_defs.h:29
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
#define IF_ENERGY(a)
Definition: pluto.h:631
double ** visc_src
Viscosity source term.
Definition: structs.h:151
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
void GLM_ExtendedSource(const State_1D *state, double dt, int beg, int end, Grid *grid)
Definition: glm.c:168
double ** rhs
Conservative right hand side.
Definition: structs.h:163
double * xr
Definition: structs.h:81
#define RHO
Definition: mod_defs.h:19
tuple scrh
Definition: configure.py:200
#define BODY_FORCE
Definition: definitions_01.h:5
#define YES
Definition: pluto.h:25
static double ** wA
#define iVPHI
Definition: mod_defs.h:67
#define g_OmegaZ
Definition: init.c:64
#define VISCOSITY
Definition: pluto.h:385
#define SPHERICAL
Definition: pluto.h:36
double * dV
Cell volume.
Definition: structs.h:86
double * dx
Definition: structs.h:83
#define iBR
Definition: mod_defs.h:150
double ** GetBackgroundField(int beg, int end, int where, Grid *grid)
Definition: bckgrnd_field.c:5
#define VX1
Definition: mod_defs.h:28
#define POTENTIAL
Definition: pluto.h:129
#define KDIR
Definition: pluto.h:195
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
Definition: structs.h:142
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
#define CELL_CENTER
Definition: pluto.h:205
double ** src
Definition: structs.h:160
#define IF_FARGO(a)
Definition: pluto.h:637
#define iVR
Definition: mod_defs.h:65
#define IDIR
Definition: pluto.h:193
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
Definition: structs.h:78
#define SB_OMEGA
Disk local orbital frequency .
Definition: shearingbox.h:81
int j
Definition: analysis.c:2
#define iVTH
Definition: mod_defs.h:66
#define MX2
Definition: mod_defs.h:21
#define iBPHI
Definition: mod_defs.h:152
#define NFLX_LOOP(n)
Definition: pluto.h:613
int k
Definition: analysis.c:2
double * x
Definition: structs.h:80
double BodyForcePotential(double, double, double)
Definition: init.c:479
#define TotBB(v, b0, a, b)
Definition: rhs_source.c:55
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
#define GEOMETRY
Definition: definitions_01.h:4
#define SB_Q
The shear parameter, .
Definition: shearingbox.h:76
double * r_1
Geometrical factor 1/r.
Definition: structs.h:88
#define CARTESIAN
Definition: pluto.h:33
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
Definition: structs.h:141
#define iBTH
Definition: mod_defs.h:151
#define s
#define BX3
Definition: mod_defs.h:27
#define VECTOR
Definition: pluto.h:128
PLUTO main header file.
int i
Definition: analysis.c:2
double ** FARGO_GetVelocity(void)
void BodyForceVector(double *, double *, double, double, double)
Definition: init.c:441
#define BX1
Definition: mod_defs.h:25
double * ct
Geometrical factor cot(theta).
Definition: structs.h:89
double * press
Upwind pressure term computed with the Riemann solver.
Definition: structs.h:164
#define BX2
Definition: mod_defs.h:26
#define JDIR
Definition: pluto.h:194
#define NVAR
Definition: pluto.h:609
#define EXPLICIT
Definition: pluto.h:67
#define ROTATING_FRAME
Definition: pluto.h:361
#define IF_DUST(a)
Definition: pluto.h:625
double * A
Right interface area, A[i] = .
Definition: structs.h:87