HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NavierStokes3DGravityField.c
Go to the documentation of this file.
1 
5 #include <stdlib.h>
6 #include <math.h>
7 #include <arrayfunctions.h>
8 #include <mathfunctions.h>
10 #include <hypar.h>
11 #include <mpivars.h>
12 
34  void *s,
35  void *m
36  )
37 {
38  HyPar *solver = (HyPar*) s;
39  MPIVariables *mpi = (MPIVariables*) m;
40  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
41 
42  double *f = param->grav_field_f;
43  double *g = param->grav_field_g;
44  int *dim = solver->dim_local;
45  int ghosts = solver->ghosts;
46  int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_],
47  offset[_MODEL_NDIMS_], d, done;
48 
49  /* set bounds for array index to include ghost points */
50  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_);
51  for (d=0; d<_MODEL_NDIMS_; d++) bounds[d] += 2*ghosts;
52  /* set offset such that index is compatible with ghost point arrangement */
53  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
54 
55  double p0 = param->p0;
56  double rho0 = param->rho0;
57  double RT = p0 / rho0;
58  double gamma= param->gamma;
59  double R = param->R;
60  double Cp = gamma*R / (gamma-1.0);
61  double T0 = p0 / (rho0 * R);
62  double gx = param->grav_x;
63  double gy = param->grav_y;
64  double gz = param->grav_z;
65  double Nbv = param->N_bv;
66 
67  /* set the value of the gravity field */
68  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
69  if (param->HB == 1) {
70  while (!done) {
71  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
72  double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
73  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
74  double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
75  f[p] = exp( (gx*xcoord+gy*ycoord+gz*zcoord)/RT);
76  g[p] = exp(-(gx*xcoord+gy*ycoord+gz*zcoord)/RT);
77  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
78  }
79  } else if (param->HB == 2) {
80  while (!done) {
81  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
82  double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
83  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
84  double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
85  f[p] = raiseto((1.0-(gx*xcoord+gy*ycoord+gz*zcoord)/(Cp*T0)), (-1.0 /(gamma-1.0)));
86  g[p] = raiseto((1.0-(gx*xcoord+gy*ycoord+gz*zcoord)/(Cp*T0)), (gamma/(gamma-1.0)));
87  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
88  }
89  } else if (param->HB == 3) {
90  while (!done) {
91  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
92  double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
93  if ((gx != 0) || (gy != 0)) {
94  if (!mpi->rank) {
95  fprintf(stderr,"Error in NavierStokes3DGravityField(): HB = 3 is implemented only for ");
96  fprintf(stderr,"gravity force along the z-coordinate.\n");
97  }
98  return(1);
99  }
100  double Pexner = 1 + ((gz*gz)/(Cp*T0*Nbv*Nbv)) * (exp(-(Nbv*Nbv/gz)*zcoord)-1.0);
101  f[p] = raiseto(Pexner, (-1.0 /(gamma-1.0))) * exp(Nbv*Nbv*zcoord/gz);
102  g[p] = raiseto(Pexner, (gamma/(gamma-1.0)));
103  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
104  }
105  } else {
106  while (!done) {
107  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
108  f[p] = g[p] = 1.0;
109  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
110  }
111  }
112 
113  /* A sensible simulation will not specify peridic boundary conditions
114  * along a direction in which gravity acts.
115  * Gravity will be zero along the dimension periodic BCs are specified,
116  * so the value of the gravity potential will be the same along those
117  * grid lines.
118  * Thus, for both these cases, extrapolate the gravity field at the
119  * boundaries */
120  int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_];
121  for (d = 0; d < _MODEL_NDIMS_; d++) {
122  /* left boundary */
123  if (!mpi->ip[d]) {
124  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
125  _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
126  done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
127  while (!done) {
128  _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = ghosts-1-indexb[d];
129  int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
130  int p2; _ArrayIndex1D_ (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
131  f[p1] = f[p2];
132  g[p1] = g[p2];
133  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
134  }
135  }
136  /* right boundary */
137  if (mpi->ip[d] == mpi->iproc[d]-1) {
138  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
139  _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
140  done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
141  while (!done) {
142  _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = dim[d]-1-indexb[d];
143  int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
144  int p2; _ArrayIndex1D_ (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
145  f[p1] = f[p2];
146  g[p1] = g[p2];
147  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
148  }
149  }
150  }
151 
152  return(0);
153 }
#define _YDIR_
Definition: euler2d.h:41
#define _ArraySetValue_(x, size, value)
#define raiseto(x, a)
Definition: math_ops.h:37
3D Navier Stokes equations (compressible flows)
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
MPI related function definitions.
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ZDIR_
int NavierStokes3DGravityField(void *s, void *m)
double * grav_field_f
Contains function definitions for common mathematical functions.
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArrayCopy1D_(x, y, size)
Contains structure definition for hypar.
double * grav_field_g
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
Contains macros and function definitions for common array operations.
Structure of MPI-related variables.
double * x
Definition: hypar.h:107
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75