HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Euler2DUpwind.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <math.h>
3 #include <basic.h>
4 #include <arrayfunctions.h>
5 #include <physicalmodels/euler2d.h>
6 #include <mathfunctions.h>
7 #include <matmult_native.h>
8 #include <hypar.h>
9 
10 int Euler2DUpwindRoe(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
11 {
12  HyPar *solver = (HyPar*) s;
13  Euler2D *param = (Euler2D*) solver->physics;
14  int done;
15 
16  int *dim = solver->dim_local;
17 
18  int bounds_outer[2], bounds_inter[2];
19  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
20  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
21  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
22  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
23  modA[_MODEL_NVARS_*_MODEL_NVARS_];
24 
25  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
26  while (!done) {
27  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
28  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
29  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
30  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
31 
32  /* Roe's upwinding scheme */
33 
34  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
35  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
36  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
37  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
38 
39  _Euler2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param);
40 
41  _Euler2DEigenvalues_(uavg,D,param,dir);
42  _Euler2DLeftEigenvectors_(uavg,L,param,dir);
43  _Euler2DRightEigenvectors_(uavg,R,param,dir);
44 
45  /* Harten's Entropy Fix - Page 362 of Leveque */
46  int k;
47  double delta = 0.000001, delta2 = delta*delta;
48  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
49  k=5; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
50  k=10; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
51  k=15; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
52 
53  MatMult4(_MODEL_NVARS_,DL,D,L);
54  MatMult4(_MODEL_NVARS_,modA,R,DL);
55  MatVecMult4(_MODEL_NVARS_,udiss,modA,udiff);
56 
57  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
58  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
59  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
60  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
61  }
62  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
63  }
64 
65  return(0);
66 }
67 
68 int Euler2DUpwindRF(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
69 {
70  HyPar *solver = (HyPar*) s;
71  Euler2D *param = (Euler2D*) solver->physics;
72  int done,k;
73 
74  int *dim = solver->dim_local;
75 
76  int bounds_outer[2], bounds_inter[2];
77  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
78  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
79  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
80  L[_MODEL_NVARS_*_MODEL_NVARS_];
81 
82  done = 0; int index_outer[2] = {0,0}, index_inter[2];
83  while (!done) {
84  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
85  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
86  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
87  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
88  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
89 
90  /* Roe-Fixed upwinding scheme */
91 
92  _Euler2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param);
93 
94  _Euler2DEigenvalues_(uavg,D,param,dir);
95  _Euler2DLeftEigenvectors_ (uavg,L,param,dir);
96  _Euler2DRightEigenvectors_(uavg,R,param,dir);
97 
98  /* calculate characteristic fluxes and variables */
99  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
100  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
101  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
102  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
103 
104  double eigL[4],eigC[4],eigR[4];
105  _Euler2DEigenvalues_((uL+_MODEL_NVARS_*p),D,param,dir);
106  eigL[0] = D[0];
107  eigL[1] = D[5];
108  eigL[2] = D[10];
109  eigL[3] = D[15];
110  _Euler2DEigenvalues_((uR+_MODEL_NVARS_*p),D,param,dir);
111  eigR[0] = D[0];
112  eigR[1] = D[5];
113  eigR[2] = D[10];
114  eigR[3] = D[15];
115  _Euler2DEigenvalues_(uavg,D,param,dir);
116  eigC[0] = D[0];
117  eigC[1] = D[5];
118  eigC[2] = D[10];
119  eigC[3] = D[15];
120 
121  for (k = 0; k < _MODEL_NVARS_; k++) {
122  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
123  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
124  else {
125  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
126  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
127  }
128  }
129 
130  /* calculate the interface flux from the characteristic flux */
131  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
132  }
133  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
134  }
135 
136  return(0);
137 }
138 
139 int Euler2DUpwindLLF(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
140 {
141  HyPar *solver = (HyPar*) s;
142  Euler2D *param = (Euler2D*) solver->physics;
143  int done;
144 
145  int *dim = solver->dim_local;
146 
147  int bounds_outer[2], bounds_inter[2];
148  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
149  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
150  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
151  L[_MODEL_NVARS_*_MODEL_NVARS_];
152 
153  done = 0; int index_outer[2] = {0,0}, index_inter[2];
154  while (!done) {
155  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
156  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
157  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
158  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
159  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
160 
161  /* Local Lax-Friedrich upwinding scheme */
162 
163  _Euler2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param);
164 
165  _Euler2DEigenvalues_(uavg,D,param,dir);
166  _Euler2DLeftEigenvectors_ (uavg,L,param,dir);
167  _Euler2DRightEigenvectors_(uavg,R,param,dir);
168 
169  /* calculate characteristic fluxes and variables */
170  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
171  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
172  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
173  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
174 
175  double eigL[4],eigC[4],eigR[4];
176  _Euler2DEigenvalues_((uL+_MODEL_NVARS_*p),D,param,dir);
177  eigL[0] = D[0];
178  eigL[1] = D[5];
179  eigL[2] = D[10];
180  eigL[3] = D[15];
181  _Euler2DEigenvalues_((uR+_MODEL_NVARS_*p),D,param,dir);
182  eigR[0] = D[0];
183  eigR[1] = D[5];
184  eigR[2] = D[10];
185  eigR[3] = D[15];
186  _Euler2DEigenvalues_(uavg,D,param,dir);
187  eigC[0] = D[0];
188  eigC[1] = D[5];
189  eigC[2] = D[10];
190  eigC[3] = D[15];
191 
192  double alpha;
193  alpha = max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
194  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
195  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
196  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
197  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
198  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
199  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
200  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
201 
202  /* calculate the interface flux from the characteristic flux */
203  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
204  }
205  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
206  }
207 
208  return(0);
209 }
210 
211 int Euler2DUpwindSWFS(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
212 {
213  HyPar *solver = (HyPar*) s;
214  Euler2D *param = (Euler2D*) solver->physics;
215  int done,k;
216  _DECLARE_IERR_;
217 
218  int ndims = solver->ndims;
219  int *dim = solver->dim_local;
220 
221  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
222  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
223  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
224  static double fp[_MODEL_NVARS_], fm[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
225 
226  done = 0; _ArraySetValue_(index_outer,ndims,0);
227  while (!done) {
228  _ArrayCopy1D_(index_outer,index_inter,ndims);
229  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
230  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
231  double rho,vx,vy,e,P,c,gamma=param->gamma,term,Mach,lp[_MODEL_NVARS_],lm[_MODEL_NVARS_];
232 
233  /* Steger Warming flux splitting */
234  _Euler2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param);
235  _Euler2DGetFlowVar_(uavg,rho,vx,vy,e,P,param);
236  Mach = (dir==_XDIR_ ? vx : vy) / sqrt(gamma*P/rho);
237 
238  if (Mach < -1.0) {
239 
240  _ArrayCopy1D_((fR+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);
241 
242  } else if (Mach < 1.0) {
243 
244  double kx = 0, ky = 0;
245  kx = (dir==_XDIR_ ? 1.0 : 0.0);
246  ky = (dir==_YDIR_ ? 1.0 : 0.0);
247 
248  _Euler2DGetFlowVar_((uL+_MODEL_NVARS_*p),rho,vx,vy,e,P,param);
249  c = sqrt(gamma*P/rho);
250  term = rho/(2.0*gamma);
251  lp[0] = lp[1] = kx*vx + ky*vy;
252  lp[2] = lp[0] + c;
253  lp[3] = lp[0] - c;
254  for (k=0; k<_MODEL_NVARS_; k++) if (lp[k] < 0.0) lp[k] = 0.0;
255 
256  fp[0] = term * (2.0*(gamma-1.0)*lp[0] + lp[2] + lp[3]);
257  fp[1] = term * (2.0*(gamma-1.0)*lp[0]*vx + lp[2]*(vx+c*kx) + lp[3]*(vx-c*kx));
258  fp[2] = term * (2.0*(gamma-1.0)*lp[0]*vy + lp[2]*(vy+c*ky) + lp[3]*(vy-c*ky));
259  fp[3] = term * ((gamma-1.0)*lp[0]*(vx*vx+vy*vy) + 0.5*lp[2]*((vx+c*kx)*(vx+c*kx) + (vy+c*ky)*(vy+c*ky))
260  + 0.5*lp[3]*((vx-c*kx)*(vx-c*kx) + (vy-c*ky)*(vy-c*ky))
261  + ((3.0-gamma)*(lp[2]+lp[3])*c*c)/(2.0*(gamma-1.0)) );
262 
263  _Euler2DGetFlowVar_((uR+_MODEL_NVARS_*p),rho,vx,vy,e,P,param);
264  c = sqrt(gamma*P/rho);
265  term = rho/(2.0*gamma);
266  lm[0] = lm[1] = kx*vx + ky*vy;
267  lm[2] = lm[0] + c;
268  lm[3] = lm[0] - c;
269  for (k=0; k<_MODEL_NVARS_; k++) if (lm[k] > 0.0) lm[k] = 0.0;
270 
271  fm[0] = term * (2.0*(gamma-1.0)*lm[0] + lm[2] + lm[3]);
272  fm[1] = term * (2.0*(gamma-1.0)*lm[0]*vx + lm[2]*(vx+c*kx) + lm[3]*(vx-c*kx));
273  fm[2] = term * (2.0*(gamma-1.0)*lm[0]*vy + lm[2]*(vy+c*ky) + lm[3]*(vy-c*ky));
274  fm[3] = term * ((gamma-1.0)*lm[0]*(vx*vx+vy*vy) + 0.5*lm[2]*((vx+c*kx)*(vx+c*kx) + (vy+c*ky)*(vy+c*ky))
275  + 0.5*lm[3]*((vx-c*kx)*(vx-c*kx) + (vy-c*ky)*(vy-c*ky))
276  + ((3.0-gamma)*(lm[2]+lm[3])*c*c)/(2.0*(gamma-1.0)) );
277 
278  _ArrayAdd1D_((fI+_MODEL_NVARS_*p),fp,fm,_MODEL_NVARS_);
279 
280  } else {
281 
282  _ArrayCopy1D_((fL+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);
283 
284  }
285 
286  }
287  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
288  }
289 
290  return(0);
291 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
max3
#define max3(a, b, c)
Definition: math_ops.h:27
Euler2DUpwindLLF
int Euler2DUpwindLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:139
_Euler2DRoeAverage_
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
_Euler2DEigenvalues_
#define _Euler2DEigenvalues_(u, D, p, dir)
Definition: euler2d.h:108
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_Euler2DGetFlowVar_
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_Euler2DRightEigenvectors_
#define _Euler2DRightEigenvectors_(u, R, p, dir)
Definition: euler2d.h:189
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
Euler2DUpwindSWFS
int Euler2DUpwindSWFS(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:211
mathfunctions.h
Contains function definitions for common mathematical functions.
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
_Euler2DLeftEigenvectors_
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
Euler2DUpwindRoe
int Euler2DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:10
hypar.h
Contains structure definition for hypar.
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
basic.h
Some basic definitions and macros.
HyPar::ndims
int ndims
Definition: hypar.h:26
arrayfunctions.h
Contains macros and function definitions for common array operations.
absolute
#define absolute(a)
Definition: math_ops.h:32
Euler2DUpwindRF
int Euler2DUpwindRF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:68
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_ArrayAdd1D_
#define _ArrayAdd1D_(x, a, b, size)
Definition: arrayfunctions.h:212
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
matmult_native.h
Contains macros and function definitions for common matrix multiplication.
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
Euler2D::gamma
double gamma
Definition: euler2d.h:245