HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Numa2DUpwind.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 <mathfunctions.h>
6 #include <physicalmodels/numa2d.h>
7 #include <hypar.h>
8 
9 int Numa2DRusanovFlux(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
10 {
11  HyPar *solver = (HyPar*) s;
12  Numa2D *param = (Numa2D*) solver->physics;
13  int done;
14 
15  int *dim = solver->dim_local;
16  int ghosts = solver->ghosts;
17 
18  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
19  _ArrayCopy1D2_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
20  _ArrayCopy1D2_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
21 
22  done = 0; int index_outer[2] = {0,0}, index_inter[2];
23  while (!done) {
24  _ArrayCopy1D2_(index_outer,index_inter,_MODEL_NDIMS_);
25  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
26  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
27  double udiff[_MODEL_NVARS_],drho,vel[2],dT,rho0,P0,T0,EP,c;
28  double ycoordL, ycoordR;
29 
30  if (dir == _YDIR_) {
31  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_]-1),dim,ghosts,solver->x,ycoordL);
32  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordR);
33  } else {
34  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordL);
35  ycoordR = ycoordL;
36  }
37 
38  /* Rusanov's upwinding scheme */
39 
40  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
41  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
42  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
43  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
44 
45  /* left of the interface */
46  param->StandardAtmosphere (param,ycoordL,&EP,&P0,&rho0,&T0);
47  _Numa2DGetFlowVars_ ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
48  _Numa2DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
49  double alphaL = c + absolute(vel[dir]);
50 
51  /* right of the interface */
52  param->StandardAtmosphere (param,ycoordR,&EP,&P0,&rho0,&T0);
53  _Numa2DGetFlowVars_ ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
54  _Numa2DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
55  double alphaR = c + absolute(vel[dir]);
56 
57  double alpha = max(alphaL,alphaR);
58  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
59  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
60  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
61  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
62  }
63  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
64  }
65 
66  return(0);
67 }
68 
69 int Numa2DRusanovLinearFlux(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
70 {
71  HyPar *solver = (HyPar*) s;
72  Numa2D *param = (Numa2D*) solver->physics;
73  int done;
74 
75  int *dim = solver->dim_local;
76  int ghosts = solver->ghosts;
77 
78  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
79  _ArrayCopy1D2_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
80  _ArrayCopy1D2_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
81 
82  done = 0; int index_outer[2] = {0,0}, index_inter[2];
83  while (!done) {
84  _ArrayCopy1D2_(index_outer,index_inter,_MODEL_NDIMS_);
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 udiff[_MODEL_NVARS_],drho,vel[2],dT,rho0,P0,T0,EP,c;
88  double ycoordL, ycoordR;
89 
90  if (dir == _YDIR_) {
91  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_]-1),dim,ghosts,solver->x,ycoordL);
92  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordR);
93  } else {
94  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordL);
95  ycoordR = ycoordL;
96  }
97 
98  /* Rusanov's upwinding scheme */
99 
100  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
101  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
102  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
103  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
104 
105  /* left of the interface */
106  param->StandardAtmosphere (param,ycoordL,&EP,&P0,&rho0,&T0);
107  _Numa2DGetFlowVars_ ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
108  _Numa2DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
109  double alphaL = c;
110 
111  /* right of the interface */
112  param->StandardAtmosphere (param,ycoordR,&EP,&P0,&rho0,&T0);
113  _Numa2DGetFlowVars_ ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
114  _Numa2DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
115  double alphaR = c;
116 
117  double alpha = max(alphaL,alphaR);
118  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
119  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
120  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
121  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
122 
123  /* some harmless statements to avoid compiler warnings */
124  vel[0] = vel[1];
125  dT = vel[0];
126  vel[1] = dT;
127  }
128  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
129  }
130 
131  return(0);
132 }
133 
_Numa2DComputeSpeedofSound_
#define _Numa2DComputeSpeedofSound_(gamma, R, T0, dT, rho0, drho, EP, c)
Definition: numa2d.h:99
absolute
#define absolute(a)
Definition: math_ops.h:32
mathfunctions.h
Contains function definitions for common mathematical functions.
Numa2DRusanovLinearFlux
int Numa2DRusanovLinearFlux(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Numa2DUpwind.c:69
basic.h
Some basic definitions and macros.
HyPar::x
double * x
Definition: hypar.h:107
Numa2D::R
double R
Definition: numa2d.h:111
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Numa2D
Definition: numa2d.h:109
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
hypar.h
Contains structure definition for hypar.
Numa2D::gamma
double gamma
Definition: numa2d.h:110
_Numa2DComputeLinearizedSpeedofSound_
#define _Numa2DComputeLinearizedSpeedofSound_(gamma, R, T0, rho0, EP, c)
Definition: numa2d.h:104
_ArrayCopy1D2_
#define _ArrayCopy1D2_(x, y, size)
Definition: arrayfunctions.h:331
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
HyPar::physics
void * physics
Definition: hypar.h:266
_Numa2DGetFlowVars_
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
Numa2DRusanovFlux
int Numa2DRusanovFlux(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Numa2DUpwind.c:9
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
Numa2D::StandardAtmosphere
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
max
#define max(a, b)
Definition: math_ops.h:18
arrayfunctions.h
Contains macros and function definitions for common array operations.