HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
navierstokes2d.h
Go to the documentation of this file.
1 
44 #include <basic.h>
45 
47 #define _NAVIER_STOKES_2D_ "navierstokes2d"
48 
49 /* define ndims and nvars for this model */
50 #undef _MODEL_NDIMS_
51 #undef _MODEL_NVARS_
52 
53 #define _MODEL_NDIMS_ 2
54 
55 #define _MODEL_NVARS_ 4
56 
57 /* choices for upwinding schemes */
59 #define _ROE_ "roe"
60 
61 #define _RF_ "rf-char"
62 
63 #define _LLF_ "llf-char"
64 
65 #define _SWFS_ "steger-warming"
66 
67 #define _RUSANOV_ "rusanov"
68 
69 /* directions */
71 #define _XDIR_ 0
72 
73 #define _YDIR_ 1
74 
81 #define _NavierStokes2DGetFlowVar_(u,rho,vx,vy,e,P,gamma) \
82  { \
83  double vsq; \
84  rho = u[0]; \
85  vx = (rho==0) ? 0 : u[1] / rho; \
86  vy = (rho==0) ? 0 : u[2] / rho; \
87  e = u[3]; \
88  vsq = (vx*vx) + (vy*vy); \
89  P = (e - 0.5*rho*vsq) * (gamma-1.0); \
90  }
91 
99 #define _NavierStokes2DSetFlux_(f,rho,vx,vy,e,P,dir) \
100  { \
101  if (dir == _XDIR_) { \
102  f[0] = rho * vx; \
103  f[1] = rho * vx * vx + P; \
104  f[2] = rho * vx * vy; \
105  f[3] = (e + P) * vx; \
106  } else if (dir == _YDIR_) { \
107  f[0] = rho * vy; \
108  f[1] = rho * vy * vx; \
109  f[2] = rho * vy * vy + P; \
110  f[3] = (e + P) * vy; \
111  } \
112  }
113 
125 #define _NavierStokes2DSetStiffFlux_(f,rho,vx,vy,e,P,dir,gamma) \
126  { \
127  double gamma_inv = 1.0/gamma; \
128  if (dir == _XDIR_) { \
129  f[0] = gamma_inv * rho * vx; \
130  f[1] = gamma_inv * rho * vx * vx + P; \
131  f[2] = gamma_inv * rho * vx * vy; \
132  f[3] = (e + P) * vx - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vx; \
133  } else if (dir == _YDIR_) { \
134  f[0] = gamma_inv * rho * vy; \
135  f[1] = gamma_inv * rho * vy * vx; \
136  f[2] = gamma_inv * rho * vy * vy + P; \
137  f[3] = (e + P) * vy - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vy; \
138  } \
139  }
140 
152 #define _NavierStokes2DSetNonStiffFlux_(f,rho,vx,vy,e,P,dir,gamma) \
153  { \
154  double gamma_inv = 1.0/gamma; \
155  if (dir == _XDIR_) { \
156  f[0] = (gamma-1.0) * gamma_inv * rho * vx; \
157  f[1] = (gamma-1.0) * gamma_inv * rho * vx * vx; \
158  f[2] = (gamma-1.0) * gamma_inv * rho * vx * vy; \
159  f[3] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vx; \
160  } else if (dir == _YDIR_) { \
161  f[0] = (gamma-1.0) * gamma_inv * rho * vy; \
162  f[1] = (gamma-1.0) * gamma_inv * rho * vy * vx; \
163  f[2] = (gamma-1.0) * gamma_inv * rho * vy * vy; \
164  f[3] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vy; \
165  } \
166  }
167 
171 #define _NavierStokes2DRoeAverage_(uavg,uL,uR,gamma) \
172  { \
173  double rho ,vx, vy, e ,P ,H ,csq, vsq; \
174  double rhoL,vxL,vyL,eL,PL,HL,cLsq,vsqL; \
175  double rhoR,vxR,vyR,eR,PR,HR,cRsq,vsqR; \
176  rhoL = uL[0]; \
177  vxL = uL[1] / rhoL; \
178  vyL = uL[2] / rhoL; \
179  eL = uL[3]; \
180  vsqL = (vxL*vxL) + (vyL*vyL); \
181  PL = (eL - 0.5*rhoL*vsqL) * (gamma-1.0); \
182  cLsq = gamma * PL/rhoL; \
183  HL = 0.5*(vxL*vxL+vyL*vyL) + cLsq / (gamma-1.0); \
184  rhoR = uR[0]; \
185  vxR = uR[1] / rhoR; \
186  vyR = uR[2] / rhoR; \
187  eR = uR[3]; \
188  vsqR = (vxR*vxR) + (vyR*vyR); \
189  PR = (eR - 0.5*rhoR*vsqR) * (gamma-1.0); \
190  cRsq = gamma * PR/rhoR; \
191  HR = 0.5*(vxR*vxR+vyR*vyR) + cRsq / (gamma-1.0); \
192  double tL = sqrt(rhoL); \
193  double tR = sqrt(rhoR); \
194  rho = tL * tR; \
195  vx = (tL*vxL + tR*vxR) / (tL + tR); \
196  vy = (tL*vyL + tR*vyR) / (tL + tR); \
197  H = (tL*HL + tR*HR) / (tL + tR); \
198  vsq = vx*vx + vy*vy; \
199  csq = (gamma-1.0) * (H-0.5*vsq); \
200  P = csq * rho / gamma; \
201  e = P/(gamma-1.0) + 0.5*rho*vsq; \
202  uavg[0] = rho; \
203  uavg[1] = rho*vx; \
204  uavg[2] = rho*vy; \
205  uavg[3] = e; \
206  }
207 
213 #define _NavierStokes2DEigenvalues_(u,D,gamma,dir) \
214  { \
215  double rho,vx,vy,e,P,c,vn,vsq; \
216  rho = u[0]; \
217  vx = u[1] / rho; \
218  vy = u[2] / rho; \
219  e = u[3]; \
220  vsq = (vx*vx) + (vy*vy); \
221  P = (e - 0.5*rho*vsq) * (gamma-1.0); \
222  c = sqrt(gamma*P/rho); \
223  if (dir == _XDIR_) vn = vx; \
224  else if (dir == _YDIR_) vn = vy; \
225  else vn = 0.0; \
226  if (dir == _XDIR_) {\
227  D[0*_MODEL_NVARS_+0] = vn-c; D[0*_MODEL_NVARS_+1] = 0; D[0*_MODEL_NVARS_+2] = 0; D[0*_MODEL_NVARS_+3] = 0; \
228  D[1*_MODEL_NVARS_+0] = 0; D[1*_MODEL_NVARS_+1] = vn+c; D[1*_MODEL_NVARS_+2] = 0; D[1*_MODEL_NVARS_+3] = 0; \
229  D[2*_MODEL_NVARS_+0] = 0; D[2*_MODEL_NVARS_+1] = 0; D[2*_MODEL_NVARS_+2] = vn; D[2*_MODEL_NVARS_+3] = 0; \
230  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vn; \
231  } else if (dir == _YDIR_) { \
232  D[0*_MODEL_NVARS_+0] = vn-c; D[0*_MODEL_NVARS_+1] = 0; D[0*_MODEL_NVARS_+2] = 0; D[0*_MODEL_NVARS_+3] = 0; \
233  D[1*_MODEL_NVARS_+0] = 0; D[1*_MODEL_NVARS_+1] = vn; D[1*_MODEL_NVARS_+2] = 0; D[1*_MODEL_NVARS_+3] = 0; \
234  D[2*_MODEL_NVARS_+0] = 0; D[2*_MODEL_NVARS_+1] = 0; D[2*_MODEL_NVARS_+2] = vn+c; D[2*_MODEL_NVARS_+3] = 0; \
235  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vn; \
236  }\
237  }
238 
247 #define _NavierStokes2DLeftEigenvectors_(u,L,ga,dir) \
248  { \
249  double ga_minus_one=ga-1.0; \
250  double rho,vx,vy,e,P,a,un,ek,vsq; \
251  double nx = 0,ny = 0; \
252  rho = u[0]; \
253  vx = u[1] / rho; \
254  vy = u[2] / rho; \
255  e = u[3]; \
256  vsq = (vx*vx) + (vy*vy); \
257  P = (e - 0.5*rho*vsq) * (ga-1.0); \
258  ek = 0.5 * (vx*vx + vy*vy); \
259  a = sqrt(ga * P / rho); \
260  if (dir == _XDIR_) { \
261  un = vx; \
262  nx = 1.0; \
263  L[0*_MODEL_NVARS_+0] = (ga_minus_one*ek + a*un) / (2*a*a); \
264  L[0*_MODEL_NVARS_+1] = ((-ga_minus_one)*vx - a*nx) / (2*a*a); \
265  L[0*_MODEL_NVARS_+2] = ((-ga_minus_one)*vy - a*ny) / (2*a*a); \
266  L[0*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
267  L[3*_MODEL_NVARS_+0] = (a*a - ga_minus_one*ek) / (a*a); \
268  L[3*_MODEL_NVARS_+1] = (ga_minus_one*vx) / (a*a); \
269  L[3*_MODEL_NVARS_+2] = (ga_minus_one*vy) / (a*a); \
270  L[3*_MODEL_NVARS_+3] = (-ga_minus_one) / (a*a); \
271  L[1*_MODEL_NVARS_+0] = (ga_minus_one*ek - a*un) / (2*a*a); \
272  L[1*_MODEL_NVARS_+1] = ((-ga_minus_one)*vx + a*nx) / (2*a*a); \
273  L[1*_MODEL_NVARS_+2] = ((-ga_minus_one)*vy + a*ny) / (2*a*a); \
274  L[1*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
275  L[2*_MODEL_NVARS_+0] = (vy - un*ny) / nx; \
276  L[2*_MODEL_NVARS_+1] = ny; \
277  L[2*_MODEL_NVARS_+2] = (ny*ny - 1.0) / nx; \
278  L[2*_MODEL_NVARS_+3] = 0.0; \
279  } else if (dir == _YDIR_) { \
280  un = vy; \
281  ny = 1.0; \
282  L[0*_MODEL_NVARS_+0] = (ga_minus_one*ek+a*un) / (2*a*a); \
283  L[0*_MODEL_NVARS_+1] = ((1.0-ga)*vx - a*nx) / (2*a*a); \
284  L[0*_MODEL_NVARS_+2] = ((1.0-ga)*vy - a*ny) / (2*a*a); \
285  L[0*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
286  L[3*_MODEL_NVARS_+0] = (a*a-ga_minus_one*ek) / (a*a); \
287  L[3*_MODEL_NVARS_+1] = ga_minus_one*vx / (a*a); \
288  L[3*_MODEL_NVARS_+2] = ga_minus_one*vy / (a*a); \
289  L[3*_MODEL_NVARS_+3] = (1.0 - ga) / (a*a); \
290  L[2*_MODEL_NVARS_+0] = (ga_minus_one*ek-a*un) / (2*a*a); \
291  L[2*_MODEL_NVARS_+1] = ((1.0-ga)*vx + a*nx) / (2*a*a); \
292  L[2*_MODEL_NVARS_+2] = ((1.0-ga)*vy + a*ny) / (2*a*a); \
293  L[2*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
294  L[1*_MODEL_NVARS_+0] = (un*nx-vx) / ny; \
295  L[1*_MODEL_NVARS_+1] = (1.0 - nx*nx) / ny; \
296  L[1*_MODEL_NVARS_+2] = - nx; \
297  L[1*_MODEL_NVARS_+3] = 0; \
298  } \
299  }
300 
309 #define _NavierStokes2DRightEigenvectors_(u,R,ga,dir) \
310  { \
311  double ga_minus_one = ga-1.0; \
312  double rho,vx,vy,e,P,un,ek,a,h0,vsq; \
313  double nx = 0,ny = 0; \
314  rho = u[0]; \
315  vx = u[1] / rho; \
316  vy = u[2] / rho; \
317  e = u[3]; \
318  vsq = (vx*vx) + (vy*vy); \
319  P = (e - 0.5*rho*vsq) * (ga-1.0); \
320  ek = 0.5 * (vx*vx + vy*vy); \
321  a = sqrt(ga * P / rho); \
322  h0 = a*a / ga_minus_one + ek; \
323  if (dir == _XDIR_) { \
324  un = vx; \
325  nx = 1.0; \
326  R[0*_MODEL_NVARS_+0] = 1.0; \
327  R[1*_MODEL_NVARS_+0] = vx - a*nx; \
328  R[2*_MODEL_NVARS_+0] = vy - a*ny; \
329  R[3*_MODEL_NVARS_+0] = h0 - a*un; \
330  R[0*_MODEL_NVARS_+3] = 1.0; \
331  R[1*_MODEL_NVARS_+3] = vx; \
332  R[2*_MODEL_NVARS_+3] = vy; \
333  R[3*_MODEL_NVARS_+3] = ek; \
334  R[0*_MODEL_NVARS_+1] = 1.0; \
335  R[1*_MODEL_NVARS_+1] = vx + a*nx; \
336  R[2*_MODEL_NVARS_+1] = vy + a*ny; \
337  R[3*_MODEL_NVARS_+1] = h0 + a*un; \
338  R[0*_MODEL_NVARS_+2] = 0.0; \
339  R[1*_MODEL_NVARS_+2] = ny; \
340  R[2*_MODEL_NVARS_+2] = -nx; \
341  R[3*_MODEL_NVARS_+2] = vx*ny - vy*nx; \
342  } else if (dir == _YDIR_) { \
343  un = vy; \
344  ny = 1.0; \
345  R[0*_MODEL_NVARS_+0] = 1.0; \
346  R[1*_MODEL_NVARS_+0] = vx - a*nx; \
347  R[2*_MODEL_NVARS_+0] = vy - a*ny; \
348  R[3*_MODEL_NVARS_+0] = h0 - a*un; \
349  R[0*_MODEL_NVARS_+3] = 1.0; \
350  R[1*_MODEL_NVARS_+3] = vx; \
351  R[2*_MODEL_NVARS_+3] = vy; \
352  R[3*_MODEL_NVARS_+3] = ek; \
353  R[0*_MODEL_NVARS_+2] = 1.0; \
354  R[1*_MODEL_NVARS_+2] = vx + a*nx; \
355  R[2*_MODEL_NVARS_+2] = vy + a*ny; \
356  R[3*_MODEL_NVARS_+2] = h0 + a*un; \
357  R[0*_MODEL_NVARS_+1] = 0; \
358  R[1*_MODEL_NVARS_+1] = ny; \
359  R[2*_MODEL_NVARS_+1] = -nx; \
360  R[3*_MODEL_NVARS_+1] = vx*ny-vy*nx; \
361  } \
362  }
363 
364 
374 typedef struct navierstokes2d_parameters {
375  double gamma;
376  char upw_choice[_MAX_STRING_SIZE_];
377  double grav_x,
378  grav_y;
379  double rho0;
380  double p0;
381  double Re;
382  double Pr;
383  double Minf;
384  double C1,
385  C2;
386  double R;
388  double *grav_field_f,
389  *grav_field_g;
391  double *fast_jac,
392  *solution;
394  /* choice of hydrostatic balance */
395  /* 1 -> isothermal */
396  /* 2 -> constant potential temperature */
397  /* 3 -> stratified atmosphere with a Brunt-Vaisala frequency */
398  int HB ;
401  double N_bv;
403 #if defined(HAVE_CUDA)
404  double *gpu_grav_field_f;
405  double *gpu_grav_field_g;
406  double *gpu_fast_jac;
407  double *gpu_solution;
408 #endif
409 } NavierStokes2D;
410 
411 int NavierStokes2DInitialize (void*,void*);
412 int NavierStokes2DCleanup (void*);
413 
NavierStokes2DCleanup
int NavierStokes2DCleanup(void *)
Definition: NavierStokes2DCleanup.c:11
basic.h
Some basic definitions and macros.
NavierStokes2D::grav_field_g
double * grav_field_g
Definition: navierstokes2d.h:388
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
NavierStokes2D::gpu_grav_field_f
double * gpu_grav_field_f
Definition: navierstokes2d.h:404
NavierStokes2D::gpu_solution
double * gpu_solution
Definition: navierstokes2d.h:407
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
NavierStokes2D::gpu_grav_field_g
double * gpu_grav_field_g
Definition: navierstokes2d.h:405
NavierStokes2DInitialize
int NavierStokes2DInitialize(void *, void *)
Definition: NavierStokes2DInitialize.c:103
NavierStokes2D::gpu_fast_jac
double * gpu_fast_jac
Definition: navierstokes2d.h:406