HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Euler2DUpwind.c File Reference
#include <stdlib.h>
#include <math.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <physicalmodels/euler2d.h>
#include <mathfunctions.h>
#include <matmult_native.h>
#include <hypar.h>

Go to the source code of this file.

Functions

int Euler2DUpwindRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int Euler2DUpwindRF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int Euler2DUpwindLLF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int Euler2DUpwindSWFS (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 

Function Documentation

◆ Euler2DUpwindRoe()

int Euler2DUpwindRoe ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Definition at line 10 of file Euler2DUpwind.c.

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]++;
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 }
#define absolute(a)
Definition: math_ops.h:32
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
#define _Euler2DRightEigenvectors_(u, R, p, dir)
Definition: euler2d.h:189
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define MatVecMult4(N, y, A, x)
#define MatMult4(N, A, X, Y)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _Euler2DEigenvalues_(u, D, p, dir)
Definition: euler2d.h:108

◆ Euler2DUpwindRF()

int Euler2DUpwindRF ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Definition at line 68 of file Euler2DUpwind.c.

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]++;
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_],
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 */
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 */
132  }
133  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
134  }
135 
136  return(0);
137 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
#define _Euler2DRightEigenvectors_(u, R, p, dir)
Definition: euler2d.h:189
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define MatVecMult4(N, y, A, x)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _Euler2DEigenvalues_(u, D, p, dir)
Definition: euler2d.h:108

◆ Euler2DUpwindLLF()

int Euler2DUpwindLLF ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Definition at line 139 of file Euler2DUpwind.c.

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]++;
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 */
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 */
204  }
205  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
206  }
207 
208  return(0);
209 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
#define _Euler2DRightEigenvectors_(u, R, p, dir)
Definition: euler2d.h:189
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define MatVecMult4(N, y, A, x)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _Euler2DEigenvalues_(u, D, p, dir)
Definition: euler2d.h:108

◆ Euler2DUpwindSWFS()

int Euler2DUpwindSWFS ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Definition at line 211 of file Euler2DUpwind.c.

212 {
213  HyPar *solver = (HyPar*) s;
214  Euler2D *param = (Euler2D*) solver->physics;
215  int done,k;
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 
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 
279 
280  } else {
281 
283 
284  }
285 
286  }
287  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
288  }
289 
290  return(0);
291 }
int ndims
Definition: hypar.h:26
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _ArrayAdd1D_(x, a, b, size)
#define _YDIR_
Definition: euler2d.h:41
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
double gamma
Definition: euler2d.h:245
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define _XDIR_
Definition: euler1d.h:75
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _DECLARE_IERR_
Definition: basic.h:17
#define _ArrayCopy1D_(x, y, size)