HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NavierStokes3DParabolicFunction.c
Go to the documentation of this file.
1 
5 #include <string.h>
6 #include <stdlib.h>
7 #include <basic.h>
8 #include <arrayfunctions.h>
9 #include <mathfunctions.h>
11 #include <mpivars.h>
12 #include <hypar.h>
13 
14 #if defined(CPU_STAT)
15 #include <time.h>
16 #endif
17 
51  double *par,
52  double *u,
53  void *s,
54  void *m,
55  double t
56  )
57 {
58  HyPar *solver = (HyPar*) s;
59  MPIVariables *mpi = (MPIVariables*) m;
60  NavierStokes3D *physics = (NavierStokes3D*) solver->physics;
61  int i,j,k,v;
63 
64  int ghosts = solver->ghosts;
65  int imax = solver->dim_local[0];
66  int jmax = solver->dim_local[1];
67  int kmax = solver->dim_local[2];
68  int *dim = solver->dim_local;
69  int size = solver->npoints_local_wghosts;
70 
71  _ArraySetValue_(par,size*_MODEL_NVARS_,0.0);
72  if (physics->Re <= 0) return(0); /* inviscid flow */
73  solver->count_par++;
74 
75  static double two_third = 2.0/3.0;
76  double inv_gamma_m1 = 1.0 / (physics->gamma-1.0);
77  double inv_Re = 1.0 / physics->Re;
78  double inv_Pr = 1.0 / physics->Pr;
79 
80 #if defined(CPU_STAT)
81  clock_t startEvent, stopEvent;
82 #endif
83 
84  double *Q; /* primitive variables */
85  Q = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
86 
87 #if defined(CPU_STAT)
88  startEvent = clock();
89 #endif
90 
91  for (i=-ghosts; i<(imax+ghosts); i++) {
92  for (j=-ghosts; j<(jmax+ghosts); j++) {
93  for (k=-ghosts; k<(kmax+ghosts); k++) {
94  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
95  double energy,pressure;
96  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
99  Q[p+0],
100  Q[p+1],
101  Q[p+2],
102  Q[p+3],
103  energy,
104  pressure,
105  physics->gamma);
106  Q[p+4] = physics->gamma*pressure/Q[p+0]; /* temperature */
107  }
108  }
109  }
110 
111 #if defined(CPU_STAT)
112  stopEvent = clock();
113  printf("%-50s CPU time (secs) = %.6f\n",
114  "NaverStokes3DParabolicFunction_Q",(double)(stopEvent-startEvent)/CLOCKS_PER_SEC);
115 #endif
116 
117  double *QDerivX = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
118  double *QDerivY = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
119  double *QDerivZ = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
120 
121  IERR solver->FirstDerivativePar(QDerivX,Q,_XDIR_,1,solver,mpi); CHECKERR(ierr);
122  IERR solver->FirstDerivativePar(QDerivY,Q,_YDIR_,1,solver,mpi); CHECKERR(ierr);
123  IERR solver->FirstDerivativePar(QDerivZ,Q,_ZDIR_,1,solver,mpi); CHECKERR(ierr);
124 
125  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
126  solver->ghosts,mpi,QDerivX); CHECKERR(ierr);
127  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
128  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
129  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
130  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
131 
132  for (i=-ghosts; i<(imax+ghosts); i++) {
133  for (j=-ghosts; j<(jmax+ghosts); j++) {
134  for (k=-ghosts; k<(kmax+ghosts); k++) {
135  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
136  double dxinv, dyinv, dzinv;
137  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
138  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
139  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
140  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
141  _ArrayScale1D_((QDerivX+p),dxinv,_MODEL_NVARS_);
142  _ArrayScale1D_((QDerivY+p),dyinv,_MODEL_NVARS_);
143  _ArrayScale1D_((QDerivZ+p),dzinv,_MODEL_NVARS_);
144  }
145  }
146  }
147 
148  double *FViscous = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
149  double *FDeriv = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
150 
151  /* Along X */
152  for (i=-ghosts; i<(imax+ghosts); i++) {
153  for (j=-ghosts; j<(jmax+ghosts); j++) {
154  for (k=-ghosts; k<(kmax+ghosts); k++) {
155  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
156  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
157 
158  double uvel, vvel, wvel, T, Tx,
159  ux, uy, uz, vx, vy, wx, wz;
160  uvel = (Q+p)[1];
161  vvel = (Q+p)[2];
162  wvel = (Q+p)[3];
163  T = (Q+p)[4];
164  Tx = (QDerivX+p)[4];
165  ux = (QDerivX+p)[1];
166  vx = (QDerivX+p)[2];
167  wx = (QDerivX+p)[3];
168  uy = (QDerivY+p)[1];
169  vy = (QDerivY+p)[2];
170  uz = (QDerivZ+p)[1];
171  wz = (QDerivZ+p)[3];
172 
173  /* calculate viscosity coeff based on Sutherland's law */
174  double mu = raiseto(T, 0.76);
175  /*
176  if (p/_MODEL_NVARS_ == 49841) {
177  printf("[CPU] 49841 mu = %.20e T = %.20e log(0.76) = %.20e T*log(0.76) = %.20e\n",
178  mu, T, log(0.76), T*log(0.76));
179  }
180  */
181 
182  double tau_xx, tau_xy, tau_xz, qx;
183  tau_xx = two_third * (mu*inv_Re) * (2*ux - vy - wz);
184  tau_xy = (mu*inv_Re) * (uy + vx);
185  tau_xz = (mu*inv_Re) * (uz + wx);
186  qx = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tx;
187 
188  (FViscous+p)[0] = 0.0;
189  (FViscous+p)[1] = tau_xx;
190  (FViscous+p)[2] = tau_xy;
191  (FViscous+p)[3] = tau_xz;
192  (FViscous+p)[4] = uvel*tau_xx + vvel*tau_xy + wvel*tau_xz + qx;
193  }
194  }
195  }
196 
197  IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
198  for (i=0; i<imax; i++) {
199  for (j=0; j<jmax; j++) {
200  for (k=0; k<kmax; k++) {
201  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
202  double dxinv;
203  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
204  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
205  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dxinv * (FDeriv+p)[v] );
206  }
207  }
208  }
209 
210  /* Along Y */
211  for (i=-ghosts; i<(imax+ghosts); i++) {
212  for (j=-ghosts; j<(jmax+ghosts); j++) {
213  for (k=-ghosts; k<(kmax+ghosts); k++) {
214  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
215  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
216 
217  double uvel, vvel, wvel, T, Ty,
218  ux, uy, vx, vy, vz, wy, wz;
219  uvel = (Q+p)[1];
220  vvel = (Q+p)[2];
221  wvel = (Q+p)[3];
222  T = (Q+p)[4];
223  Ty = (QDerivY+p)[4];
224  ux = (QDerivX+p)[1];
225  vx = (QDerivX+p)[2];
226  uy = (QDerivY+p)[1];
227  vy = (QDerivY+p)[2];
228  wy = (QDerivY+p)[3];
229  vz = (QDerivZ+p)[2];
230  wz = (QDerivZ+p)[3];
231 
232  /* calculate viscosity coeff based on Sutherland's law */
233  double mu = raiseto(T, 0.76);
234 
235  double tau_yx, tau_yy, tau_yz, qy;
236  tau_yx = (mu*inv_Re) * (uy + vx);
237  tau_yy = two_third * (mu*inv_Re) * (-ux + 2*vy - wz);
238  tau_yz = (mu*inv_Re) * (vz + wy);
239  qy = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Ty;
240 
241  (FViscous+p)[0] = 0.0;
242  (FViscous+p)[1] = tau_yx;
243  (FViscous+p)[2] = tau_yy;
244  (FViscous+p)[3] = tau_yz;
245  (FViscous+p)[4] = uvel*tau_yx + vvel*tau_yy + wvel*tau_yz + qy;
246  }
247  }
248  }
249 
250  IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
251  for (i=0; i<imax; i++) {
252  for (j=0; j<jmax; j++) {
253  for (k=0; k<kmax; k++) {
254  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
255  double dyinv;
256  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
257  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
258  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
259  }
260  }
261  }
262 
263  /* Along Z */
264  for (i=-ghosts; i<(imax+ghosts); i++) {
265  for (j=-ghosts; j<(jmax+ghosts); j++) {
266  for (k=-ghosts; k<(kmax+ghosts); k++) {
267  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
268  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
269 
270  double uvel, vvel, wvel, T, Tz,
271  ux, uz, vy, vz, wx, wy, wz;
272  uvel = (Q+p)[1];
273  vvel = (Q+p)[2];
274  wvel = (Q+p)[3];
275  T = (Q+p)[4];
276  Tz = (QDerivZ+p)[4];
277  ux = (QDerivX+p)[1];
278  wx = (QDerivX+p)[3];
279  vy = (QDerivY+p)[2];
280  wy = (QDerivY+p)[3];
281  uz = (QDerivZ+p)[1];
282  vz = (QDerivZ+p)[2];
283  wz = (QDerivZ+p)[3];
284 
285  /* calculate viscosity coeff based on Sutherland's law */
286  double mu = raiseto(T,0.76);
287 
288  double tau_zx, tau_zy, tau_zz, qz;
289  tau_zx = (mu*inv_Re) * (uz + wx);
290  tau_zy = (mu*inv_Re) * (vz + wy);
291  tau_zz = two_third * (mu*inv_Re) * (-ux - vy + 2*wz);
292  qz = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tz;
293 
294  (FViscous+p)[0] = 0.0;
295  (FViscous+p)[1] = tau_zx;
296  (FViscous+p)[2] = tau_zy;
297  (FViscous+p)[3] = tau_zz;
298  (FViscous+p)[4] = uvel*tau_zx + vvel*tau_zy + wvel*tau_zz + qz;
299  }
300  }
301  }
302 
303  IERR solver->FirstDerivativePar(FDeriv,FViscous,_ZDIR_,-1,solver,mpi); CHECKERR(ierr);
304  for (i=0; i<imax; i++) {
305  for (j=0; j<jmax; j++) {
306  for (k=0; k<kmax; k++) {
307  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
308  double dzinv;
309  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
310  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
311  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dzinv * (FDeriv+p)[v] );
312  }
313  }
314  }
315 
316  free(Q);
317  free(QDerivX);
318  free(QDerivY);
319  free(QDerivZ);
320  free(FViscous);
321  free(FDeriv);
322 
323  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,_MODEL_NVARS_);
324  return(0);
325 }
#define _YDIR_
Definition: euler2d.h:41
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
#define raiseto(x, a)
Definition: math_ops.h:37
3D Navier Stokes equations (compressible flows)
double * iblank
Definition: hypar.h:436
#define _MODEL_NVARS_
Definition: euler1d.h:58
void * physics
Definition: hypar.h:266
int flag_ib
Definition: hypar.h:441
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
#define _ArrayScale1D_(x, a, size)
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ZDIR_
int count_par
Definition: hypar.h:418
Contains function definitions for common mathematical functions.
#define _ArrayBlockMultiply_(x, a, n, bs)
int(* FirstDerivativePar)(double *, double *, int, int, void *, void *)
Definition: hypar.h:243
#define CHECKERR(ierr)
Definition: basic.h:18
Contains structure definition for hypar.
double * dxinv
Definition: hypar.h:110
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.
Some basic definitions and macros.
Contains macros and function definitions for common array operations.
int NavierStokes3DParabolicFunction(double *, double *, void *, void *, double)
#define IERR
Definition: basic.h:16
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75
#define _DECLARE_IERR_
Definition: basic.h:17
static const int _NavierStokes3D_stride_