HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
NavierStokes2DInitialize
int NavierStokes2DInitialize(void *, void *)
Definition: NavierStokes2DInitialize.c:103
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
NavierStokes2D::gpu_grav_field_g
double * gpu_grav_field_g
Definition: navierstokes2d.h:405
NavierStokes2D::grav_field_g
double * grav_field_g
Definition: navierstokes2d.h:388
NavierStokes2D::gpu_solution
double * gpu_solution
Definition: navierstokes2d.h:407
NavierStokes2D::gpu_fast_jac
double * gpu_fast_jac
Definition: navierstokes2d.h:406
basic.h
Some basic definitions and macros.
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_f
double * gpu_grav_field_f
Definition: navierstokes2d.h:404
NavierStokes2DCleanup
int NavierStokes2DCleanup(void *)
Definition: NavierStokes2DCleanup.c:11