HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NavierStokes3DUpwind.c
Go to the documentation of this file.
1 
5 #include <stdlib.h>
6 #include <math.h>
7 #include <basic.h>
8 #include <arrayfunctions.h>
9 #include <mathfunctions.h>
10 #include <matmult_native.h>
11 #include <physicalmodels/navierstokes3d.h>
12 #include <hypar.h>
13 
14 #if defined(CPU_STAT)
15 #include <time.h>
16 #endif
17 
18 static const int dummy = 1;
19 
40 int NavierStokes3DUpwindRoe(
41  double *fI,
42  double *fL,
43  double *fR,
44  double *uL,
45  double *uR,
46  double *u,
47  int dir,
48  void *s,
49  double t
50  )
51 {
52  HyPar *solver = (HyPar*) s;
53  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
54  int done;
55 
56  int *dim = solver->dim_local;
57 
58  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
59  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
60  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
61  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
62  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
63  modA[_MODEL_NVARS_*_MODEL_NVARS_];
64 
65  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
66 
67 #if defined(CPU_STAT)
68  clock_t startEvent, stopEvent;
69  startEvent = clock();
70 #endif
71 
72  while (!done) {
73  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
74  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
75  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
76  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
77  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
78  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
79  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
80  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
81 
82  /* Roe's upwinding scheme */
83 
84  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
85  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
86  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
87  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
88  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
89 
90  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
91  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
92  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
93  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
94 
95  /* Harten's Entropy Fix - Page 362 of Leveque */
96  int k;
97  double delta = 0.000001, delta2 = delta*delta;
98  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
99  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
100  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
101  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
102  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
103 
104  MatMult5(_MODEL_NVARS_,DL,D,L);
105  MatMult5(_MODEL_NVARS_,modA,R,DL);
106  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
107 
108  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
109  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
110  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
111  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
112  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
113  }
114  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
115  }
116 
117 #if defined(CPU_STAT)
118  stopEvent = clock();
119  printf("%-50s CPU time (secs) = %.6f dir = %d\n",
120  "NavierStokes3DUpwindRoe", (double)(stopEvent-startEvent)/CLOCKS_PER_SEC, dir);
121 #endif
122 
123  return(0);
124 }
125 
140 int NavierStokes3DUpwindRF(
141  double *fI,
142  double *fL,
143  double *fR,
144  double *uL,
145  double *uR,
146  double *u,
147  int dir,
148  void *s,
149  double t
150  )
151 {
152  HyPar *solver = (HyPar*) s;
153  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
154  int done,k;
155 
156  int *dim = solver->dim_local;
157 
158  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
159  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
160  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
161  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
162 
163  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
164  while (!done) {
165  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
166  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
167  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
168  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
169  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
170 
171  /* Roe-Fixed upwinding scheme */
172 
173  _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
174 
175  _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
176  _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
177 
178  /* calculate characteristic fluxes and variables */
179  MatVecMult5(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
180  MatVecMult5(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
181  MatVecMult5(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
182  MatVecMult5(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
183 
184  double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
185  _NavierStokes3DEigenvalues_((uL+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
186  eigL[0] = D[0];
187  eigL[1] = D[6];
188  eigL[2] = D[12];
189  eigL[3] = D[18];
190  eigL[4] = D[24];
191  _NavierStokes3DEigenvalues_((uR+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
192  eigR[0] = D[0];
193  eigR[1] = D[6];
194  eigR[2] = D[12];
195  eigR[3] = D[18];
196  eigR[4] = D[24];
197  _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
198  eigC[0] = D[0];
199  eigC[1] = D[6];
200  eigC[2] = D[12];
201  eigC[3] = D[18];
202  eigC[4] = D[24];
203 
204  for (k = 0; k < _MODEL_NVARS_; k++) {
205  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
206  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
207  else {
208  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
209  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
210  }
211  }
212 
213  /* calculate the interface flux from the characteristic flux */
214  MatVecMult5(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
215  }
216  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
217  }
218 
219  return(0);
220 }
221 
244 int NavierStokes3DUpwindLLF(
245  double *fI,
246  double *fL,
247  double *fR,
248  double *uL,
249  double *uR,
250  double *u,
251  int dir,
252  void *s,
253  double t
254  )
255 {
256  HyPar *solver = (HyPar*) s;
257  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
258  int done;
259 
260  int *dim = solver->dim_local;
261 
262  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
263  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
264  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
265  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
266 
267  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
268  while (!done) {
269  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
270  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
271  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
272  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
273  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
274 
275  /* Local Lax-Friedrich upwinding scheme */
276 
277  _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
278 
279  _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
280  _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
281 
282  /* calculate characteristic fluxes and variables */
283  MatVecMult5(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
284  MatVecMult5(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
285  MatVecMult5(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
286  MatVecMult5(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
287 
288  double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
289  _NavierStokes3DEigenvalues_((uL+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
290  eigL[0] = D[0];
291  eigL[1] = D[6];
292  eigL[2] = D[12];
293  eigL[3] = D[18];
294  eigL[4] = D[24];
295  _NavierStokes3DEigenvalues_((uR+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
296  eigR[0] = D[0];
297  eigR[1] = D[6];
298  eigR[2] = D[12];
299  eigR[3] = D[18];
300  eigR[4] = D[24];
301  _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
302  eigC[0] = D[0];
303  eigC[1] = D[6];
304  eigC[2] = D[12];
305  eigC[3] = D[18];
306  eigC[4] = D[24];
307 
308  double alpha;
309  alpha = max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
310  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
311  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
312  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
313  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
314  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
315  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
316  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
317  alpha = max3(absolute(eigL[4]),absolute(eigC[4]),absolute(eigR[4]));
318  fc[4] = 0.5 * (fcL[4] + fcR[4] + alpha * (ucL[4]-ucR[4]));
319 
320  /* calculate the interface flux from the characteristic flux */
321  MatVecMult5(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
322  }
323  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
324  }
325 
326  return(0);
327 }
328 
349 int NavierStokes3DUpwindRusanov(
350  double *fI,
351  double *fL,
352  double *fR,
353  double *uL,
354  double *uR,
355  double *u,
356  int dir,
357  void *s,
358  double t
359  )
360 {
361  HyPar *solver = (HyPar*) s;
362  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
363  int *dim = solver->dim_local, done;
364 
365  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
366  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
367  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
368 
369  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
370 
371  static int index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
372  indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_];
373 
374  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
375  while (!done) {
376  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
377  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
378 
379  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
380  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
381  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
382  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
383  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
384  int q = p*_MODEL_NVARS_;
385 
386  /* Modified Rusanov's upwinding scheme */
387 
388  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
389  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
390  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
391  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
392  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
393 
394  _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
395 
396  double c, vel[_MODEL_NDIMS_], rho,E,P;
397  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
398  c = sqrt(param->gamma*P/rho);
399  double alphaL = c + absolute(vel[dir]);
400  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
401  c = sqrt(param->gamma*P/rho);
402  double alphaR = c + absolute(vel[dir]);
403  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
404  c = sqrt(param->gamma*P/rho);
405  double alphaavg = c + absolute(vel[dir]);
406 
407  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
408  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
409 
410  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-alpha*udiff[0];
411  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-alpha*udiff[1];
412  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-alpha*udiff[2];
413  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-alpha*udiff[3];
414  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-alpha*udiff[4];
415  }
416  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
417  }
418 
419  return(0);
420 }
421 
430 int NavierStokes3DUpwinddFRoe(
431  double *fI,
432  double *fL,
433  double *fR,
434  double *uL,
435  double *uR,
436  double *u,
437  int dir,
438  void *s,
439  double t
440  )
441 {
442  HyPar *solver = (HyPar*) s;
443  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
444  int done;
445 
446  int *dim = solver->dim_local;
447  double *uref = param->solution;
448 
449  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
450  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
451  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
452  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
453  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
454  modA[_MODEL_NVARS_*_MODEL_NVARS_];
455 
456  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
457  while (!done) {
458  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
459  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
460  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
461  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
462  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
463  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
464  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
465  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
466 
467  /* Roe's upwinding scheme */
468 
469  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
470  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
471  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
472  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
473  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
474 
475  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
476  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
477  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
478  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
479 
480  /* Harten's Entropy Fix - Page 362 of Leveque */
481  int k;
482  double delta = 0.000001, delta2 = delta*delta;
483  if (dir == _XDIR_) {
484  k=0; D[k] = 0.0;
485  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
486  k=12; D[k] = 0.0;
487  k=18; D[k] = 0.0;
488  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
489  } else if (dir == _YDIR_) {
490  k=0; D[k] = 0.0;
491  k=6; D[k] = 0.0;
492  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
493  k=18; D[k] = 0.0;
494  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
495  } else if (dir == _ZDIR_) {
496  k=0; D[k] = 0.0;
497  k=6; D[k] = 0.0;
498  k=12; D[k] = 0.0;
499  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
500  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
501  }
502 
503  MatMult5(_MODEL_NVARS_,DL,D,L);
504  MatMult5(_MODEL_NVARS_,modA,R,DL);
505  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
506 
507  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
508  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
509  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
510  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
511  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
512  }
513  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
514  }
515 
516  return(0);
517 }
518 
527 int NavierStokes3DUpwindFdFRoe(
528  double *fI,
529  double *fL,
530  double *fR,
531  double *uL,
532  double *uR,
533  double *u,
534  int dir,
535  void *s,
536  double t
537  )
538 {
539  HyPar *solver = (HyPar*) s;
540  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
541  int done;
542 
543  int *dim = solver->dim_local;
544  double *uref = param->solution;
545 
546  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
547  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
548  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
549  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
550  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
551  modA[_MODEL_NVARS_*_MODEL_NVARS_];
552 
553  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
554  while (!done) {
555  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
556  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
557  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
558  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
559  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
560  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
561  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
562  int k;
563  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_],
564  udiss_total[_MODEL_NVARS_],udiss_stiff[_MODEL_NVARS_];
565  double delta = 0.000001, delta2 = delta*delta;
566 
567  /* Roe's upwinding scheme */
568 
569  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
570  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
571  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
572  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
573  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
574 
575  /* Compute total dissipation */
576  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
577  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
578  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
579  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
580  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
581  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
582  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
583  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
584  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
585  MatMult5(_MODEL_NVARS_,DL,D,L);
586  MatMult5(_MODEL_NVARS_,modA,R,DL);
587  MatVecMult5(_MODEL_NVARS_,udiss_total,modA,udiff);
588 
589  /* Compute dissipation corresponding to acoustic modes */
590  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
591  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
592  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
593  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
594  if (dir == _XDIR_) {
595  k=0; D[k] = 0.0;
596  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
597  k=12; D[k] = 0.0;
598  k=18; D[k] = 0.0;
599  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
600  } else if (dir == _YDIR_) {
601  k=0; D[k] = 0.0;
602  k=6; D[k] = 0.0;
603  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
604  k=18; D[k] = 0.0;
605  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
606  } else if (dir == _ZDIR_) {
607  k=0; D[k] = 0.0;
608  k=6; D[k] = 0.0;
609  k=12; D[k] = 0.0;
610  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
611  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
612  }
613  MatMult5(_MODEL_NVARS_,DL,D,L);
614  MatMult5(_MODEL_NVARS_,modA,R,DL);
615  MatVecMult5(_MODEL_NVARS_,udiss_stiff,modA,udiff);
616 
617  /* Compute the dissipation term for the entropy modes */
618  _ArraySubtract1D_(udiss,udiss_total,udiss_stiff,_MODEL_NVARS_);
619 
620  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
621  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
622  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
623  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
624  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
625  }
626  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
627  }
628 
629  return(0);
630 }
631 
638 int NavierStokes3DUpwindRusanovModified(
639  double *fI,
640  double *fL,
641  double *fR,
642  double *uL,
643  double *uR,
644  double *u,
645  int dir,
646  void *s,
647  double t
648  )
649 {
650  HyPar *solver = (HyPar*) s;
651  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
652  int *dim = solver->dim_local, done;
653 
654  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
655  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
656  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
657 
658  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
659  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
660  modA[_MODEL_NVARS_*_MODEL_NVARS_];
661 
662  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
663  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
664 
665  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
666 
667  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
668  while (!done) {
669  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
670  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
671 
672  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
673  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
674  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
675  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
676  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
677  int q = p*_MODEL_NVARS_;
678 
679  /* Modified Rusanov's upwinding scheme */
680 
681  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
682  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
683  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
684  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
685  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
686 
687  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
688  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
689  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
690 
691  double c, vel[_MODEL_NDIMS_], rho,E,P;
692 
693  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
694  c = sqrt(param->gamma*P/rho);
695  double alphaL = c + absolute(vel[dir]);
696  double betaL = absolute(vel[dir]);
697 
698  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
699  c = sqrt(param->gamma*P/rho);
700  double alphaR = c + absolute(vel[dir]);
701  double betaR = absolute(vel[dir]);
702 
703  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
704  c = sqrt(param->gamma*P/rho);
705  double alphaavg = c + absolute(vel[dir]);
706  double betaavg = absolute(vel[dir]);
707 
708  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
709  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
710  double beta = kappa*max3(betaL,betaR,betaavg);
711 
712  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
713  if (dir == _XDIR_) {
714  D[0] = beta;
715  D[6] = alpha;
716  D[12] = beta;
717  D[18] = beta;
718  D[24] = alpha;
719  } else if (dir == _YDIR_) {
720  D[0] = beta;
721  D[6] = beta;
722  D[12] = alpha;
723  D[18] = beta;
724  D[24] = alpha;
725  } else if (dir == _ZDIR_) {
726  D[0] = beta;
727  D[6] = beta;
728  D[12] = beta;
729  D[18] = alpha;
730  D[24] = alpha;
731  }
732  MatMult5 (_MODEL_NVARS_,DL,D,L);
733  MatMult5 (_MODEL_NVARS_,modA,R,DL);
734  MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
735 
736  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
737  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
738  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
739  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
740  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
741  }
742  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
743  }
744 
745  return(0);
746 }
747 
756 int NavierStokes3DUpwinddFRusanovModified(
757  double *fI,
758  double *fL,
759  double *fR,
760  double *uL,
761  double *uR,
762  double *u,
763  int dir,
764  void *s,
765  double t
766  )
767 {
768  HyPar *solver = (HyPar*) s;
769  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
770  int *dim = solver->dim_local, done;
771  double *uref = param->solution;
772 
773  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
774  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
775  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
776 
777  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
778  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
779  modA[_MODEL_NVARS_*_MODEL_NVARS_];
780 
781  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
782  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
783 
784  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
785 
786  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
787  while (!done) {
788  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
789  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
790 
791  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
792  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
793  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
794  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
795  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
796  int q = p*_MODEL_NVARS_;
797 
798  /* Modified Rusanov's upwinding scheme */
799 
800  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
801  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
802  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
803  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
804  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
805 
806  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
807  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
808  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
809 
810  double c, vel[_MODEL_NDIMS_], rho,E,P;
811 
812  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
813  c = sqrt(param->gamma*P/rho);
814  double alphaL = c + absolute(vel[dir]);
815 
816  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
817  c = sqrt(param->gamma*P/rho);
818  double alphaR = c + absolute(vel[dir]);
819 
820  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
821  c = sqrt(param->gamma*P/rho);
822  double alphaavg = c + absolute(vel[dir]);
823 
824  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
825  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
826 
827  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
828  if (dir == _XDIR_) {
829  D[6] = alpha;
830  D[24] = alpha;
831  } else if (dir == _YDIR_) {
832  D[12] = alpha;
833  D[24] = alpha;
834  } else if (dir == _ZDIR_) {
835  D[18] = alpha;
836  D[24] = alpha;
837  }
838  MatMult5 (_MODEL_NVARS_,DL,D,L);
839  MatMult5 (_MODEL_NVARS_,modA,R,DL);
840  MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
841 
842  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
843  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
844  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
845  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
846  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
847  }
848  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
849  }
850 
851  return(0);
852 }
853 
862 int NavierStokes3DUpwindFdFRusanovModified(
863  double *fI,
864  double *fL,
865  double *fR,
866  double *uL,
867  double *uR,
868  double *u,
869  int dir,
870  void *s,
871  double t
872  )
873 {
874  HyPar *solver = (HyPar*) s;
875  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
876  int *dim = solver->dim_local, done;
877  double *uref = param->solution;
878 
879  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
880  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
881  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
882 
883  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
884  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
885  modA[_MODEL_NVARS_*_MODEL_NVARS_];
886 
887  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
888  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
889 
890  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
891 
892  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
893  while (!done) {
894  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
895  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
896 
897  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
898  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
899  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
900  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
901  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
902  int q = p*_MODEL_NVARS_;
903  double udiff[_MODEL_NVARS_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_],
904  udiss_total[_MODEL_NVARS_],udiss_acoustic[_MODEL_NVARS_];
905 
906  /* Modified Rusanov's upwinding scheme */
907  /* Modified Rusanov's upwinding scheme */
908  double c, vel[_MODEL_NDIMS_], rho,E,P, alphaL, alphaR, alphaavg,
909  betaL, betaR, betaavg, alpha, beta,
910  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
911 
912  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
913  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
914  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
915  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
916  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
917 
918  /* Compute total dissipation */
919  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
920  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
921  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
922 
923  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
924  c = sqrt(param->gamma*P/rho);
925  alphaL = c + absolute(vel[dir]);
926  betaL = absolute(vel[dir]);
927 
928  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
929  c = sqrt(param->gamma*P/rho);
930  alphaR = c + absolute(vel[dir]);
931  betaR = absolute(vel[dir]);
932 
933  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
934  c = sqrt(param->gamma*P/rho);
935  alphaavg = c + absolute(vel[dir]);
936  betaavg = absolute(vel[dir]);
937 
938  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
939  alpha = kappa*max3(alphaL,alphaR,alphaavg);
940  beta = kappa*max3(betaL,betaR,betaavg);
941 
942  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
943  if (dir == _XDIR_) {
944  D[0] = beta;
945  D[6] = alpha;
946  D[12] = beta;
947  D[18] = beta;
948  D[24] = alpha;
949  } else if (dir == _YDIR_) {
950  D[0] = beta;
951  D[6] = beta;
952  D[12] = alpha;
953  D[18] = beta;
954  D[24] = alpha;
955  } else if (dir == _ZDIR_) {
956  D[0] = beta;
957  D[6] = beta;
958  D[12] = beta;
959  D[18] = alpha;
960  D[24] = alpha;
961  }
962  MatMult5 (_MODEL_NVARS_,DL,D,L);
963  MatMult5 (_MODEL_NVARS_,modA,R,DL);
964  MatVecMult5 (_MODEL_NVARS_,udiss_total,modA,udiff);
965 
966  /* Compute dissipation for the linearized acoustic modes */
967  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
968  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
969  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
970 
971  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
972  c = sqrt(param->gamma*P/rho);
973  alphaL = c + absolute(vel[dir]);
974 
975  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
976  c = sqrt(param->gamma*P/rho);
977  alphaR = c + absolute(vel[dir]);
978 
979  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
980  c = sqrt(param->gamma*P/rho);
981  alphaavg = c + absolute(vel[dir]);
982 
983  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
984  alpha = kappa*max3(alphaL,alphaR,alphaavg);
985 
986  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
987  if (dir == _XDIR_) {
988  D[6] = alpha;
989  D[24] = alpha;
990  } else if (dir == _YDIR_) {
991  D[12] = alpha;
992  D[24] = alpha;
993  } else if (dir == _ZDIR_) {
994  D[18] = alpha;
995  D[24] = alpha;
996  }
997  MatMult5 (_MODEL_NVARS_,DL,D,L);
998  MatMult5 (_MODEL_NVARS_,modA,R,DL);
999  MatVecMult5 (_MODEL_NVARS_,udiss_acoustic,modA,udiff);
1000 
1001  /* Compute dissipation for the entropy modes */
1002  _ArraySubtract1D_(udiss,udiss_total,udiss_acoustic,_MODEL_NVARS_);
1003 
1004  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
1005  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
1006  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
1007  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
1008  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
1009  }
1010  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1011  }
1012 
1013  return(0);
1014 }
NavierStokes3DUpwindRoe
int NavierStokes3DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:40
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_NavierStokes3DGetFlowVar_
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
Definition: navierstokes3d.h:98
_ArrayIndex1D3_
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:69
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
max3
#define max3(a, b, c)
Definition: math_ops.h:27
navierstokes3d.h
3D Navier Stokes equations (compressible flows)
_NavierStokes3DLeftEigenvectors_
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
Definition: navierstokes3d.h:288
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
dummy
static const int dummy
Definition: NavierStokes3DUpwind.c:18
NavierStokes3DUpwindRusanov
int NavierStokes3DUpwindRusanov(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:349
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
NavierStokes3DUpwindFdFRusanovModified
int NavierStokes3DUpwindFdFRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:862
_ArrayCopy1D3_
#define _ArrayCopy1D3_(x, y, size)
Definition: arrayfunctions.h:343
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
MatVecMult5
#define MatVecMult5(N, y, A, x)
Definition: matmult_native.h:158
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
NavierStokes3DUpwinddFRoe
int NavierStokes3DUpwinddFRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:430
_NavierStokes3DEigenvalues_
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
Definition: navierstokes3d.h:254
NavierStokes3DUpwindRusanovModified
int NavierStokes3DUpwindRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:638
mathfunctions.h
Contains function definitions for common mathematical functions.
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
NavierStokes3DUpwindLLF
int NavierStokes3DUpwindLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:244
hypar.h
Contains structure definition for hypar.
MatMult5
#define MatMult5(N, A, X, Y)
Definition: matmult_native.h:82
NavierStokes3DUpwindRF
int NavierStokes3DUpwindRF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:140
NavierStokes3D::grav_field_g
double * grav_field_g
Definition: navierstokes3d.h:497
NavierStokes3D
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.
Definition: navierstokes3d.h:482
_NavierStokes3DRightEigenvectors_
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
Definition: navierstokes3d.h:384
basic.h
Some basic definitions and macros.
arrayfunctions.h
Contains macros and function definitions for common array operations.
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes3DRoeAverage_
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Definition: navierstokes3d.h:144
NavierStokes3DUpwinddFRusanovModified
int NavierStokes3DUpwinddFRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:756
max
#define max(a, b)
Definition: math_ops.h:18
NavierStokes3DUpwindFdFRoe
int NavierStokes3DUpwindFdFRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DUpwind.c:527
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
matmult_native.h
Contains macros and function definitions for common matrix multiplication.
_NavierStokes3D_stride_
static const int _NavierStokes3D_stride_
Definition: navierstokes3d.h:559