PLUTO
Functions
prim_eqn.c File Reference

Compute the right hand side of the relativistic hydro (RHD) 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 *w, double *dw, double cs2, double h, double *Adw)
 
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 relativistic hydro (RHD) equations in primitive form.

Implements the right hand side of the quasi-linear form of the relativistic 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:

  • "The Piecewise Parabolic Method for Multidimensional Relativistic Fluid Dynamics", Mignone, Plewa and Bodo, ApJS (2005) 160,199.

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
July 03, 2015

Definition in file prim_eqn.c.

Function Documentation

void PrimRHS ( double *  w,
double *  dw,
double  cs2,
double  h,
double *  Adw 
)

Compute the matrix-vector multiplication $ A(\mathbf{v})\cdot \Delta\mathbf{v} $ where A is the matrix of the quasi-linear form of the RHD equations.

Parameters
[in]wvector of primitive variables;
[in]dwlimited (linear) slopes;
[in]cs2local sound speed;
[in]hlocal enthalpy;
[out]Adwmatrix-vector product.
Note
Returns
This function has no return value.
Note
In the 7-wave and 8-wave formulations we use the the same matrix being decomposed into right and left eigenvectors during the Characteristic Tracing step. Note, however, that it DOES NOT include two additional terms (-vy*dV[BX] for By, -vz*dv[BX] for Bz) that are needed in the

7-wave form and are added using source terms.

Definition at line 30 of file prim_eqn.c.

45 {
46  int nv;
47  double rho, u1, u2, u3;
48  double d_2, g2, scrh;
49  double g, v1, v2, v3, gx2;
50 
51 #if RECONSTRUCT_4VEL
52  rho = w[RHO];
53  EXPAND(u1 = w[VXn]; ,
54  u2 = w[VXt]; ,
55  u3 = w[VXb];)
56 
57  scrh = EXPAND(u1*u1, + u2*u2, + u3*u3);
58  g2 = 1.0 + scrh;
59  g = sqrt(g2);
60  d_2 = g/(g2 - scrh*cs2);
61 
62 /* Get 3vel */
63 
64  EXPAND(v1 = u1/g; ,
65  v2 = u2/g; ,
66  v3 = u3/g;)
67 
68  gx2 = 1.0/(1.0 - v1*v1);
69 
70  scrh = EXPAND(v1*dw[VXn], + v2*dw[VXt], + v3*dw[VXb]);
71 
72  Adw[PRS] = d_2*(rho*h*cs2*(dw[VXn] - v1*scrh)
73  + u1*(1.0 - cs2)*dw[PRS]);
74 
75  Adw[RHO] = v1*dw[RHO] - (v1*dw[PRS] - Adw[PRS])/(h*cs2);
76 
77  scrh = 1.0/(g*rho*h);
78  d_2 = u1*dw[PRS] - g*Adw[PRS];
79 
80  EXPAND(Adw[VXn] = v1*dw[VXn] + scrh*(dw[PRS] + u1*d_2); ,
81  Adw[VXt] = v1*dw[VXt] + scrh*u2*d_2; ,
82  Adw[VXb] = v1*dw[VXb] + scrh*u3*d_2;)
83 
84 #else
85  rho = w[RHO];
86  EXPAND(v1 = w[VXn]; ,
87  v2 = w[VXt]; ,
88  v3 = w[VXb];)
89 
90  g2 = EXPAND(v1*v1, + v2*v2, + v3*v3);
91  d_2 = 1.0/(1.0 - g2*cs2);
92  g2 = 1.0/(1.0 - g2);
93 
94  Adw[PRS] = d_2*(cs2*rho*h*dw[VXn]
95  + v1*(1.0 - cs2)*dw[PRS]);
96 
97  Adw[RHO] = v1*dw[RHO] - (v1*dw[PRS] - Adw[PRS])/(h*cs2);
98 
99  scrh = 1.0/(g2*rho*h);
100  EXPAND(Adw[VXn] = v1*dw[VXn]
101  + scrh*(dw[PRS] - v1*Adw[PRS]); ,
102  Adw[VXt] = v1*dw[VXt] - scrh*v2*Adw[PRS]; ,
103  Adw[VXb] = v1*dw[VXb] - scrh*v3*Adw[PRS];)
104 #endif
105 
106 #if NSCL > 0
107  NSCL_LOOP(nv) Adw[nv] = v1*dw[nv];
108 #endif
109 
110 }
NSCL_LOOP
#define NSCL_LOOP(n)
Definition: pluto.h:616
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
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 RHD equations in primitive variables.

  • Geometrical sources;
  • Gravity;

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.

Definition at line 113 of file prim_eqn.c.

137 {
138  int nv, i;
139  double r_1, scrh, alpha;
140  double vel2, delta;
141  double *q, *x1, *x2, *x3, g[3];
142 
143  #if GEOMETRY == CYLINDRICAL
144  x1 = grid[IDIR].xgc;
145  x2 = grid[JDIR].xgc;
146  x3 = grid[KDIR].xgc;
147  #else
148  x1 = grid[IDIR].x;
149  x2 = grid[JDIR].x;
150  x3 = grid[KDIR].x;
151  #endif
152 
153  for (i = beg; i <= end; i++){
154  for (nv = NVAR; nv--; ){
155  src[i][nv] = 0.0;
156  }}
157 
158 #if GEOMETRY == CARTESIAN
159 
160 #elif GEOMETRY == CYLINDRICAL
161 
162  if (g_dir == IDIR){
163  for (i = beg; i <= end; i++) {
164 
165  r_1 = 1.0/x1[i];
166  q = state->v[i];
167  vel2 = EXPAND(q[VX1]*q[VX1], +q[VX2]*q[VX2], +q[VX3]*q[VX3]);
168 
169  #if RECONSTRUCT_4VEL
170  scrh = sqrt(1.0 + vel2);
171  alpha = q[VXn]*r_1*scrh/(1.0 + vel2*(1.0 - a2[i]));
172  scrh = a2[i]*alpha;
173  print1 ("! Primitive source terms not yet implemented\n");
174  print1 ("! with 4-vel. Please try 3-vel\n");
175  QUIT_PLUTO(1);
176  #else
177  alpha = q[VXn]*r_1/(1.0 - a2[i]*vel2);
178  scrh = a2[i]*(1.0 - vel2)*alpha;
179  #endif
180 
181  src[i][RHO] = -q[RHO]*alpha;
182  EXPAND (src[i][VX1] = scrh*q[VX1]; ,
183  src[i][VX2] = scrh*q[VX2]; ,
184  src[i][VX3] = scrh*q[VX3];)
185 
186  #if COMPONENTS == 3
187  EXPAND(src[i][iVR] += q[iVPHI]*q[iVPHI]*r_1; ,
188  ,
189  src[i][iVPHI] += -q[iVPHI]*q[iVR]*r_1;)
190  #endif
191 
192  src[i][PRS] = -a2[i]*q[RHO]*h[i]*alpha;
193 
194  }
195  }
196 
197 #elif GEOMETRY == SPHERICAL && DIMENSIONS == 1 && COMPONENTS == 1
198 
199  if (g_dir == IDIR){
200  for (i = beg; i <= end; i++) {
201 
202  r_1 = 1.0/x1[i];
203  q = state->v[i];
204  vel2 = EXPAND(q[VX1]*q[VX1], +q[VX2]*q[VX2], +q[VX3]*u[VX3]);
205 
206  #if RECONSTRUCT_4VEL
207  print1 ("! PrimRHSD(): primitive source terms not yet implemented\n");
208  print1 ("! with 4-vel. Please try 3-vel\n");
209  QUIT_PLUTO(1);
210  #else
211  delta = 1.0 - vel2*a2[i];
212  scrh = 2.0*q[iVR]/delta*r_1;
213  #endif
214 
215  src[i][RHO] = -q[RHO]*scrh;
216  EXPAND (src[i][VX1] = scrh*q[VX1]*a2[i]*(1.0 - vel2); ,
217  src[i][VX2] = scrh*q[VX2]*a2[i]*(1.0 - vel2); ,
218  src[i][VX3] = scrh*q[VX3]*a2[i]*(1.0 - vel2);)
219 
220  src[i][PRS] = src[i][RHO]*h[i]*a2[i];
221  }
222  }
223 
224 #else
225 
226  print1 ("! PrimRHS(): primitive source terms not available for this geometry\n");
227  print1 ("! Use RK integrators\n");
228  QUIT_PLUTO(1);
229 
230 #endif
231 
232 /* -----------------------------------------------------------
233  Add body force
234  ----------------------------------------------------------- */
235 
236  #if (BODY_FORCE != NO)
237  i = beg - 1;
238  if (g_dir == IDIR) {
239  #if BODY_FORCE & POTENTIAL
240  gPhi[i] = BodyForcePotential(x1p[i], x2[g_j], x3[g_k]);
241  #endif
242  for (i = beg; i <= end; i++){
243  #if BODY_FORCE & VECTOR
244  v = state->v[i];
245  BodyForceVector(v, g, x1[i], x2[g_j], x3[g_k]);
246  src[i][VX1] += g[g_dir];
247  #endif
248  #if BODY_FORCE & POTENTIAL
249  gPhi[i] = BodyForcePotential(x1p[i], x2[g_j], x3[g_k]);
250  src[i][VX1] -= (gPhi[i] - gPhi[i-1])/(hscale*dx1[i]);
251  #endif
252 
253  #if DIMENSIONS == 1
254  EXPAND( ,
255  src[i][VX2] += g[JDIR]; ,
256  src[i][VX3] += g[KDIR];)
257  #endif
258  }
259  }else if (g_dir == JDIR){
260  #if BODY_FORCE & POTENTIAL
261  gPhi[i] = BodyForcePotential(x1[g_i], x2p[i], x3[g_k]);
262  #endif
263  for (i = beg; i <= end; i++){
264  #if BODY_FORCE & VECTOR
265  v = state->v[i];
266  BodyForceVector(v, g, x1[g_i], x2[i], x3[g_k]);
267  src[i][VX2] += g[g_dir];
268  #endif
269  #if BODY_FORCE & POTENTIAL
270  gPhi[i] = BodyForcePotential(x1[g_i], x2p[i], x3[g_k]);
271  src[i][VX2] -= (gPhi[i] - gPhi[i-1])/(hscale*dx2[i]);
272  #endif
273 
274  #if DIMENSIONS == 2 && COMPONENTS == 3
275  src[i][VX3] += g[KDIR];
276  #endif
277  }
278  }else if (g_dir == KDIR){
279  #if BODY_FORCE & POTENTIAL
280  gPhi[i] = BodyForcePotential(x1[g_i], x2[g_j], x3p[i]);
281  #endif
282  for (i = beg; i <= end; i++){
283  #if BODY_FORCE & VECTOR
284  v = state->v[i];
285  BodyForceVector(v, g, x1[g_i], x2[g_j], x3[i]);
286  src[i][VX3] += g[g_dir];
287  #endif
288  #if BODY_FORCE & POTENTIAL
289  gPhi[i] = BodyForcePotential(x1[g_i], x2[g_j], x3p[i]);
290  src[i][VX3] -= (gPhi[i] - gPhi[i-1])/(hscale*dx3[i]);
291  #endif
292  }
293  }
294  #endif
295 }
alpha
static double alpha
Definition: init.c:3
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
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
iVPHI
#define iVPHI
Definition: mod_defs.h:67
func_param::v
double v[NVAR]
Definition: eos.h:106
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
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
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
COMPONENTS
#define COMPONENTS
Definition: definitions_01.h:3
i
int i
Definition: analysis.c:2
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
q
static Runtime q
Definition: runtime_setup.c:490
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function: