HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
NavierStokes2DUpwind.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 <physicalmodels/navierstokes2d.h>
10 #include <mathfunctions.h>
11 #include <matmult_native.h>
12 #include <hypar.h>
13 
34 int NavierStokes2DUpwindRoe(
35  double *fI,
36  double *fL,
37  double *fR,
38  double *uL,
39  double *uR,
40  double *u,
41  int dir,
42  void *s,
43  double t
44  )
45 {
46  HyPar *solver = (HyPar*) s;
47  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
48  int done;
49 
50  int *dim = solver->dim_local;
51 
52  int bounds_outer[2], bounds_inter[2];
53  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
54  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
55  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
56  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
57  modA[_MODEL_NVARS_*_MODEL_NVARS_];
58 
59  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
60  while (!done) {
61  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
62  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
63  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
64  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
65  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
66  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
67  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
68  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
69 
70  /* Roe's upwinding scheme */
71 
72  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
73  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
74  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
75  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
76 
77  _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
78  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
79  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
80  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
81 
82  /* Harten's Entropy Fix - Page 362 of Leveque */
83  int k;
84  double delta = 0.000001, delta2 = delta*delta;
85  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
86  k=0; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
87  k=5; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
88  k=10; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
89  k=15; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
90 
91  MatMult4(_MODEL_NVARS_,DL,D,L);
92  MatMult4(_MODEL_NVARS_,modA,R,DL);
93  MatVecMult4(_MODEL_NVARS_,udiss,modA,udiff);
94 
95  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
96  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
97  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
98  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
99  }
100  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
101  }
102 
103  return(0);
104 }
105 
120 int NavierStokes2DUpwindRF(
121  double *fI,
122  double *fL,
123  double *fR,
124  double *uL,
125  double *uR,
126  double *u,
127  int dir,
128  void *s,
129  double t
130  )
131 {
132  HyPar *solver = (HyPar*) s;
133  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
134  int done,k;
135 
136  int *dim = solver->dim_local;
137 
138  int bounds_outer[2], bounds_inter[2];
139  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
140  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
141  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
142  L[_MODEL_NVARS_*_MODEL_NVARS_];
143 
144  done = 0; int index_outer[2] = {0,0}, index_inter[2];
145  while (!done) {
146  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
147  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
148  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
149  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
150  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
151 
152  /* Roe-Fixed upwinding scheme */
153 
154  _NavierStokes2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
155 
156  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
157  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
158  _NavierStokes2DRightEigenvectors_(uavg,R,param->gamma,dir);
159 
160  /* calculate characteristic fluxes and variables */
161  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
162  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
163  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
164  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
165 
166  double eigL[4],eigC[4],eigR[4];
167  _NavierStokes2DEigenvalues_((uL+_MODEL_NVARS_*p),D,param->gamma,dir);
168  eigL[0] = D[0];
169  eigL[1] = D[5];
170  eigL[2] = D[10];
171  eigL[3] = D[15];
172  _NavierStokes2DEigenvalues_((uR+_MODEL_NVARS_*p),D,param->gamma,dir);
173  eigR[0] = D[0];
174  eigR[1] = D[5];
175  eigR[2] = D[10];
176  eigR[3] = D[15];
177  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
178  eigC[0] = D[0];
179  eigC[1] = D[5];
180  eigC[2] = D[10];
181  eigC[3] = D[15];
182 
183  for (k = 0; k < _MODEL_NVARS_; k++) {
184  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
185  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
186  else {
187  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
188  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
189  }
190  }
191 
192  /* calculate the interface flux from the characteristic flux */
193  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
194  }
195  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
196  }
197 
198  return(0);
199 }
200 
223 int NavierStokes2DUpwindLLF(
224  double *fI,
225  double *fL,
226  double *fR,
227  double *uL,
228  double *uR,
229  double *u,
230  int dir,
231  void *s,
232  double t
233  )
234 {
235  HyPar *solver = (HyPar*) s;
236  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
237  int done;
238 
239  int *dim = solver->dim_local;
240 
241  int bounds_outer[2], bounds_inter[2];
242  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
243  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
244  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
245  L[_MODEL_NVARS_*_MODEL_NVARS_];
246 
247  done = 0; int index_outer[2] = {0,0}, index_inter[2];
248  while (!done) {
249  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
250  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
251  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
252  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
253  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
254  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
255  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
256  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
257  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
258  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
259 
260  /* Local Lax-Friedrich upwinding scheme */
261 
262  _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
263  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
264  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
265  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
266 
267  /* calculate characteristic fluxes and variables */
268  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
269  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
270  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
271  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
272 
273  double eigL[4],eigC[4],eigR[4];
274  _NavierStokes2DEigenvalues_((u+_MODEL_NVARS_*pL),D,param->gamma,dir);
275  eigL[0] = D[0];
276  eigL[1] = D[5];
277  eigL[2] = D[10];
278  eigL[3] = D[15];
279  _NavierStokes2DEigenvalues_((u+_MODEL_NVARS_*pR),D,param->gamma,dir);
280  eigR[0] = D[0];
281  eigR[1] = D[5];
282  eigR[2] = D[10];
283  eigR[3] = D[15];
284  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
285  eigC[0] = D[0];
286  eigC[1] = D[5];
287  eigC[2] = D[10];
288  eigC[3] = D[15];
289 
290  double alpha;
291  alpha = kappa * max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
292  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
293  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
294  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
295  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
296  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
297  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
298  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
299 
300  /* calculate the interface flux from the characteristic flux */
301  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
302  }
303  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
304  }
305 
306  return(0);
307 }
308 
316 int NavierStokes2DUpwindSWFS(
317  double *fI,
318  double *fL,
319  double *fR,
320  double *uL,
321  double *uR,
322  double *u,
323  int dir,
324  void *s,
325  double t
326  )
327 {
328  HyPar *solver = (HyPar*) s;
329  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
330  int done,k;
331  _DECLARE_IERR_;
332 
333  int ndims = solver->ndims;
334  int *dim = solver->dim_local;
335 
336  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
337  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
338  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
339  static double fp[_MODEL_NVARS_], fm[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
340 
341  done = 0; _ArraySetValue_(index_outer,ndims,0);
342  while (!done) {
343  _ArrayCopy1D_(index_outer,index_inter,ndims);
344  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
345  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
346  double rho,vx,vy,e,P,c,gamma=param->gamma,term,Mach,lp[_MODEL_NVARS_],lm[_MODEL_NVARS_];
347 
348  /* Steger Warming flux splitting */
349  _NavierStokes2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
350  _NavierStokes2DGetFlowVar_(uavg,rho,vx,vy,e,P,param->gamma);
351  Mach = (dir==_XDIR_ ? vx : vy) / sqrt(gamma*P/rho);
352 
353  if (Mach < -1.0) {
354 
355  _ArrayCopy1D_((fR+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);
356 
357  } else if (Mach < 1.0) {
358 
359  double kx = 0, ky = 0;
360  kx = (dir==_XDIR_ ? 1.0 : 0.0);
361  ky = (dir==_YDIR_ ? 1.0 : 0.0);
362 
363  _NavierStokes2DGetFlowVar_((uL+_MODEL_NVARS_*p),rho,vx,vy,e,P,param->gamma);
364  c = sqrt(gamma*P/rho);
365  term = rho/(2.0*gamma);
366  lp[0] = lp[1] = kx*vx + ky*vy;
367  lp[2] = lp[0] + c;
368  lp[3] = lp[0] - c;
369  for (k=0; k<_MODEL_NVARS_; k++) if (lp[k] < 0.0) lp[k] = 0.0;
370 
371  fp[0] = term * (2.0*(gamma-1.0)*lp[0] + lp[2] + lp[3]);
372  fp[1] = term * (2.0*(gamma-1.0)*lp[0]*vx + lp[2]*(vx+c*kx) + lp[3]*(vx-c*kx));
373  fp[2] = term * (2.0*(gamma-1.0)*lp[0]*vy + lp[2]*(vy+c*ky) + lp[3]*(vy-c*ky));
374  fp[3] = term * ((gamma-1.0)*lp[0]*(vx*vx+vy*vy) + 0.5*lp[2]*((vx+c*kx)*(vx+c*kx) + (vy+c*ky)*(vy+c*ky))
375  + 0.5*lp[3]*((vx-c*kx)*(vx-c*kx) + (vy-c*ky)*(vy-c*ky))
376  + ((3.0-gamma)*(lp[2]+lp[3])*c*c)/(2.0*(gamma-1.0)) );
377 
378  _NavierStokes2DGetFlowVar_((uR+_MODEL_NVARS_*p),rho,vx,vy,e,P,param->gamma);
379  c = sqrt(gamma*P/rho);
380  term = rho/(2.0*gamma);
381  lm[0] = lm[1] = kx*vx + ky*vy;
382  lm[2] = lm[0] + c;
383  lm[3] = lm[0] - c;
384  for (k=0; k<_MODEL_NVARS_; k++) if (lm[k] > 0.0) lm[k] = 0.0;
385 
386  fm[0] = term * (2.0*(gamma-1.0)*lm[0] + lm[2] + lm[3]);
387  fm[1] = term * (2.0*(gamma-1.0)*lm[0]*vx + lm[2]*(vx+c*kx) + lm[3]*(vx-c*kx));
388  fm[2] = term * (2.0*(gamma-1.0)*lm[0]*vy + lm[2]*(vy+c*ky) + lm[3]*(vy-c*ky));
389  fm[3] = term * ((gamma-1.0)*lm[0]*(vx*vx+vy*vy) + 0.5*lm[2]*((vx+c*kx)*(vx+c*kx) + (vy+c*ky)*(vy+c*ky))
390  + 0.5*lm[3]*((vx-c*kx)*(vx-c*kx) + (vy-c*ky)*(vy-c*ky))
391  + ((3.0-gamma)*(lm[2]+lm[3])*c*c)/(2.0*(gamma-1.0)) );
392 
393  _ArrayAdd1D_((fI+_MODEL_NVARS_*p),fp,fm,_MODEL_NVARS_);
394 
395  } else {
396 
397  _ArrayCopy1D_((fL+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);
398 
399  }
400 
401  }
402  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
403  }
404 
405  return(0);
406 }
407 
428 int NavierStokes2DUpwindRusanov(
429  double *fI,
430  double *fL,
431  double *fR,
432  double *uL,
433  double *uR,
434  double *u,
435  int dir,
436  void *s,
437  double t
438  )
439 {
440  HyPar *solver = (HyPar*) s;
441  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
442  int *dim = solver->dim_local, done;
443 
444  int bounds_outer[2], bounds_inter[2];
445  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
446  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
447 
448  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
449  while (!done) {
450  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
451  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
452  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
453  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
454  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
455  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
456  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
457  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
458 
459  /* Modified Rusanov's upwinding scheme */
460 
461  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
462  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
463  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
464  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
465 
466  _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
467 
468  double c, vel[_MODEL_NDIMS_], rho,E,P;
469  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
470  c = sqrt(param->gamma*P/rho);
471  double alphaL = c + absolute(vel[dir]), betaL = absolute(vel[dir]);
472  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
473  c = sqrt(param->gamma*P/rho);
474  double alphaR = c + absolute(vel[dir]), betaR = absolute(vel[dir]);
475  _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
476  c = sqrt(param->gamma*P/rho);
477  double alphaavg = c + absolute(vel[dir]), betaavg = absolute(vel[dir]);
478 
479  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
480  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
481 
482  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
483  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
484  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
485  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
486  }
487  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
488  }
489 
490  return(0);
491 }
492 
499 int NavierStokes2DUpwindRusanovModified(
500  double *fI,
501  double *fL,
502  double *fR,
503  double *uL,
504  double *uR,
505  double *u,
506  int dir,
507  void *s,
508  double t
509  )
510 {
511  HyPar *solver = (HyPar*) s;
512  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
513  int done;
514 
515  int *dim = solver->dim_local;
516 
517  int bounds_outer[2], bounds_inter[2];
518  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
519  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
520  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
521  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
522  modA[_MODEL_NVARS_*_MODEL_NVARS_];
523 
524  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
525  while (!done) {
526  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
527  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
528  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
529  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
530  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
531  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
532  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
533  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
534 
535  /* Modified Rusanov's upwinding scheme */
536 
537  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
538  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
539  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
540  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
541 
542  _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
543  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
544  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
545 
546  double c, vel[_MODEL_NDIMS_], rho,E,P;
547  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
548  c = sqrt(param->gamma*P/rho);
549  double alphaL = c + absolute(vel[dir]), betaL = absolute(vel[dir]);
550  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
551  c = sqrt(param->gamma*P/rho);
552  double alphaR = c + absolute(vel[dir]), betaR = absolute(vel[dir]);
553  _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
554  c = sqrt(param->gamma*P/rho);
555  double alphaavg = c + absolute(vel[dir]), betaavg = absolute(vel[dir]);
556 
557  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
558  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
559  double beta = kappa*max3(betaL,betaR,betaavg);
560 
561  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
562  D[0] = alpha;
563  D[5] = (dir == _XDIR_ ? alpha : beta);
564  D[10] = (dir == _YDIR_ ? alpha : beta);
565  D[15] = beta;
566  MatMult4(_MODEL_NVARS_,DL,D,L);
567  MatMult4(_MODEL_NVARS_,modA,R,DL);
568  MatVecMult4(_MODEL_NVARS_,udiss,modA,udiff);
569 
570  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-udiss[0];
571  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-udiss[1];
572  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-udiss[2];
573  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-udiss[3];
574  }
575  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
576  }
577 
578  return(0);
579 }
580 
589 int NavierStokes2DUpwinddFRoe(
590  double *fI,
591  double *fL,
592  double *fR,
593  double *uL,
594  double *uR,
595  double *u,
596  int dir,
597  void *s,
598  double t
599  )
600 {
601  HyPar *solver = (HyPar*) s;
602  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
603  int done;
604 
605  int *dim = solver->dim_local;
606  double *uref = param->solution;
607 
608  int bounds_outer[2], bounds_inter[2];
609  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
610  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
611  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
612  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
613  modA[_MODEL_NVARS_*_MODEL_NVARS_];
614 
615  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
616  while (!done) {
617  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
618  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
619  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
620  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
621  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
622  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
623  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
624  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
625 
626  /* Roe's upwinding scheme */
627 
628  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
629  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
630  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
631  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
632 
633  _NavierStokes2DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
634  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
635  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
636  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
637 
638  /* Harten's Entropy Fix - Page 362 of Leveque */
639  int k;
640  double delta = 0.000001, delta2 = delta*delta;
641  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
642  k=0; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
643  k=5; D[k] = (dir == _YDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
644  k=10; D[k] = (dir == _XDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
645  k=15; D[k] = 0.0;
646 
647  MatMult4(_MODEL_NVARS_,DL,D,L);
648  MatMult4(_MODEL_NVARS_,modA,R,DL);
649  MatVecMult4(_MODEL_NVARS_,udiss,modA,udiff);
650 
651  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
652  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
653  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
654  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
655  }
656  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
657  }
658 
659  return(0);
660 }
661 
670 int NavierStokes2DUpwinddFRF(
671  double *fI,
672  double *fL,
673  double *fR,
674  double *uL,
675  double *uR,
676  double *u,
677  int dir,
678  void *s,
679  double t
680  )
681 {
682  HyPar *solver = (HyPar*) s;
683  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
684  int done,k;
685 
686  int *dim = solver->dim_local;
687  double *uref = param->solution;
688 
689  int bounds_outer[2], bounds_inter[2];
690  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
691  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
692  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
693  L[_MODEL_NVARS_*_MODEL_NVARS_];
694 
695  done = 0; int index_outer[2] = {0,0}, index_inter[2];
696  while (!done) {
697  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
698  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
699  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
700  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
701  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
702  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
703  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
704  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
705  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
706 
707  /* Roe-Fixed upwinding scheme */
708 
709  _NavierStokes2DRoeAverage_(uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
710 
711  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
712  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
713  _NavierStokes2DRightEigenvectors_(uavg,R,param->gamma,dir);
714 
715  /* calculate characteristic fluxes and variables */
716  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
717  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
718  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
719  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
720 
721  double eigL[4],eigC[4],eigR[4];
722  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
723  eigL[0] = D[0];
724  eigL[1] = (dir == _YDIR_ ? 0.0 : D[5]);
725  eigL[2] = (dir == _XDIR_ ? 0.0 : D[10]);
726  eigL[3] = 0.0;
727  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
728  eigR[0] = D[0];
729  eigR[1] = (dir == _YDIR_ ? 0.0 : D[5]);
730  eigR[2] = (dir == _XDIR_ ? 0.0 : D[10]);
731  eigR[3] = 0.0;
732  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
733  eigC[0] = D[0];
734  eigC[1] = (dir == _YDIR_ ? 0.0 : D[5]);
735  eigC[2] = (dir == _XDIR_ ? 0.0 : D[10]);
736  eigC[3] = 0.0;
737 
738  for (k = 0; k < _MODEL_NVARS_; k++) {
739  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
740  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
741  else {
742  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
743  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
744  }
745  }
746 
747  /* calculate the interface flux from the characteristic flux */
748  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
749  }
750  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
751  }
752 
753  return(0);
754 }
755 
764 int NavierStokes2DUpwinddFLLF(
765  double *fI,
766  double *fL,
767  double *fR,
768  double *uL,
769  double *uR,
770  double *u,
771  int dir,
772  void *s,
773  double t
774  )
775 {
776  HyPar *solver = (HyPar*) s;
777  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
778  int done;
779 
780  int *dim = solver->dim_local;
781  double *uref = param->solution;
782 
783  int bounds_outer[2], bounds_inter[2];
784  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
785  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
786  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
787  L[_MODEL_NVARS_*_MODEL_NVARS_];
788 
789  done = 0; int index_outer[2] = {0,0}, index_inter[2];
790  while (!done) {
791  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
792  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
793  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
794  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
795  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
796  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
797  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
798  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
799  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
800  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
801 
802  /* Local Lax-Friedrich upwinding scheme */
803 
804  _NavierStokes2DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
805  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
806  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
807  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
808 
809  /* calculate characteristic fluxes and variables */
810  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
811  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
812  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
813  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
814 
815  double eigL[4],eigC[4],eigR[4];
816  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
817  eigL[0] = D[0];
818  eigL[1] = (dir == _YDIR_ ? 0.0 : D[5]);
819  eigL[2] = (dir == _XDIR_ ? 0.0 : D[10]);
820  eigL[3] = 0.0;
821  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
822  eigR[0] = D[0];
823  eigR[1] = (dir == _YDIR_ ? 0.0 : D[5]);
824  eigR[2] = (dir == _XDIR_ ? 0.0 : D[10]);
825  eigR[3] = 0.0;
826  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
827  eigC[0] = D[0];
828  eigC[1] = (dir == _YDIR_ ? 0.0 : D[5]);
829  eigC[2] = (dir == _XDIR_ ? 0.0 : D[10]);
830  eigC[3] = 0.0;
831 
832  double alpha;
833  alpha = kappa * max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
834  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
835  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
836  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
837  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
838  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
839  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
840  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
841 
842  /* calculate the interface flux from the characteristic flux */
843  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
844  }
845  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
846  }
847 
848  return(0);
849 }
850 
859 int NavierStokes2DUpwinddFRusanovModified(
860  double *fI,
861  double *fL,
862  double *fR,
863  double *uL,
864  double *uR,
865  double *u,
866  int dir,
867  void *s,
868  double t
869  )
870 {
871  HyPar *solver = (HyPar*) s;
872  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
873  int done;
874 
875  int *dim = solver->dim_local;
876  double *uref = param->solution;
877 
878  int bounds_outer[2], bounds_inter[2];
879  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
880  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
881  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
882  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
883  modA[_MODEL_NVARS_*_MODEL_NVARS_];
884 
885  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
886  while (!done) {
887  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
888  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
889  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
890  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
891  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
892  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
893  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
894  double udiff[_MODEL_NVARS_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_];
895 
896  /* Rusanov's upwinding scheme */
897 
898  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
899  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
900  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
901  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
902 
903  _NavierStokes2DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
904  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
905  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
906 
907  double c, vel[_MODEL_NDIMS_], rho,E,P;
908  _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
909  c = sqrt(param->gamma*P/rho);
910  double alphaL = c + absolute(vel[dir]);
911  _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
912  c = sqrt(param->gamma*P/rho);
913  double alphaR = c + absolute(vel[dir]);
914  _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
915  c = sqrt(param->gamma*P/rho);
916  double alphaavg = c + absolute(vel[dir]);
917 
918  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
919  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
920 
921  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
922  D[0] = alpha;
923  D[5] = (dir == _YDIR_ ? 0.0 : alpha);
924  D[10] = (dir == _XDIR_ ? 0.0 : alpha);
925  D[15] = 0.0;
926 
927  MatMult4(_MODEL_NVARS_,DL,D,L);
928  MatMult4(_MODEL_NVARS_,modA,R,DL);
929  MatVecMult4(_MODEL_NVARS_,udiss,modA,udiff);
930 
931  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
932  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
933  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
934  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
935  }
936  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
937  }
938 
939  return(0);
940 }
941 
950 int NavierStokes2DUpwindFdFRoe(
951  double *fI,
952  double *fL,
953  double *fR,
954  double *uL,
955  double *uR,
956  double *u,
957  int dir,
958  void *s,
959  double t
960  )
961 {
962  HyPar *solver = (HyPar*) s;
963  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
964  int done;
965 
966  int *dim = solver->dim_local;
967  double *uref = param->solution;
968 
969  int bounds_outer[2], bounds_inter[2];
970  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
971  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
972  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
973  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
974  modA[_MODEL_NVARS_*_MODEL_NVARS_];
975 
976  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
977  while (!done) {
978  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
979  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
980  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
981  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
982  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
983  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
984  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
985  int k;
986  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_],
987  udiss_total[_MODEL_NVARS_],udiss_stiff[_MODEL_NVARS_];
988  double delta = 0.000001, delta2 = delta*delta;
989  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
990 
991  /* Roe's upwinding scheme */
992 
993  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
994  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
995  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
996  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
997 
998  /* Compute total dissipation */
999  _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
1000  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
1001  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
1002  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
1003  k=0; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
1004  k=5; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
1005  k=10; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
1006  k=15; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
1007  MatMult4(_MODEL_NVARS_,DL,D,L);
1008  MatMult4(_MODEL_NVARS_,modA,R,DL);
1009  MatVecMult4(_MODEL_NVARS_,udiss_total,modA,udiff);
1010 
1011  /* Compute dissipation corresponding to acoustic modes */
1012  _NavierStokes2DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
1013  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
1014  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
1015  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
1016  k=0; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
1017  k=5; D[k] = (dir == _YDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
1018  k=10; D[k] = (dir == _XDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
1019  k=15; D[k] = 0.0;
1020  MatMult4(_MODEL_NVARS_,DL,D,L);
1021  MatMult4(_MODEL_NVARS_,modA,R,DL);
1022  MatVecMult4(_MODEL_NVARS_,udiss_stiff,modA,udiff);
1023 
1024  /* Compute the dissipation term for the entropy modes */
1025  _ArraySubtract1D_(udiss,udiss_total,udiss_stiff,_MODEL_NVARS_);
1026 
1027  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
1028  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
1029  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
1030  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
1031  }
1032  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1033  }
1034 
1035  return(0);
1036 }
1037 
1046 int NavierStokes2DUpwindFdFRF(
1047  double *fI,
1048  double *fL,
1049  double *fR,
1050  double *uL,
1051  double *uR,
1052  double *u,
1053  int dir,
1054  void *s,
1055  double t
1056  )
1057 {
1058  HyPar *solver = (HyPar*) s;
1059  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
1060  int done,k;
1061 
1062  int *dim = solver->dim_local;
1063  double *uref = param->solution;
1064 
1065  int bounds_outer[2], bounds_inter[2];
1066  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
1067  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
1068  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
1069  L[_MODEL_NVARS_*_MODEL_NVARS_];
1070 
1071  done = 0; int index_outer[2] = {0,0}, index_inter[2];
1072  while (!done) {
1073  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
1074  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
1075  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
1076  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
1077  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
1078  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
1079  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
1080  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
1081  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
1082 
1083  /* Roe-Fixed upwinding scheme */
1084 
1085  _NavierStokes2DRoeAverage_(uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
1086 
1087  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
1088  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
1089  _NavierStokes2DRightEigenvectors_(uavg,R,param->gamma,dir);
1090 
1091  /* calculate characteristic fluxes and variables */
1092  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
1093  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
1094  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
1095  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
1096 
1097  double eigL[4],eigC[4],eigR[4];
1098  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
1099  eigL[0] = 0.0;
1100  eigL[1] = (dir == _XDIR_ ? 0.0 : D[5]);
1101  eigL[2] = (dir == _YDIR_ ? 0.0 : D[10]);
1102  eigL[3] = D[15];
1103  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
1104  eigR[0] = 0.0;
1105  eigR[1] = (dir == _XDIR_ ? 0.0 : D[5]);
1106  eigR[2] = (dir == _YDIR_ ? 0.0 : D[10]);
1107  eigR[3] = D[15];
1108  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
1109  eigC[0] = 0.0;
1110  eigC[1] = (dir == _XDIR_ ? 0.0 : D[5]);
1111  eigC[2] = (dir == _YDIR_ ? 0.0 : D[10]);
1112  eigC[3] = D[15];
1113 
1114  for (k = 0; k < _MODEL_NVARS_; k++) {
1115  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
1116  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
1117  else {
1118  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
1119  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
1120  }
1121  }
1122 
1123  /* calculate the interface flux from the characteristic flux */
1124  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
1125  }
1126  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1127  }
1128 
1129  return(0);
1130 }
1131 
1140 int NavierStokes2DUpwindFdFLLF(
1141  double *fI,
1142  double *fL,
1143  double *fR,
1144  double *uL,
1145  double *uR,
1146  double *u,
1147  int dir,
1148  void *s,
1149  double t
1150  )
1151 {
1152  HyPar *solver = (HyPar*) s;
1153  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
1154  int done;
1155 
1156  int *dim = solver->dim_local;
1157  double *uref = param->solution;
1158 
1159  int bounds_outer[2], bounds_inter[2];
1160  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
1161  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
1162  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
1163  L[_MODEL_NVARS_*_MODEL_NVARS_];
1164 
1165  done = 0; int index_outer[2] = {0,0}, index_inter[2];
1166  while (!done) {
1167  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
1168  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
1169  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
1170  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
1171  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
1172  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
1173  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
1174  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
1175  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
1176  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
1177 
1178  /* Local Lax-Friedrich upwinding scheme */
1179 
1180  _NavierStokes2DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
1181  _NavierStokes2DEigenvalues_ (uavg,D,param->gamma,dir);
1182  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
1183  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
1184 
1185  /* calculate characteristic fluxes and variables */
1186  MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
1187  MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
1188  MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
1189  MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
1190 
1191  double eigL[4],eigC[4],eigR[4];
1192  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
1193  eigL[0] = 0.0;
1194  eigL[1] = (dir == _XDIR_ ? 0.0 : D[5]);
1195  eigL[2] = (dir == _YDIR_ ? 0.0 : D[10]);
1196  eigL[3] = D[15];
1197  _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
1198  eigR[0] = 0.0;
1199  eigR[1] = (dir == _XDIR_ ? 0.0 : D[5]);
1200  eigR[2] = (dir == _YDIR_ ? 0.0 : D[10]);
1201  eigR[3] = D[15];
1202  _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
1203  eigC[0] = 0.0;
1204  eigC[1] = (dir == _XDIR_ ? 0.0 : D[5]);
1205  eigC[2] = (dir == _YDIR_ ? 0.0 : D[10]);
1206  eigC[3] = D[15];
1207 
1208  double alpha;
1209  alpha = kappa * max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
1210  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
1211  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
1212  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
1213  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
1214  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
1215  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
1216  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
1217 
1218  /* calculate the interface flux from the characteristic flux */
1219  MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
1220  }
1221  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1222  }
1223 
1224  return(0);
1225 }
1226 
1235 int NavierStokes2DUpwindFdFRusanovModified(
1236  double *fI,
1237  double *fL,
1238  double *fR,
1239  double *uL,
1240  double *uR,
1241  double *u,
1242  int dir,
1243  void *s,
1244  double t
1245  )
1246 {
1247  HyPar *solver = (HyPar*) s;
1248  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
1249  int done;
1250 
1251  int *dim = solver->dim_local;
1252  double *uref = param->solution;
1253 
1254  int bounds_outer[2], bounds_inter[2];
1255  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
1256  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
1257  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
1258  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
1259  modA[_MODEL_NVARS_*_MODEL_NVARS_];
1260 
1261  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
1262  while (!done) {
1263  index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
1264  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
1265  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
1266  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
1267  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
1268  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
1269  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
1270  double udiff[_MODEL_NVARS_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_],
1271  udiss_total[_MODEL_NVARS_],udiss_acoustic[_MODEL_NVARS_];
1272 
1273  /* Modified Rusanov's upwinding scheme */
1274  double c, vel[_MODEL_NDIMS_], rho,E,P, alphaL, alphaR, alphaavg,
1275  betaL, betaR, betaavg, alpha, beta,
1276  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
1277 
1278 
1279  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
1280  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
1281  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
1282  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
1283 
1284  /* Compute total dissipation */
1285  _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
1286  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
1287  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
1288  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
1289  c = sqrt(param->gamma*P/rho);
1290  alphaL = c + absolute(vel[dir]);
1291  betaL = absolute(vel[dir]);
1292  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
1293  c = sqrt(param->gamma*P/rho);
1294  alphaR = c + absolute(vel[dir]);
1295  betaR = absolute(vel[dir]);
1296  _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
1297  c = sqrt(param->gamma*P/rho);
1298  alphaavg = c + absolute(vel[dir]);
1299  betaavg = absolute(vel[dir]);
1300  alpha = kappa*max3(alphaL,alphaR,alphaavg);
1301  beta = kappa*max3(betaL,betaR,betaavg);
1302  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
1303  D[0] = alpha;
1304  D[5] = (dir == _XDIR_ ? alpha : beta);
1305  D[10] = (dir == _YDIR_ ? alpha : beta);
1306  D[15] = beta;
1307  MatMult4(_MODEL_NVARS_,DL,D,L);
1308  MatMult4(_MODEL_NVARS_,modA,R,DL);
1309  MatVecMult4(_MODEL_NVARS_,udiss_total,modA,udiff);
1310 
1311  /* Compute dissipation for the linearized acoustic modes */
1312  _NavierStokes2DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
1313  _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
1314  _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
1315  _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
1316  c = sqrt(param->gamma*P/rho);
1317  alphaL = c + absolute(vel[dir]);
1318  _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
1319  c = sqrt(param->gamma*P/rho);
1320  alphaR = c + absolute(vel[dir]);
1321  _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
1322  c = sqrt(param->gamma*P/rho);
1323  alphaavg = c + absolute(vel[dir]);
1324  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
1325  alpha = kappa*max3(alphaL,alphaR,alphaavg);
1326  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
1327  D[0] = alpha;
1328  D[5] = (dir == _YDIR_ ? 0.0 : alpha);
1329  D[10] = (dir == _XDIR_ ? 0.0 : alpha);
1330  D[15] = 0.0;
1331  MatMult4(_MODEL_NVARS_,DL,D,L);
1332  MatMult4(_MODEL_NVARS_,modA,R,DL);
1333  MatVecMult4(_MODEL_NVARS_,udiss_acoustic,modA,udiff);
1334 
1335  /* Compute dissipation for the entropy modes */
1336  _ArraySubtract1D_(udiss,udiss_total,udiss_acoustic,_MODEL_NVARS_);
1337 
1338  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
1339  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
1340  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
1341  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
1342  }
1343  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1344  }
1345 
1346  return(0);
1347 }
NavierStokes2DUpwindRF
int NavierStokes2DUpwindRF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:120
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
mathfunctions.h
Contains function definitions for common mathematical functions.
_NavierStokes2DGetFlowVar_
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Definition: navierstokes2d.h:81
basic.h
Some basic definitions and macros.
NavierStokes2DUpwindFdFLLF
int NavierStokes2DUpwindFdFLLF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:1140
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
NavierStokes2DUpwindRoe
int NavierStokes2DUpwindRoe(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:34
HyPar::ndims
int ndims
Definition: hypar.h:26
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
NavierStokes2DUpwindLLF
int NavierStokes2DUpwindLLF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:223
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_ArrayAdd1D_
#define _ArrayAdd1D_(x, a, b, size)
Definition: arrayfunctions.h:212
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
hypar.h
Contains structure definition for hypar.
navierstokes2d.h
2D Navier Stokes equations (compressible flows)
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2DUpwinddFRoe
int NavierStokes2DUpwinddFRoe(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:589
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
NavierStokes2DUpwindRusanovModified
int NavierStokes2DUpwindRusanovModified(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:499
HyPar::physics
void * physics
Definition: hypar.h:266
NavierStokes2DUpwindFdFRoe
int NavierStokes2DUpwindFdFRoe(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:950
NavierStokes2DUpwinddFRusanovModified
int NavierStokes2DUpwinddFRusanovModified(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:859
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
NavierStokes2DUpwinddFLLF
int NavierStokes2DUpwinddFLLF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:764
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
NavierStokes2DUpwindRusanov
int NavierStokes2DUpwindRusanov(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:428
matmult_native.h
Contains macros and function definitions for common matrix multiplication.
NavierStokes2DUpwindFdFRusanovModified
int NavierStokes2DUpwindFdFRusanovModified(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:1235
NavierStokes2DUpwinddFRF
int NavierStokes2DUpwinddFRF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:670
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
NavierStokes2DUpwindSWFS
int NavierStokes2DUpwindSWFS(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:316
NavierStokes2DUpwindFdFRF
int NavierStokes2DUpwindFdFRF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: NavierStokes2DUpwind.c:1046
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18
arrayfunctions.h
Contains macros and function definitions for common array operations.