HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
navierstokes3d.h
Go to the documentation of this file.
1 
47 #include <basic.h>
48 
50 #define _NAVIER_STOKES_3D_ "navierstokes3d"
51 
52 /* define ndims and nvars for this model */
53 #undef _MODEL_NDIMS_
54 #undef _MODEL_NVARS_
55 
56 #define _MODEL_NDIMS_ 3
57 
58 #define _MODEL_NVARS_ 5
59 
60 /* choices for upwinding schemes */
62 #define _ROE_ "roe"
63 
64 #define _RF_ "rf-char"
65 
66 #define _LLF_ "llf-char"
67 
68 #define _RUSANOV_ "rusanov"
69 
70 /* grid directions */
72 #define _XDIR_ 0
73 
74 #define _YDIR_ 1
75 
76 #define _ZDIR_ 2
77 
78 /* immersed boundary wall types */
80 #define _IB_ADIABATIC_ "adiabatic"
81 
82 #define _IB_ISOTHERMAL_ "isothermal"
83 
84 /* types of immersed boundary application ramps */
86 #define _IB_RAMP_LINEAR_ "linear"
87 
88 #define _IB_RAMP_SMOOTHEDSLAB_ "smoothed_slab"
89 
90 #define _IB_RAMP_DISABLE_ "no_ib"
91 
98 #define _NavierStokes3DGetFlowVar_(u,stride,rho,vx,vy,vz,e,P,gamma) \
99  { \
100  double vsq; \
101  rho = u[0]; \
102  vx = (rho==0) ? 0 : u[stride] / rho; \
103  vy = (rho==0) ? 0 : u[2*stride] / rho; \
104  vz = (rho==0) ? 0 : u[3*stride] / rho; \
105  e = u[4*stride]; \
106  vsq = vx*vx + vy*vy + vz*vz; \
107  P = (e - 0.5*rho*vsq) * (gamma-1.0); \
108  }
109 
118 #define _NavierStokes3DSetFlux_(f,stride,rho,vx,vy,vz,e,P,dir) \
119  { \
120  if (dir == _XDIR_) { \
121  f[0] = rho * vx; \
122  f[stride] = rho * vx * vx + P; \
123  f[2*stride] = rho * vx * vy; \
124  f[3*stride] = rho * vx * vz; \
125  f[4*stride] = (e + P) * vx; \
126  } else if (dir == _YDIR_) { \
127  f[0] = rho * vy; \
128  f[stride] = rho * vy * vx; \
129  f[2*stride] = rho * vy * vy + P; \
130  f[3*stride] = rho * vy * vz; \
131  f[4*stride] = (e + P) * vy; \
132  } else if (dir == _ZDIR_) { \
133  f[0] = rho * vz; \
134  f[stride] = rho * vz * vx; \
135  f[2*stride] = rho * vz * vy; \
136  f[3*stride] = rho * vz * vz + P; \
137  f[4*stride] = (e + P) * vz; \
138  } \
139  }
140 
144 #define _NavierStokes3DRoeAverage_(uavg,stride,uL,uR,gamma) \
145  { \
146  double rho ,vx, vy, vz, e ,P ,H; \
147  double rhoL,vxL,vyL,vzL,eL,PL,HL,cLsq; \
148  double rhoR,vxR,vyR,vzR,eR,PR,HR,cRsq; \
149  _NavierStokes3DGetFlowVar_(uL,stride,rhoL,vxL,vyL,vzL,eL,PL,gamma); \
150  cLsq = gamma * PL/rhoL; \
151  HL = 0.5*(vxL*vxL+vyL*vyL+vzL*vzL) + cLsq / (gamma-1.0); \
152  _NavierStokes3DGetFlowVar_(uR,stride,rhoR,vxR,vyR,vzR,eR,PR,gamma); \
153  cRsq = gamma * PR/rhoR; \
154  HR = 0.5*(vxR*vxR+vyR*vyR+vzR*vzR) + cRsq / (gamma-1.0); \
155  double tL = sqrt(rhoL); \
156  double tR = sqrt(rhoR); \
157  rho = tL * tR; \
158  vx = (tL*vxL + tR*vxR) / (tL + tR); \
159  vy = (tL*vyL + tR*vyR) / (tL + tR); \
160  vz = (tL*vzL + tR*vzR) / (tL + tR); \
161  H = (tL*HL + tR*HR ) / (tL + tR); \
162  P = (H - 0.5* (vx*vx+vy*vy+vz*vz)) * (rho*(gamma-1.0))/gamma; \
163  e = P/(gamma-1.0) + 0.5*rho*(vx*vx+vy*vy+vz*vz); \
164  uavg[0] = rho; \
165  uavg[1] = rho*vx; \
166  uavg[2] = rho*vy; \
167  uavg[3] = rho*vz; \
168  uavg[4] = e; \
169  }
170 
183 #define _NavierStokes3DSetStiffFlux_(f,stride,rho,vx,vy,vz,e,P,dir,gamma) \
184  { \
185  double gamma_inv = 1.0/gamma; \
186  if (dir == _XDIR_) { \
187  f[0*stride] = gamma_inv * rho * vx; \
188  f[1*stride] = gamma_inv * rho * vx * vx + P; \
189  f[2*stride] = gamma_inv * rho * vx * vy; \
190  f[3*stride] = gamma_inv * rho * vx * vz; \
191  f[4*stride] = (e + P) * vx - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy+vz*vz) * vx; \
192  } else if (dir == _YDIR_) { \
193  f[0*stride] = gamma_inv * rho * vy; \
194  f[1*stride] = gamma_inv * rho * vy * vx; \
195  f[2*stride] = gamma_inv * rho * vy * vy + P; \
196  f[3*stride] = gamma_inv * rho * vy * vz; \
197  f[4*stride] = (e + P) * vy - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy+vz*vz) * vy; \
198  } else if (dir == _ZDIR_) { \
199  f[0*stride] = gamma_inv * rho * vz; \
200  f[1*stride] = gamma_inv * rho * vz * vx; \
201  f[2*stride] = gamma_inv * rho * vz * vy; \
202  f[3*stride] = gamma_inv * rho * vz * vz + P; \
203  f[4*stride] = (e + P) * vz - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy+vz*vz) * vz; \
204  } \
205  }
206 
219 #define _NavierStokes3DSetNonStiffFlux_(f,rho,vx,vy,vz,e,P,dir,gamma) \
220  { \
221  double gamma_inv = 1.0/gamma; \
222  if (dir == _XDIR_) { \
223  f[0*stride] = (gamma-1.0) * gamma_inv * rho * vx; \
224  f[1*stride] = (gamma-1.0) * gamma_inv * rho * vx * vx; \
225  f[2*stride] = (gamma-1.0) * gamma_inv * rho * vx * vy; \
226  f[3*stride] = (gamma-1.0) * gamma_inv * rho * vx * vz; \
227  f[4*stride] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy+vz*vz) * vx; \
228  } else if (dir == _YDIR_) { \
229  f[0*stride] = (gamma-1.0) * gamma_inv * rho * vy; \
230  f[1*stride] = (gamma-1.0) * gamma_inv * rho * vy * vx; \
231  f[2*stride] = (gamma-1.0) * gamma_inv * rho * vy * vy; \
232  f[3*stride] = (gamma-1.0) * gamma_inv * rho * vy * vz; \
233  f[4*stride] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy+vz*vz) * vy; \
234  } else if (dir == _ZDIR_) { \
235  f[0*stride] = (gamma-1.0) * gamma_inv * rho * vz; \
236  f[1*stride] = (gamma-1.0) * gamma_inv * rho * vz * vx; \
237  f[2*stride] = (gamma-1.0) * gamma_inv * rho * vz * vy; \
238  f[3*stride] = (gamma-1.0) * gamma_inv * rho * vz * vz; \
239  f[4*stride] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy+vz*vz) * vz; \
240  } \
241  }
242 
254 #define _NavierStokes3DEigenvalues_(u,stride,D,gamma,dir) \
255  { \
256  double rho,vx,vy,vz,e,P,c; \
257  _NavierStokes3DGetFlowVar_(u,stride,rho,vx,vy,vz,e,P,gamma); \
258  c = sqrt(gamma*P/rho); \
259  if (dir == _XDIR_) { \
260  D[0*_MODEL_NVARS_+0] = vx; D[0*_MODEL_NVARS_+1] = 0; D[0*_MODEL_NVARS_+2] = 0; D[0*_MODEL_NVARS_+3] = 0; D[0*_MODEL_NVARS_+4] = 0; \
261  D[1*_MODEL_NVARS_+0] = 0; D[1*_MODEL_NVARS_+1] = vx-c; D[1*_MODEL_NVARS_+2] = 0; D[1*_MODEL_NVARS_+3] = 0; D[1*_MODEL_NVARS_+4] = 0; \
262  D[2*_MODEL_NVARS_+0] = 0; D[2*_MODEL_NVARS_+1] = 0; D[2*_MODEL_NVARS_+2] = vx; D[2*_MODEL_NVARS_+3] = 0; D[2*_MODEL_NVARS_+4] = 0; \
263  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vx; D[3*_MODEL_NVARS_+4] = 0; \
264  D[4*_MODEL_NVARS_+0] = 0; D[4*_MODEL_NVARS_+1] = 0; D[4*_MODEL_NVARS_+2] = 0; D[4*_MODEL_NVARS_+3] = 0; D[4*_MODEL_NVARS_+4] = vx+c;\
265  } else if (dir == _YDIR_) { \
266  D[0*_MODEL_NVARS_+0] = vy; D[0*_MODEL_NVARS_+1] = 0; D[0*_MODEL_NVARS_+2] = 0; D[0*_MODEL_NVARS_+3] = 0; D[0*_MODEL_NVARS_+4] = 0; \
267  D[1*_MODEL_NVARS_+0] = 0; D[1*_MODEL_NVARS_+1] = vy; D[1*_MODEL_NVARS_+2] = 0; D[1*_MODEL_NVARS_+3] = 0; D[1*_MODEL_NVARS_+4] = 0; \
268  D[2*_MODEL_NVARS_+0] = 0; D[2*_MODEL_NVARS_+1] = 0; D[2*_MODEL_NVARS_+2] = vy-c; D[2*_MODEL_NVARS_+3] = 0; D[2*_MODEL_NVARS_+4] = 0; \
269  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vy; D[3*_MODEL_NVARS_+4] = 0; \
270  D[4*_MODEL_NVARS_+0] = 0; D[4*_MODEL_NVARS_+1] = 0; D[4*_MODEL_NVARS_+2] = 0; D[4*_MODEL_NVARS_+3] = 0; D[4*_MODEL_NVARS_+4] = vy+c;\
271  } else if (dir == _ZDIR_) { \
272  D[0*_MODEL_NVARS_+0] = vz; D[0*_MODEL_NVARS_+1] = 0; D[0*_MODEL_NVARS_+2] = 0; D[0*_MODEL_NVARS_+3] = 0; D[0*_MODEL_NVARS_+4] = 0; \
273  D[1*_MODEL_NVARS_+0] = 0; D[1*_MODEL_NVARS_+1] = vz; D[1*_MODEL_NVARS_+2] = 0; D[1*_MODEL_NVARS_+3] = 0; D[1*_MODEL_NVARS_+4] = 0; \
274  D[2*_MODEL_NVARS_+0] = 0; D[2*_MODEL_NVARS_+1] = 0; D[2*_MODEL_NVARS_+2] = vz; D[2*_MODEL_NVARS_+3] = 0; D[2*_MODEL_NVARS_+4] = 0; \
275  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vz-c; D[3*_MODEL_NVARS_+4] = 0; \
276  D[4*_MODEL_NVARS_+0] = 0; D[4*_MODEL_NVARS_+1] = 0; D[4*_MODEL_NVARS_+2] = 0; D[4*_MODEL_NVARS_+3] = 0; D[4*_MODEL_NVARS_+4] = vz+c;\
277  } \
278  }
279 
288 #define _NavierStokes3DLeftEigenvectors_(u,stride,L,ga,dir) \
289  { \
290  double ga_minus_one=ga-1.0; \
291  double rho,vx,vy,vz,e,P,a,ek; \
292  _NavierStokes3DGetFlowVar_(u,stride,rho,vx,vy,vz,e,P,ga); \
293  ek = 0.5 * (vx*vx + vy*vy + vz*vz); \
294  a = sqrt(ga * P / rho); \
295  if (dir == _XDIR_) { \
296  L[1*_MODEL_NVARS_+0] = (ga_minus_one*ek + a*vx) / (2*a*a); \
297  L[1*_MODEL_NVARS_+1] = ((-ga_minus_one)*vx-a) / (2*a*a); \
298  L[1*_MODEL_NVARS_+2] = ((-ga_minus_one)*vy) / (2*a*a); \
299  L[1*_MODEL_NVARS_+3] = ((-ga_minus_one)*vz) / (2*a*a); \
300  L[1*_MODEL_NVARS_+4] = ga_minus_one / (2*a*a); \
301  L[0*_MODEL_NVARS_+0] = (a*a - ga_minus_one*ek) / (a*a); \
302  L[0*_MODEL_NVARS_+1] = (ga_minus_one*vx) / (a*a); \
303  L[0*_MODEL_NVARS_+2] = (ga_minus_one*vy) / (a*a); \
304  L[0*_MODEL_NVARS_+3] = (ga_minus_one*vz) / (a*a); \
305  L[0*_MODEL_NVARS_+4] = (-ga_minus_one) / (a*a); \
306  L[4*_MODEL_NVARS_+0] = (ga_minus_one*ek - a*vx) / (2*a*a); \
307  L[4*_MODEL_NVARS_+1] = ((-ga_minus_one)*vx+a) / (2*a*a); \
308  L[4*_MODEL_NVARS_+2] = ((-ga_minus_one)*vy) / (2*a*a); \
309  L[4*_MODEL_NVARS_+3] = ((-ga_minus_one)*vz) / (2*a*a); \
310  L[4*_MODEL_NVARS_+4] = ga_minus_one / (2*a*a); \
311  L[2*_MODEL_NVARS_+0] = vy; \
312  L[2*_MODEL_NVARS_+1] = 0.0; \
313  L[2*_MODEL_NVARS_+2] = -1.0; \
314  L[2*_MODEL_NVARS_+3] = 0.0; \
315  L[2*_MODEL_NVARS_+4] = 0.0; \
316  L[3*_MODEL_NVARS_+0] = -vz; \
317  L[3*_MODEL_NVARS_+1] = 0.0; \
318  L[3*_MODEL_NVARS_+2] = 0.0; \
319  L[3*_MODEL_NVARS_+3] = 1.0; \
320  L[3*_MODEL_NVARS_+4] = 0.0; \
321  } else if (dir == _YDIR_) { \
322  L[2*_MODEL_NVARS_+0] = (ga_minus_one*ek+a*vy) / (2*a*a); \
323  L[2*_MODEL_NVARS_+1] = ((1.0-ga)*vx) / (2*a*a); \
324  L[2*_MODEL_NVARS_+2] = ((1.0-ga)*vy-a) / (2*a*a); \
325  L[2*_MODEL_NVARS_+3] = ((1.0-ga)*vz) / (2*a*a); \
326  L[2*_MODEL_NVARS_+4] = ga_minus_one / (2*a*a); \
327  L[0*_MODEL_NVARS_+0] = (a*a-ga_minus_one*ek) / (a*a); \
328  L[0*_MODEL_NVARS_+1] = ga_minus_one*vx / (a*a); \
329  L[0*_MODEL_NVARS_+2] = ga_minus_one*vy / (a*a); \
330  L[0*_MODEL_NVARS_+3] = ga_minus_one*vz / (a*a); \
331  L[0*_MODEL_NVARS_+4] = (1.0 - ga) / (a*a); \
332  L[4*_MODEL_NVARS_+0] = (ga_minus_one*ek-a*vy) / (2*a*a); \
333  L[4*_MODEL_NVARS_+1] = ((1.0-ga)*vx) / (2*a*a); \
334  L[4*_MODEL_NVARS_+2] = ((1.0-ga)*vy+a) / (2*a*a); \
335  L[4*_MODEL_NVARS_+3] = ((1.0-ga)*vz) / (2*a*a); \
336  L[4*_MODEL_NVARS_+4] = ga_minus_one / (2*a*a); \
337  L[1*_MODEL_NVARS_+0] = -vx; \
338  L[1*_MODEL_NVARS_+1] = 1.0; \
339  L[1*_MODEL_NVARS_+2] = 0.0; \
340  L[1*_MODEL_NVARS_+3] = 0.0; \
341  L[1*_MODEL_NVARS_+4] = 0; \
342  L[3*_MODEL_NVARS_+0] = vz; \
343  L[3*_MODEL_NVARS_+1] = 0.0; \
344  L[3*_MODEL_NVARS_+2] = 0.0; \
345  L[3*_MODEL_NVARS_+3] = -1.0; \
346  L[3*_MODEL_NVARS_+4] = 0; \
347  } else if (dir == _ZDIR_) { \
348  L[3*_MODEL_NVARS_+0] = (ga_minus_one*ek+a*vz) / (2*a*a); \
349  L[3*_MODEL_NVARS_+1] = ((1.0-ga)*vx) / (2*a*a); \
350  L[3*_MODEL_NVARS_+2] = ((1.0-ga)*vy) / (2*a*a); \
351  L[3*_MODEL_NVARS_+3] = ((1.0-ga)*vz-a) / (2*a*a); \
352  L[3*_MODEL_NVARS_+4] = ga_minus_one / (2*a*a); \
353  L[0*_MODEL_NVARS_+0] = (a*a-ga_minus_one*ek) / (a*a); \
354  L[0*_MODEL_NVARS_+1] = ga_minus_one*vx / (a*a); \
355  L[0*_MODEL_NVARS_+2] = ga_minus_one*vy / (a*a); \
356  L[0*_MODEL_NVARS_+3] = ga_minus_one*vz / (a*a); \
357  L[0*_MODEL_NVARS_+4] = (1.0-ga) / (a*a); \
358  L[4*_MODEL_NVARS_+0] = (ga_minus_one*ek-a*vz) / (2*a*a); \
359  L[4*_MODEL_NVARS_+1] = ((1.0-ga)*vx) / (2*a*a); \
360  L[4*_MODEL_NVARS_+2] = ((1.0-ga)*vy) / (2*a*a); \
361  L[4*_MODEL_NVARS_+3] = ((1.0-ga)*vz+a) / (2*a*a); \
362  L[4*_MODEL_NVARS_+4] = ga_minus_one / (2*a*a); \
363  L[1*_MODEL_NVARS_+0] = vx; \
364  L[1*_MODEL_NVARS_+1] = -1.0; \
365  L[1*_MODEL_NVARS_+2] = 0.0; \
366  L[1*_MODEL_NVARS_+3] = 0.0; \
367  L[1*_MODEL_NVARS_+4] = 0; \
368  L[2*_MODEL_NVARS_+0] = -vy; \
369  L[2*_MODEL_NVARS_+1] = 0.0; \
370  L[2*_MODEL_NVARS_+2] = 1.0; \
371  L[2*_MODEL_NVARS_+3] = 0.0; \
372  L[2*_MODEL_NVARS_+4] = 0; \
373  } \
374  }
375 
384 #define _NavierStokes3DRightEigenvectors_(u,stride,R,ga,dir) \
385  { \
386  double ga_minus_one = ga-1.0; \
387  double rho,vx,vy,vz,e,P,ek,a,h0; \
388  _NavierStokes3DGetFlowVar_(u,stride,rho,vx,vy,vz,e,P,ga); \
389  ek = 0.5 * (vx*vx + vy*vy + vz*vz); \
390  a = sqrt(ga * P / rho); \
391  h0 = a*a / ga_minus_one + ek; \
392  if (dir == _XDIR_) { \
393  R[0*_MODEL_NVARS_+1] = 1.0; \
394  R[1*_MODEL_NVARS_+1] = vx-a; \
395  R[2*_MODEL_NVARS_+1] = vy; \
396  R[3*_MODEL_NVARS_+1] = vz; \
397  R[4*_MODEL_NVARS_+1] = h0 - a*vx; \
398  R[0*_MODEL_NVARS_+0] = 1.0; \
399  R[1*_MODEL_NVARS_+0] = vx; \
400  R[2*_MODEL_NVARS_+0] = vy; \
401  R[3*_MODEL_NVARS_+0] = vz; \
402  R[4*_MODEL_NVARS_+0] = ek; \
403  R[0*_MODEL_NVARS_+4] = 1.0; \
404  R[1*_MODEL_NVARS_+4] = vx+a; \
405  R[2*_MODEL_NVARS_+4] = vy; \
406  R[3*_MODEL_NVARS_+4] = vz; \
407  R[4*_MODEL_NVARS_+4] = h0 + a*vx; \
408  R[0*_MODEL_NVARS_+2] = 0.0; \
409  R[1*_MODEL_NVARS_+2] = 0.0; \
410  R[2*_MODEL_NVARS_+2] = -1.0; \
411  R[3*_MODEL_NVARS_+2] = 0.0; \
412  R[4*_MODEL_NVARS_+2] = -vy; \
413  R[0*_MODEL_NVARS_+3] = 0.0; \
414  R[1*_MODEL_NVARS_+3] = 0.0; \
415  R[2*_MODEL_NVARS_+3] = 0.0; \
416  R[3*_MODEL_NVARS_+3] = 1.0; \
417  R[4*_MODEL_NVARS_+3] = vz; \
418  } else if (dir == _YDIR_) { \
419  R[0*_MODEL_NVARS_+2] = 1.0; \
420  R[1*_MODEL_NVARS_+2] = vx; \
421  R[2*_MODEL_NVARS_+2] = vy-a; \
422  R[3*_MODEL_NVARS_+2] = vz; \
423  R[4*_MODEL_NVARS_+2] = h0 - a*vy; \
424  R[0*_MODEL_NVARS_+0] = 1.0; \
425  R[1*_MODEL_NVARS_+0] = vx; \
426  R[2*_MODEL_NVARS_+0] = vy; \
427  R[3*_MODEL_NVARS_+0] = vz; \
428  R[4*_MODEL_NVARS_+0] = ek; \
429  R[0*_MODEL_NVARS_+4] = 1.0; \
430  R[1*_MODEL_NVARS_+4] = vx; \
431  R[2*_MODEL_NVARS_+4] = vy+a; \
432  R[3*_MODEL_NVARS_+4] = vz; \
433  R[4*_MODEL_NVARS_+4] = h0 + a*vy; \
434  R[0*_MODEL_NVARS_+1] = 0.0; \
435  R[1*_MODEL_NVARS_+1] = 1.0; \
436  R[2*_MODEL_NVARS_+1] = 0.0; \
437  R[3*_MODEL_NVARS_+1] = 0.0; \
438  R[4*_MODEL_NVARS_+1] = vx; \
439  R[0*_MODEL_NVARS_+3] = 0.0; \
440  R[1*_MODEL_NVARS_+3] = 0.0; \
441  R[2*_MODEL_NVARS_+3] = 0.0; \
442  R[3*_MODEL_NVARS_+3] = -1.0; \
443  R[4*_MODEL_NVARS_+3] = -vz; \
444  } else if (dir == _ZDIR_) { \
445  R[0*_MODEL_NVARS_+3] = 1.0; \
446  R[1*_MODEL_NVARS_+3] = vx; \
447  R[2*_MODEL_NVARS_+3] = vy; \
448  R[3*_MODEL_NVARS_+3] = vz-a; \
449  R[4*_MODEL_NVARS_+3] = h0-a*vz; \
450  R[0*_MODEL_NVARS_+0] = 1.0; \
451  R[1*_MODEL_NVARS_+0] = vx; \
452  R[2*_MODEL_NVARS_+0] = vy; \
453  R[3*_MODEL_NVARS_+0] = vz; \
454  R[4*_MODEL_NVARS_+0] = ek; \
455  R[0*_MODEL_NVARS_+4] = 1.0; \
456  R[1*_MODEL_NVARS_+4] = vx; \
457  R[2*_MODEL_NVARS_+4] = vy; \
458  R[3*_MODEL_NVARS_+4] = vz+a; \
459  R[4*_MODEL_NVARS_+4] = h0+a*vz; \
460  R[0*_MODEL_NVARS_+1] = 0.0; \
461  R[1*_MODEL_NVARS_+1] = -1.0; \
462  R[2*_MODEL_NVARS_+1] = 0.0; \
463  R[3*_MODEL_NVARS_+1] = 0.0; \
464  R[4*_MODEL_NVARS_+1] = -vx; \
465  R[0*_MODEL_NVARS_+2] = 0.0; \
466  R[1*_MODEL_NVARS_+2] = 0.0; \
467  R[2*_MODEL_NVARS_+2] = 1.0; \
468  R[3*_MODEL_NVARS_+2] = 0.0; \
469  R[4*_MODEL_NVARS_+2] = vy; \
470  } \
471  }
472 
482 typedef struct navierstokes3d_parameters {
483  double gamma;
484  double Re;
485  double Pr;
486  double Minf;
487  double C1,
488  C2;
489  double grav_x,
490  grav_y,
491  grav_z;
492  double rho0;
493  double p0;
494  double R;
495  char upw_choice[_MAX_STRING_SIZE_];
497  double *grav_field_f,
498  *grav_field_g;
500  double *fast_jac,
501  *solution;
504  char ib_wall_type[_MAX_STRING_SIZE_];
506  double T_ib_wall;
507 
508  /* choice of hydrostatic balance */
509  /* 1 -> isothermal */
510  /* 2 -> constant potential temperature */
511  /* 3 -> stratified atmosphere with a Brunt-Vaisala frequency */
512  int HB ;
515  double N_bv;
517  char ib_write_surface_data[_MAX_STRING_SIZE_];
522  double t_ib_ramp;
523 
526  double t_ib_width;
527 
531  char ib_ramp_type[_MAX_STRING_SIZE_];
532 
536  double ib_T_tol;
537 
538 #if defined(HAVE_CUDA)
539  double *gpu_Q;
540  double *gpu_QDerivX;
541  double *gpu_QDerivY;
542  double *gpu_QDerivZ;
543  double *gpu_FViscous;
544  double *gpu_FDeriv;
545  double *gpu_grav_field_f;
546  double *gpu_grav_field_g;
547  double *gpu_fast_jac;
548  double *gpu_solution;
549 #endif
550 } NavierStokes3D;
551 
552 int NavierStokes3DInitialize (void*,void*);
553 int NavierStokes3DCleanup (void*);
554 
555 #if defined(HAVE_CUDA)
556 int gpuNavierStokes3DCleanup (void*);
557 #endif
558 
559 static const int _NavierStokes3D_stride_ = 1;
NavierStokes3D::gpu_QDerivX
double * gpu_QDerivX
Definition: navierstokes3d.h:540
NavierStokes3DInitialize
int NavierStokes3DInitialize(void *, void *)
Definition: NavierStokes3DInitialize.c:116
NavierStokes3DCleanup
int NavierStokes3DCleanup(void *)
Definition: NavierStokes3DCleanup.c:11
NavierStokes3D::gpu_QDerivY
double * gpu_QDerivY
Definition: navierstokes3d.h:541
NavierStokes3D::gpu_grav_field_f
double * gpu_grav_field_f
Definition: navierstokes3d.h:545
NavierStokes3D::gpu_solution
double * gpu_solution
Definition: navierstokes3d.h:548
NavierStokes3D
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
Definition: navierstokes3d.h:482
basic.h
Some basic definitions and macros.
NavierStokes3D::grav_field_g
double * grav_field_g
Definition: navierstokes3d.h:497
gpuNavierStokes3DCleanup
int gpuNavierStokes3DCleanup(void *)
Definition: NavierStokes3DCleanup_GPU.c:13
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
_NavierStokes3D_stride_
static const int _NavierStokes3D_stride_
Definition: navierstokes3d.h:559
NavierStokes3D::gpu_fast_jac
double * gpu_fast_jac
Definition: navierstokes3d.h:547
NavierStokes3D::gpu_grav_field_g
double * gpu_grav_field_g
Definition: navierstokes3d.h:546
NavierStokes3D::gpu_FViscous
double * gpu_FViscous
Definition: navierstokes3d.h:543
NavierStokes3D::gpu_QDerivZ
double * gpu_QDerivZ
Definition: navierstokes3d.h:542
NavierStokes3D::gpu_FDeriv
double * gpu_FDeriv
Definition: navierstokes3d.h:544