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
Numa3DUpwind.c File Reference
#include <stdlib.h>
#include <math.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <physicalmodels/numa3d.h>
#include <hypar.h>

Go to the source code of this file.

Functions

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

Function Documentation

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

Definition at line 9 of file Numa3DUpwind.c.

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 }
_ArrayIndex1D3_
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:69
_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
Numa3D::R
double R
Definition: numa3d.h:130
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayCopy1D3_
#define _ArrayCopy1D3_(x, y, size)
Definition: arrayfunctions.h:343
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
Numa3D
Definition: numa3d.h:128
_Numa3DComputeSpeedofSound_
#define _Numa3DComputeSpeedofSound_(gamma, R, T0, dT, rho0, drho, EP, c)
Definition: numa3d.h:118
Numa3D::StandardAtmosphere
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa3d.h:139
Numa3D::gamma
double gamma
Definition: numa3d.h:129
absolute
#define absolute(a)
Definition: math_ops.h:32
max
#define max(a, b)
Definition: math_ops.h:18
_Numa3DGetFlowVars_
#define _Numa3DGetFlowVars_(u, drho, uvel, vvel, wvel, dT, rho0)
Definition: numa3d.h:39
HyPar::x
double * x
Definition: hypar.h:107
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int Numa3DRusanovLinearFlux ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Definition at line 71 of file Numa3DUpwind.c.

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 }
_ArrayIndex1D3_
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:69
_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
Numa3D::R
double R
Definition: numa3d.h:130
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayCopy1D3_
#define _ArrayCopy1D3_(x, y, size)
Definition: arrayfunctions.h:343
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
Numa3D
Definition: numa3d.h:128
Numa3D::StandardAtmosphere
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa3d.h:139
Numa3D::gamma
double gamma
Definition: numa3d.h:129
_Numa3DComputeLinearizedSpeedofSound_
#define _Numa3DComputeLinearizedSpeedofSound_(gamma, R, T0, rho0, EP, c)
Definition: numa3d.h:123
absolute
#define absolute(a)
Definition: math_ops.h:32
max
#define max(a, b)
Definition: math_ops.h:18
_Numa3DGetFlowVars_
#define _Numa3DGetFlowVars_(u, drho, uvel, vvel, wvel, dT, rho0)
Definition: numa3d.h:39
HyPar::x
double * x
Definition: hypar.h:107
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23