HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
euler2d.h
Go to the documentation of this file.
1 /*
2 
3  2D Euler Equations for Inviscid, Compressible Flows
4 
5 
6  d [ rho ] d [ rho*u ] d [ rho*v ]
7  -- [ rho*u ] + -- [rho*u*u + p ] + -- [ rho*u*v ] = 0
8  dt [ rho*v ] dx [ rho*u*v ] dy [rho*v*v + p]
9  [ e ] [ (e+p)*u ] [ (e+p)v ]
10 
11  Equation of state:
12  p 1
13  e = ------- + - rho * (u^2 + v^2)
14  gamma-1 2
15 
16  Choices for upwinding:
17  "roe" Roe upwinding
18  "rf-char" Roe-fixed
19  "llf-char" Local Lax-Friedrich
20 
21 */
22 
23 #include <basic.h>
24 
25 #define _EULER_2D_ "euler2d"
26 
27 /* define ndims and nvars for this model */
28 #undef _MODEL_NDIMS_
29 #undef _MODEL_NVARS_
30 #define _MODEL_NDIMS_ 2
31 #define _MODEL_NVARS_ 4
32 
33 /* choices for upwinding schemes */
34 #define _ROE_ "roe"
35 #define _RF_ "rf-char"
36 #define _LLF_ "llf-char"
37 #define _SWFS_ "steger-warming"
38 
39 /* directions */
40 #define _XDIR_ 0
41 #define _YDIR_ 1
42 
43 
44 #define _Euler2DGetFlowVar_(u,rho,vx,vy,e,P,p) \
45  { \
46  double gamma = p->gamma, vsq; \
47  rho = u[0]; \
48  vx = u[1] / rho; \
49  vy = u[2] / rho; \
50  e = u[3]; \
51  vsq = (vx*vx) + (vy*vy); \
52  P = (e - 0.5*rho*vsq) * (gamma-1.0); \
53  }
54 
55 #define _Euler2DSetFlux_(f,rho,vx,vy,e,P,p,dir) \
56  { \
57  if (dir == _XDIR_) { \
58  f[0] = rho * vx; \
59  f[1] = rho * vx * vx + P; \
60  f[2] = rho * vx * vy; \
61  f[3] = (e + P) * vx; \
62  } else if (dir == _YDIR_) { \
63  f[0] = rho * vy; \
64  f[1] = rho * vy * vx; \
65  f[2] = rho * vy * vy + P; \
66  f[3] = (e + P) * vy; \
67  } \
68  }
69 
70 #define _Euler2DRoeAverage_(uavg,uL,uR,p) \
71  { \
72  double rho ,vx, vy, e ,P ,H ,csq, vsq; \
73  double rhoL,vxL,vyL,eL,PL,HL,cLsq,vsqL; \
74  double rhoR,vxR,vyR,eR,PR,HR,cRsq,vsqR; \
75  double gamma = p->gamma; \
76  rhoL = uL[0]; \
77  vxL = uL[1] / rhoL; \
78  vyL = uL[2] / rhoL; \
79  eL = uL[3]; \
80  vsqL = (vxL*vxL) + (vyL*vyL); \
81  PL = (eL - 0.5*rhoL*vsqL) * (gamma-1.0); \
82  cLsq = gamma * PL/rhoL; \
83  HL = 0.5*(vxL*vxL+vyL*vyL) + cLsq / (gamma-1.0); \
84  rhoR = uR[0]; \
85  vxR = uR[1] / rhoR; \
86  vyR = uR[2] / rhoR; \
87  eR = uR[3]; \
88  vsqR = (vxR*vxR) + (vyR*vyR); \
89  PR = (eR - 0.5*rhoR*vsqR) * (gamma-1.0); \
90  cRsq = gamma * PR/rhoR; \
91  HR = 0.5*(vxR*vxR+vyR*vyR) + cRsq / (gamma-1.0); \
92  double tL = sqrt(rhoL); \
93  double tR = sqrt(rhoR); \
94  rho = tL * tR; \
95  vx = (tL*vxL + tR*vxR) / (tL + tR); \
96  vy = (tL*vyL + tR*vyR) / (tL + tR); \
97  H = (tL*HL + tR*HR) / (tL + tR); \
98  vsq = vx*vx + vy*vy; \
99  csq = (gamma-1.0) * (H-0.5*vsq); \
100  P = csq * rho / gamma; \
101  e = P/(gamma-1.0) + 0.5*rho*vsq; \
102  uavg[0] = rho; \
103  uavg[1] = rho*vx; \
104  uavg[2] = rho*vy; \
105  uavg[3] = e; \
106  }
107 
108 #define _Euler2DEigenvalues_(u,D,p,dir) \
109  { \
110  double gamma = p->gamma; \
111  double rho,vx,vy,e,P,c,vn,vsq; \
112  rho = u[0]; \
113  vx = u[1] / rho; \
114  vy = u[2] / rho; \
115  e = u[3]; \
116  vsq = (vx*vx) + (vy*vy); \
117  P = (e - 0.5*rho*vsq) * (gamma-1.0); \
118  c = sqrt(gamma*P/rho); \
119  if (dir == _XDIR_) vn = vx; \
120  else if (dir == _YDIR_) vn = vy; \
121  else vn = 0.0; \
122  if (dir == _XDIR_) {\
123  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; \
124  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; \
125  D[2*_MODEL_NVARS_+0] = 0; D[2*_MODEL_NVARS_+1] = 0; D[2*_MODEL_NVARS_+2] = vn; D[2*_MODEL_NVARS_+3] = 0; \
126  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vn; \
127  } else if (dir == _YDIR_) { \
128  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; \
129  D[1*_MODEL_NVARS_+0] = 0; D[1*_MODEL_NVARS_+1] = vn; D[1*_MODEL_NVARS_+2] = 0; D[1*_MODEL_NVARS_+3] = 0; \
130  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; \
131  D[3*_MODEL_NVARS_+0] = 0; D[3*_MODEL_NVARS_+1] = 0; D[3*_MODEL_NVARS_+2] = 0; D[3*_MODEL_NVARS_+3] = vn; \
132  }\
133  }
134 
135 #define _Euler2DLeftEigenvectors_(u,L,p,dir) \
136  { \
137  double ga = param->gamma, ga_minus_one=ga-1.0; \
138  double rho,vx,vy,e,P,a,un,ek,vsq; \
139  double nx = 0,ny = 0; \
140  rho = u[0]; \
141  vx = u[1] / rho; \
142  vy = u[2] / rho; \
143  e = u[3]; \
144  vsq = (vx*vx) + (vy*vy); \
145  P = (e - 0.5*rho*vsq) * (ga-1.0); \
146  ek = 0.5 * (vx*vx + vy*vy); \
147  a = sqrt(ga * P / rho); \
148  if (dir == _XDIR_) { \
149  un = vx; \
150  nx = 1.0; \
151  L[0*_MODEL_NVARS_+0] = (ga_minus_one*ek + a*un) / (2*a*a); \
152  L[0*_MODEL_NVARS_+1] = ((-ga_minus_one)*vx - a*nx) / (2*a*a); \
153  L[0*_MODEL_NVARS_+2] = ((-ga_minus_one)*vy - a*ny) / (2*a*a); \
154  L[0*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
155  L[3*_MODEL_NVARS_+0] = (a*a - ga_minus_one*ek) / (a*a); \
156  L[3*_MODEL_NVARS_+1] = (ga_minus_one*vx) / (a*a); \
157  L[3*_MODEL_NVARS_+2] = (ga_minus_one*vy) / (a*a); \
158  L[3*_MODEL_NVARS_+3] = (-ga_minus_one) / (a*a); \
159  L[1*_MODEL_NVARS_+0] = (ga_minus_one*ek - a*un) / (2*a*a); \
160  L[1*_MODEL_NVARS_+1] = ((-ga_minus_one)*vx + a*nx) / (2*a*a); \
161  L[1*_MODEL_NVARS_+2] = ((-ga_minus_one)*vy + a*ny) / (2*a*a); \
162  L[1*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
163  L[2*_MODEL_NVARS_+0] = (vy - un*ny) / nx; \
164  L[2*_MODEL_NVARS_+1] = ny; \
165  L[2*_MODEL_NVARS_+2] = (ny*ny - 1.0) / nx; \
166  L[2*_MODEL_NVARS_+3] = 0.0; \
167  } else if (dir == _YDIR_) { \
168  un = vy; \
169  ny = 1.0; \
170  L[0*_MODEL_NVARS_+0] = (ga_minus_one*ek+a*un) / (2*a*a); \
171  L[0*_MODEL_NVARS_+1] = ((1.0-ga)*vx - a*nx) / (2*a*a); \
172  L[0*_MODEL_NVARS_+2] = ((1.0-ga)*vy - a*ny) / (2*a*a); \
173  L[0*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
174  L[3*_MODEL_NVARS_+0] = (a*a-ga_minus_one*ek) / (a*a); \
175  L[3*_MODEL_NVARS_+1] = ga_minus_one*vx / (a*a); \
176  L[3*_MODEL_NVARS_+2] = ga_minus_one*vy / (a*a); \
177  L[3*_MODEL_NVARS_+3] = (1.0 - ga) / (a*a); \
178  L[2*_MODEL_NVARS_+0] = (ga_minus_one*ek-a*un) / (2*a*a); \
179  L[2*_MODEL_NVARS_+1] = ((1.0-ga)*vx + a*nx) / (2*a*a); \
180  L[2*_MODEL_NVARS_+2] = ((1.0-ga)*vy + a*ny) / (2*a*a); \
181  L[2*_MODEL_NVARS_+3] = ga_minus_one / (2*a*a); \
182  L[1*_MODEL_NVARS_+0] = (un*nx-vx) / ny; \
183  L[1*_MODEL_NVARS_+1] = (1.0 - nx*nx) / ny; \
184  L[1*_MODEL_NVARS_+2] = - nx; \
185  L[1*_MODEL_NVARS_+3] = 0; \
186  } \
187  }
188 
189 #define _Euler2DRightEigenvectors_(u,R,p,dir) \
190  { \
191  double ga = param->gamma, ga_minus_one = ga-1.0; \
192  double rho,vx,vy,e,P,un,ek,a,h0,vsq; \
193  double nx = 0,ny = 0; \
194  rho = u[0]; \
195  vx = u[1] / rho; \
196  vy = u[2] / rho; \
197  e = u[3]; \
198  vsq = (vx*vx) + (vy*vy); \
199  P = (e - 0.5*rho*vsq) * (ga-1.0); \
200  ek = 0.5 * (vx*vx + vy*vy); \
201  a = sqrt(ga * P / rho); \
202  h0 = a*a / ga_minus_one + ek; \
203  if (dir == _XDIR_) { \
204  un = vx; \
205  nx = 1.0; \
206  R[0*_MODEL_NVARS_+0] = 1.0; \
207  R[1*_MODEL_NVARS_+0] = vx - a*nx; \
208  R[2*_MODEL_NVARS_+0] = vy - a*ny; \
209  R[3*_MODEL_NVARS_+0] = h0 - a*un; \
210  R[0*_MODEL_NVARS_+3] = 1.0; \
211  R[1*_MODEL_NVARS_+3] = vx; \
212  R[2*_MODEL_NVARS_+3] = vy; \
213  R[3*_MODEL_NVARS_+3] = ek; \
214  R[0*_MODEL_NVARS_+1] = 1.0; \
215  R[1*_MODEL_NVARS_+1] = vx + a*nx; \
216  R[2*_MODEL_NVARS_+1] = vy + a*ny; \
217  R[3*_MODEL_NVARS_+1] = h0 + a*un; \
218  R[0*_MODEL_NVARS_+2] = 0.0; \
219  R[1*_MODEL_NVARS_+2] = ny; \
220  R[2*_MODEL_NVARS_+2] = -nx; \
221  R[3*_MODEL_NVARS_+2] = vx*ny - vy*nx; \
222  } else if (dir == _YDIR_) { \
223  un = vy; \
224  ny = 1.0; \
225  R[0*_MODEL_NVARS_+0] = 1.0; \
226  R[1*_MODEL_NVARS_+0] = vx - a*nx; \
227  R[2*_MODEL_NVARS_+0] = vy - a*ny; \
228  R[3*_MODEL_NVARS_+0] = h0 - a*un; \
229  R[0*_MODEL_NVARS_+3] = 1.0; \
230  R[1*_MODEL_NVARS_+3] = vx; \
231  R[2*_MODEL_NVARS_+3] = vy; \
232  R[3*_MODEL_NVARS_+3] = ek; \
233  R[0*_MODEL_NVARS_+2] = 1.0; \
234  R[1*_MODEL_NVARS_+2] = vx + a*nx; \
235  R[2*_MODEL_NVARS_+2] = vy + a*ny; \
236  R[3*_MODEL_NVARS_+2] = h0 + a*un; \
237  R[0*_MODEL_NVARS_+1] = 0; \
238  R[1*_MODEL_NVARS_+1] = ny; \
239  R[2*_MODEL_NVARS_+1] = -nx; \
240  R[3*_MODEL_NVARS_+1] = vx*ny-vy*nx; \
241  } \
242  }
243 
244 typedef struct euler2d_parameters {
245  double gamma; /* Ratio of heat capacities */
246  char upw_choice[_MAX_STRING_SIZE_]; /* choice of upwinding */
247 } Euler2D;
248 
249 int Euler2DInitialize (void*,void*);
250 int Euler2DCleanup (void*);
251 
Euler2DInitialize
int Euler2DInitialize(void *, void *)
Definition: Euler2DInitialize.c:21
Euler2DCleanup
int Euler2DCleanup(void *)
Definition: Euler2DCleanup.c:4
basic.h
Some basic definitions and macros.
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
Euler2D::gamma
double gamma
Definition: euler2d.h:245