HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
Euler2DInitialize.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <boundaryconditions.h>
#include <physicalmodels/euler2d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions

double Euler2DComputeCFL (void *, void *, double, double)
 
int Euler2DFlux (double *, double *, int, void *, double)
 
int Euler2DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int Euler2DUpwindRF (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int Euler2DUpwindLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int Euler2DUpwindSWFS (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int Euler2DRoeAverage (double *, double *, double *, void *)
 
int Euler2DLeftEigenvectors (double *, double *, void *, int)
 
int Euler2DRightEigenvectors (double *, double *, void *, int)
 
int Euler2DInitialize (void *s, void *m)
 

Function Documentation

double Euler2DComputeCFL ( void *  ,
void *  ,
double  ,
double   
)

Definition at line 9 of file Euler2DComputeCFL.c.

10 {
11  HyPar *solver = (HyPar*) s;
12  Euler2D *param = (Euler2D*) solver->physics;
13  _DECLARE_IERR_;
14 
15  int *dim = solver->dim_local;
16  int ghosts = solver->ghosts;
17  int ndims = solver->ndims;
18  int index[ndims];
19  double *u = solver->u;
20 
21  double max_cfl = 0;
22  int done = 0; _ArraySetValue_(index,ndims,0);
23  while (!done) {
24  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
25  double rho,vx,vy,e,P,c,dxinv,dyinv,local_cfl[2];
26  _Euler2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,vx,vy,e,P,param);
27 
28  c = sqrt(param->gamma*P/rho); /* speed of sound */
29  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
30  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv); /* 1/dy */
31 
32  local_cfl[_XDIR_] = (absolute(vx)+c)*dt*dxinv; /* local cfl for this grid point (x) */
33  local_cfl[_YDIR_] = (absolute(vy)+c)*dt*dyinv; /* local cfl for this grid point (y) */
34  if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
35  if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];
36 
37  _ArrayIncrementIndex_(ndims,dim,index,done);
38  }
39 
40  return(max_cfl);
41 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::u
double * u
Definition: hypar.h:116
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
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
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
HyPar::dxinv
double * dxinv
Definition: hypar.h:110
HyPar::ndims
int ndims
Definition: hypar.h:26
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
Euler2D::gamma
double gamma
Definition: euler2d.h:245
int Euler2DFlux ( double *  ,
double *  ,
int  ,
void *  ,
double   
)

Definition at line 7 of file Euler2DFlux.c.

8 {
9  HyPar *solver = (HyPar*) s;
10  Euler2D *param = (Euler2D*) solver->physics;
11  int i;
12 
13  int *dim = solver->dim_local;
14  int ghosts = solver->ghosts;
15  int ndims = solver->ndims;
16 
17  int index[ndims], bounds[ndims], offset[ndims];
18 
19  /* set bounds for array index to include ghost points */
20  _ArrayCopy1D_(dim,bounds,ndims);
21  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
22 
23  /* set offset such that index is compatible with ghost point arrangement */
24  _ArraySetValue_(offset,ndims,-ghosts);
25 
26  int done = 0; _ArraySetValue_(index,ndims,0);
27  while (!done) {
28  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
29  double rho, vx, vy, e, P;
30  _Euler2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,vx,vy,e,P,param);
31  _Euler2DSetFlux_((f+_MODEL_NVARS_*p),rho,vx,vy,e,P,param,dir);
32  _ArrayIncrementIndex_(ndims,bounds,index,done);
33  }
34 
35  return(0);
36 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
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
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
_Euler2DSetFlux_
#define _Euler2DSetFlux_(f, rho, vx, vy, e, P, p, dir)
Definition: euler2d.h:55
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int Euler2DUpwindRoe ( double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int  ,
void *  ,
double   
)

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]++;
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 }
_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
_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
_Euler2DLeftEigenvectors_
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int Euler2DUpwindRF ( double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int  ,
void *  ,
double   
)

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]++;
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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
_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
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_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
_Euler2DLeftEigenvectors_
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int Euler2DUpwindLLF ( double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int  ,
void *  ,
double   
)

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]++;
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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
_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
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_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
_Euler2DLeftEigenvectors_
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int Euler2DUpwindSWFS ( double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int  ,
void *  ,
double   
)

Definition at line 211 of file Euler2DUpwind.c.

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
_Euler2DRoeAverage_
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
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
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::ndims
int ndims
Definition: hypar.h:26
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
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
Euler2D::gamma
double gamma
Definition: euler2d.h:245
int Euler2DRoeAverage ( double *  ,
double *  ,
double *  ,
void *   
)

Definition at line 5 of file Euler2DFunctions.c.

6 {
7  Euler2D *param = (Euler2D*) p;
8  _Euler2DRoeAverage_(uavg,uL,uR,param);
9  return(0);
10 }
_Euler2DRoeAverage_
#define _Euler2DRoeAverage_(uavg, uL, uR, p)
Definition: euler2d.h:70
int Euler2DLeftEigenvectors ( double *  ,
double *  ,
void *  ,
int   
)

Definition at line 16 of file Euler2DEigen.c.

17 {
18  Euler2D *param = (Euler2D*) p;
19  _Euler2DLeftEigenvectors_(u,L,param,dir);
20  return(0);
21 }
_Euler2DLeftEigenvectors_
#define _Euler2DLeftEigenvectors_(u, L, p, dir)
Definition: euler2d.h:135
int Euler2DRightEigenvectors ( double *  ,
double *  ,
void *  ,
int   
)

Definition at line 23 of file Euler2DEigen.c.

24 {
25  Euler2D *param = (Euler2D*) p;
26  _Euler2DRightEigenvectors_(u,R,param,dir);
27  return(0);
28 }
_Euler2DRightEigenvectors_
#define _Euler2DRightEigenvectors_(u, R, p, dir)
Definition: euler2d.h:189
int Euler2DInitialize ( void *  s,
void *  m 
)

Definition at line 21 of file Euler2DInitialize.c.

22 {
23  HyPar *solver = (HyPar*) s;
24  MPIVariables *mpi = (MPIVariables*) m;
25  Euler2D *physics = (Euler2D*) solver->physics;
26  int ferr = 0;
27 
28  static int count = 0;
29 
30  if (solver->nvars != _MODEL_NVARS_) {
31  fprintf(stderr,"Error in Euler2DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
32  return(1);
33  }
34  if (solver->ndims != _MODEL_NDIMS_) {
35  fprintf(stderr,"Error in Euler2DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
36  return(1);
37  }
38 
39  /* default values */
40  physics->gamma = 1.4;
41  strcpy(physics->upw_choice,"roe");
42 
43  /* reading physical model specific inputs - all processes */
44  if (!mpi->rank) {
45  FILE *in;
46  if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
47  in = fopen("physics.inp","r");
48  if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
49  else {
50  char word[_MAX_STRING_SIZE_];
51  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
52  if (!strcmp(word, "begin")){
53  while (strcmp(word, "end")){
54  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
55  if (!strcmp(word, "gamma")) {
56  ferr = fscanf(in,"%lf",&physics->gamma); if (ferr != 1) return(1);
57  } else if (!strcmp(word,"upwinding")) {
58  ferr = fscanf(in,"%s",physics->upw_choice); if (ferr != 1) return(1);
59  } else if (strcmp(word,"end")) {
60  char useless[_MAX_STRING_SIZE_];
61  ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
62  printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
63  printf("recognized or extraneous. Ignoring.\n");
64  }
65  }
66  } else {
67  fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
68  return(1);
69  }
70  }
71  fclose(in);
72  }
73 
74 #ifndef serial
75  IERR MPIBroadcast_double (&physics->gamma,1,0,&mpi->world); CHECKERR(ierr);
76  IERR MPIBroadcast_character (physics->upw_choice,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
77 #endif
78 
79  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
80  if (!mpi->rank) {
81  fprintf(stderr,"Error in Euler2DInitialize: This physical model does not have a splitting ");
82  fprintf(stderr,"of the hyperbolic term defined.\n");
83  }
84  return(1);
85  }
86 
87  /* initializing physical model-specific functions */
88  solver->ComputeCFL = Euler2DComputeCFL;
89  solver->FFunction = Euler2DFlux;
90  if (!strcmp(physics->upw_choice,_ROE_ )) solver->Upwind = Euler2DUpwindRoe;
91  else if (!strcmp(physics->upw_choice,_RF_ )) solver->Upwind = Euler2DUpwindRF;
92  else if (!strcmp(physics->upw_choice,_LLF_ )) solver->Upwind = Euler2DUpwindLLF;
93  else if (!strcmp(physics->upw_choice,_SWFS_)) solver->Upwind = Euler2DUpwindSWFS;
94  else {
95  fprintf(stderr,"Error in Euler2DInitialize(): %s is not a valid upwinding scheme.\n",
96  physics->upw_choice);
97  return(1);
98  }
99  solver->AveragingFunction = Euler2DRoeAverage;
100  solver->GetLeftEigenvectors = Euler2DLeftEigenvectors;
101  solver->GetRightEigenvectors = Euler2DRightEigenvectors;
102 
103  /* set the value of gamma in all the boundary objects */
104  int n;
105  DomainBoundary *boundary = (DomainBoundary*) solver->boundary;
106  for (n = 0; n < solver->nBoundaryZones; n++) boundary[n].gamma = physics->gamma;
107 
108  count++;
109  return(0);
110 }
HyPar::ComputeCFL
double(* ComputeCFL)(void *, void *, double, double)
Definition: hypar.h:269
Euler2DUpwindLLF
int Euler2DUpwindLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:139
MPIBroadcast_double
int MPIBroadcast_double(double *, int, int, void *)
Definition: MPIBroadcast.c:9
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::Upwind
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
HyPar::GetRightEigenvectors
int(* GetRightEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:359
_LLF_
#define _LLF_
Definition: euler1d.h:66
HyPar::FFunction
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_SWFS_
#define _SWFS_
Definition: euler1d.h:68
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
Euler2DLeftEigenvectors
int Euler2DLeftEigenvectors(double *u, double *L, void *p, int dir)
Definition: Euler2DEigen.c:16
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
HyPar::AveragingFunction
int(* AveragingFunction)(double *, double *, double *, void *)
Definition: hypar.h:354
HyPar::boundary
void * boundary
Definition: hypar.h:159
DomainBoundary
Structure containing the variables and function pointers defining a boundary.
Definition: boundaryconditions.h:68
Euler2DComputeCFL
double Euler2DComputeCFL(void *s, void *m, double dt, double t)
Definition: Euler2DComputeCFL.c:9
Euler2DUpwindSWFS
int Euler2DUpwindSWFS(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:211
Euler2DRightEigenvectors
int Euler2DRightEigenvectors(double *u, double *R, void *p, int dir)
Definition: Euler2DEigen.c:23
HyPar::GetLeftEigenvectors
int(* GetLeftEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:357
_ROE_
#define _ROE_
Definition: euler1d.h:62
HyPar::nBoundaryZones
int nBoundaryZones
Definition: hypar.h:157
HyPar::nvars
int nvars
Definition: hypar.h:29
Euler2DUpwindRoe
int Euler2DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:10
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
Euler2DRoeAverage
int Euler2DRoeAverage(double *uavg, double *uL, double *uR, void *p)
Definition: Euler2DFunctions.c:5
MPIBroadcast_character
int MPIBroadcast_character(char *, int, int, void *)
Definition: MPIBroadcast.c:37
HyPar::SplitHyperbolicFlux
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
HyPar::ndims
int ndims
Definition: hypar.h:26
_RF_
#define _RF_
Definition: euler1d.h:64
IERR
#define IERR
Definition: basic.h:16
Euler2DUpwindRF
int Euler2DUpwindRF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Euler2DUpwind.c:68
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
Euler2D::upw_choice
char upw_choice[_MAX_STRING_SIZE_]
Definition: euler2d.h:246
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Euler2DFlux
int Euler2DFlux(double *f, double *u, int dir, void *s, double t)
Definition: Euler2DFlux.c:7
Euler2D::gamma
double gamma
Definition: euler2d.h:245