PLUTO
Functions
prim_eqn.c File Reference

Compute the right hand side of the HD equations in primitive form. 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 HD equations in primitive form.

Implements the right hand side of the quasi-linear form of the hydro 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
April 02, 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 \Delta\mathbf{v} $ where A is the matrix of the quasi-linear form of the HD equations.

Parameters
[in]wvector of primitive variables;
[in]dwlimited (linear) slopes;
[in]cs2local sound speed;
[in]hlocal enthalpy;
[out]Adwmatrix-vector product.
Returns
This function has no return value.

Definition at line 31 of file prim_eqn.c.

44 {
45  int nv;
46  double u;
47 
48  u = v[VXn];
49 
50  Adv[RHO] = u*dv[RHO] + v[RHO]*dv[VXn];
51  #if EOS == IDEAL || EOS == PVTE_LAW
52  EXPAND(Adv[VXn] = u*dv[VXn] + dv[PRS]/v[RHO]; ,
53  Adv[VXt] = u*dv[VXt]; ,
54  Adv[VXb] = u*dv[VXb];)
55  Adv[PRS] = cs2*v[RHO]*dv[VXn] + u*dv[PRS];
56  #elif EOS == ISOTHERMAL
57  EXPAND(Adv[VXn] = u*dv[VXn] + cs2*dv[RHO]/v[RHO]; ,
58  Adv[VXt] = u*dv[VXt]; ,
59  Adv[VXb] = u*dv[VXb];)
60  #endif
61 
62  /* ---- scalars ---- */
63 
64 #if NSCL > 0
65  NSCL_LOOP(nv) Adv[nv] = u*dv[nv];
66 #endif
67 
68 }
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
RHO
#define RHO
Definition: mod_defs.h:19
func_param::v
double v[NVAR]
Definition: eos.h:106
VXb
int VXb
Definition: globals.h:73
NSCL
#define NSCL
Definition: pluto.h:572
VXt
int VXt
Definition: globals.h:73
VXn
int VXn
Definition: globals.h:73
void PrimSource ( const State_1D *  state,
int  beg,
int  end,
double *  a2,
double *  h,
double **  src,
Grid *  grid 
)

Compute source terms of the HD equations in primitive variables.

  • Geometrical sources;
  • Gravity;
  • 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.
Returns
This function has no return value.
Note
This function does not work in spherical coordinates yet.

Definition at line 71 of file prim_eqn.c.

97 {
98  int nv, i, j, k;
99  double tau, dA_dV, th;
100  double hscale; /* scale factor */
101  double *v, *vp, *A, *dV, r_inv, ct;
102  double *x1, *x2, *x3;
103  double *x1p, *x2p, *x3p;
104  double *dx1, *dx2, *dx3;
105  static double *phi_p;
106  double g[3], scrh;
107 
108 #if ROTATING_FRAME == YES
109  print1 ("! PrimSource(): does not work with rotations\n");
110  QUIT_PLUTO(1);
111 #endif
112 
113 /* ----------------------------------------------------------
114  1. Memory allocation and pointer shortcuts
115  ---------------------------------------------------------- */
116 
117  if (phi_p == NULL) phi_p = ARRAY_1D(NMAX_POINT, double);
118 
119  #if GEOMETRY == CYLINDRICAL
120  x1 = grid[IDIR].xgc; x1p = grid[IDIR].xr; dx1 = grid[IDIR].dx;
121  x2 = grid[JDIR].xgc; x2p = grid[JDIR].xr; dx2 = grid[JDIR].dx;
122  x3 = grid[KDIR].xgc; x3p = grid[KDIR].xr; dx3 = grid[KDIR].dx;
123  #else
124  x1 = grid[IDIR].x; x1p = grid[IDIR].xr; dx1 = grid[IDIR].dx;
125  x2 = grid[JDIR].x; x2p = grid[JDIR].xr; dx2 = grid[JDIR].dx;
126  x3 = grid[KDIR].x; x3p = grid[KDIR].xr; dx3 = grid[KDIR].dx;
127  #endif
128 
129  A = grid[g_dir].A;
130  dV = grid[g_dir].dV;
131  hscale = 1.0;
132 
133  i = g_i; j = g_j; k = g_k;
134 
135 /* ----------------------------------------------------------
136  initialize all elements of src to zero
137  ---------------------------------------------------------- */
138 
139  memset((void *)src[0], '\0',NMAX_POINT*NVAR*sizeof(double));
140 
141 /* ----------------------------------------------------------
142  2. Compute geometrical source terms
143  ---------------------------------------------------------- */
144 
145 #if GEOMETRY == CYLINDRICAL
146 
147  if (g_dir == IDIR) {
148  for (i = beg; i <= end; i++){
149  v = state->v[i];
150  tau = 1.0/v[RHO];
151  dA_dV = 1.0/x1[i];
152 
153  src[i][RHO] = -v[RHO]*v[VXn]*dA_dV;
154  #if COMPONENTS == 3
155  src[i][iVR] = v[iVPHI]*v[iVPHI]*dA_dV;
156  src[i][iVPHI] = -v[iVR]*v[iVPHI]*dA_dV;
157  #endif
158  #if EOS == IDEAL
159  src[i][PRS] = a2[i]*src[i][RHO];
160  #endif
161  }
162  }
163 
164 #elif GEOMETRY == POLAR
165 
166  if (g_dir == IDIR) {
167  for (i = beg; i <= end; i++){
168  v = state->v[i];
169 
170  tau = 1.0/v[RHO];
171  dA_dV = 1.0/x1[i];
172  src[i][RHO] = -v[RHO]*v[VXn]*dA_dV;
173 
174  #if COMPONENTS >= 2
175  src[i][iVR] = v[iVPHI]*v[iVPHI]*dA_dV;
176  #endif
177 
178  /* -- in 1D, all sources should be included during this sweep -- */
179 
180  #if DIMENSIONS == 1 && COMPONENTS >= 2
181  src[i][iVPHI] = -v[iVR]*v[iVPHI]*dA_dV;
182  #endif
183 
184  #if EOS == IDEAL
185  src[i][PRS] = a2[i]*src[i][RHO];
186  #endif
187  }
188  } else if (g_dir == JDIR) {
189  dA_dV = 1.0/x1[i];
190  hscale = x1[i];
191  for (j = beg; j <= end; j++){
192  v = state->v[j];
193  tau = 1.0/v[RHO];
194  src[j][iVPHI] = -v[iVR]*v[iVPHI]*dA_dV;
195  }
196  }
197 
198 #elif GEOMETRY == SPHERICAL
199 
200  print1 ("! PrimSource: not implemented in Spherical geometry\n");
201  QUIT_PLUTO(1);
202 
203 #endif
204 
205 /* ----------------------------------------------------------
206  3. Add body forces. This includes:
207  - Coriolis terms for the shearing box module
208  - Body forces
209  ---------------------------------------------------------- */
210 
211 #ifdef SHEARINGBOX
212 
213  if (g_dir == IDIR){
214  for (i = beg; i <= end; i++) {
215  src[i][VX1] = 2.0*state->v[i][VX2]*SB_OMEGA;
216  }
217  }else if (g_dir == JDIR){
218  for (j = beg; j <= end; j++) {
219  #ifdef FARGO
220  src[j][VX2] = (SB_Q - 2.0)*state->v[j][VX1]*SB_OMEGA;
221  #else
222  src[j][VX2] = -2.0*state->v[j][VX1]*SB_OMEGA;
223  #endif
224  }
225  }
226 
227 #endif
228 
229  #if (BODY_FORCE != NO)
230  if (g_dir == IDIR) {
231  i = beg-1;
232  #if BODY_FORCE & POTENTIAL
233  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
234  #endif
235  for (i = beg; i <= end; i++){
236  #if BODY_FORCE & VECTOR
237  v = state->v[i];
238  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
239  src[i][VX1] += g[IDIR];
240  #endif
241  #if BODY_FORCE & POTENTIAL
242  phi_p[i] = BodyForcePotential(x1p[i], x2[j], x3[k]);
243  src[i][VX1] -= (phi_p[i] - phi_p[i-1])/(hscale*dx1[i]);
244  #endif
245 
246  /* -- Add tangential components in 1D -- */
247 
248  #if DIMENSIONS == 1
249  EXPAND( ,
250  src[i][VX2] += g[JDIR]; ,
251  src[i][VX3] += g[KDIR];)
252  #endif
253  }
254  }else if (g_dir == JDIR){
255  j = beg - 1;
256  #if BODY_FORCE & POTENTIAL
257  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
258  #endif
259  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
260  hscale = x1[i];
261  #endif
262  for (j = beg; j <= end; j++){
263  #if BODY_FORCE & VECTOR
264  v = state->v[j];
265  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
266  src[j][VX2] += g[JDIR];
267  #endif
268  #if BODY_FORCE & POTENTIAL
269  phi_p[j] = BodyForcePotential(x1[i], x2p[j], x3[k]);
270  src[j][VX2] -= (phi_p[j] - phi_p[j-1])/(hscale*dx2[j]);
271  #endif
272 
273  /* -- Add 3rd component in 2D -- */
274 
275  #if DIMENSIONS == 2 && COMPONENTS == 3
276  src[j][VX3] += g[KDIR];
277  #endif
278  }
279  }else if (g_dir == KDIR){
280  k = beg - 1;
281  #if BODY_FORCE & POTENTIAL
282  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
283  #endif
284  #if GEOMETRY == SPHERICAL
285  th = x2[j];
286  hscale = x1[i]*sin(th);
287  #endif
288  for (k = beg; k <= end; k++){
289  #if BODY_FORCE & VECTOR
290  v = state->v[k];
291  BodyForceVector(v, g, x1[i], x2[j], x3[k]);
292  src[k][VX3] += g[KDIR];
293  #endif
294  #if BODY_FORCE & POTENTIAL
295  phi_p[k] = BodyForcePotential(x1[i], x2[j], x3p[k]);
296  src[k][VX3] -= (phi_p[k] - phi_p[k-1])/(hscale*dx3[k]);
297  #endif
298  }
299  }
300  #endif
301 
302 /* ---------------------------------------------------------------
303  5. Add FARGO source terms (except for SHEARINGBOX).
304  When SHEARINGBOX is also enabled, we do not include
305  these source terms since they're all provided by body_force)
306  --------------------------------------------------------------- */
307 
308  #if (defined FARGO) && !(defined SHEARINGBOX)
309  #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
310  print1 ("! Time Stepping works only in Cartesian or cylindrical coords\n");
311  print1 ("! Use RK instead\n");
312  QUIT_PLUTO(1);
313  #endif
314 
315  double **wA, *dx, *dz;
316  wA = FARGO_GetVelocity();
317  if (g_dir == IDIR){
318  dx = grid[IDIR].dx;
319  for (i = beg; i <= end; i++){
320  v = state->v[i];
321  src[i][VX2] -= 0.5*v[VX1]*(wA[k][i+1] - wA[k][i-1])/dx[i];
322  }
323  }else if (g_dir == KDIR){
324  dz = grid[KDIR].dx;
325  for (k = beg; k <= end; k++){
326  v = state->v[k];
327  src[k][VX2] -= 0.5*v[VX3]*(wA[k+1][i] - wA[k-1][i])/dz[k];
328  }
329  }
330  #endif
331 }
STATE_1D::v
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
Definition: structs.h:134
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
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
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
VX1
#define VX1
Definition: mod_defs.h:28
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
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
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
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87