HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Numa3DUpwind.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/numa3d.h>
7 #include <hypar.h>
8 
9 int Numa3DRusanovFlux(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  Numa3D *param = (Numa3D*) 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  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
20  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
21 
22  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
23  while (!done) {
24  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
25  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
26  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
27  double udiff[_MODEL_NVARS_],drho,vel[3],dT,rho0,P0,T0,EP,c;
28  double zcoordL, zcoordR;
29 
30  if (dir == _ZDIR_) {
31  _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]-1),dim,ghosts,solver->x,zcoordL);
32  _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_] ),dim,ghosts,solver->x,zcoordR);
33  } else {
34  _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_] ),dim,ghosts,solver->x,zcoordL);
35  zcoordR = zcoordL;
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  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
45 
46  /* left of the interface */
47  param->StandardAtmosphere(param,zcoordL,&EP,&P0,&rho0,&T0);
48  _Numa3DGetFlowVars_ ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
49  _Numa3DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
50  double alphaL = c + absolute(vel[dir]);
51 
52  /* right of the interface */
53  param->StandardAtmosphere(param,zcoordR,&EP,&P0,&rho0,&T0);
54  _Numa3DGetFlowVars_ ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
55  _Numa3DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
56  double alphaR = c + absolute(vel[dir]);
57 
58  double alpha = max(alphaL,alphaR);
59  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
60  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
61  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
62  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
63  fI[_MODEL_NVARS_*p+4] = 0.5*(fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4])-alpha*udiff[4];
64  }
65  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
66  }
67 
68  return(0);
69 }
70 
71 int Numa3DRusanovLinearFlux(double *fI,double *fL,double *fR,double *uL,double *uR,double *u,int dir,void *s,double t)
72 {
73  HyPar *solver = (HyPar*) s;
74  Numa3D *param = (Numa3D*) solver->physics;
75  int done;
76 
77  int *dim = solver->dim_local;
78  int ghosts = solver->ghosts;
79 
80  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
81  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
82  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
83 
84  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
85  while (!done) {
86  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
87  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
88  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
89  double udiff[_MODEL_NVARS_],drho,vel[3],dT,rho0,P0,T0,EP,c;
90  double zcoordL, zcoordR;
91 
92  if (dir == _ZDIR_) {
93  _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_]-1),dim,ghosts,solver->x,zcoordL);
94  _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_] ),dim,ghosts,solver->x,zcoordR);
95  } else {
96  _GetCoordinate_(_ZDIR_,(index_inter[_ZDIR_] ),dim,ghosts,solver->x,zcoordL);
97  zcoordR = zcoordL;
98  }
99 
100  /* Rusanov's upwinding scheme */
101 
102  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
103  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
104  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
105  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
106  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
107 
108  /* left of the interface */
109  param->StandardAtmosphere (param,zcoordL,&EP,&P0,&rho0,&T0);
110  _Numa3DGetFlowVars_ ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
111  _Numa3DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
112  double alphaL = c + absolute(vel[dir]);
113 
114  /* right of the interface */
115  param->StandardAtmosphere(param,zcoordR,&EP,&P0,&rho0,&T0);
116  _Numa3DGetFlowVars_ ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],vel[2],dT,rho0);
117  _Numa3DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
118  double alphaR = c + absolute(vel[dir]);
119 
120  double alpha = max(alphaL,alphaR);
121  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
122  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
123  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
124  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
125  fI[_MODEL_NVARS_*p+4] = 0.5*(fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4])-alpha*udiff[4];
126 
127  vel[0] = dT; /* useless statement to avoid compiler warning */
128  }
129  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
130  }
131 
132  return(0);
133 }
absolute
#define absolute(a)
Definition: math_ops.h:32
mathfunctions.h
Contains function definitions for common mathematical functions.
_Numa3DComputeLinearizedSpeedofSound_
#define _Numa3DComputeLinearizedSpeedofSound_(gamma, R, T0, rho0, EP, c)
Definition: numa3d.h:123
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
_Numa3DGetFlowVars_
#define _Numa3DGetFlowVars_(u, drho, uvel, vvel, wvel, dT, rho0)
Definition: numa3d.h:39
basic.h
Some basic definitions and macros.
HyPar::x
double * x
Definition: hypar.h:107
Numa3DRusanovLinearFlux
int Numa3DRusanovLinearFlux(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Numa3DUpwind.c:71
Numa3D
Definition: numa3d.h:128
_ArrayIndex1D3_
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:69
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Numa3DRusanovFlux
int Numa3DRusanovFlux(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Numa3DUpwind.c:9
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
hypar.h
Contains structure definition for hypar.
_ArrayCopy1D3_
#define _ArrayCopy1D3_(x, y, size)
Definition: arrayfunctions.h:343
Numa3D::R
double R
Definition: numa3d.h:130
_Numa3DComputeSpeedofSound_
#define _Numa3DComputeSpeedofSound_(gamma, R, T0, dT, rho0, drho, EP, c)
Definition: numa3d.h:118
Numa3D::gamma
double gamma
Definition: numa3d.h:129
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
Numa3D::StandardAtmosphere
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa3d.h:139
_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
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
max
#define max(a, b)
Definition: math_ops.h:18
arrayfunctions.h
Contains macros and function definitions for common array operations.