HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
NavierStokes2DUpwindRF
int NavierStokes2DUpwindRF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:120
max3
#define max3(a, b, c)
Definition: math_ops.h:27
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
NavierStokes2DUpwinddFRusanovModified
int NavierStokes2DUpwinddFRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:859
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
NavierStokes2DUpwinddFLLF
int NavierStokes2DUpwinddFLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:764
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2DUpwindRusanov
int NavierStokes2DUpwindRusanov(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:428
NavierStokes2DUpwindFdFLLF
int NavierStokes2DUpwindFdFLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:1140
navierstokes2d.h
2D Navier Stokes equations (compressible flows)
NavierStokes2DUpwindFdFRF
int NavierStokes2DUpwindFdFRF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:1046
NavierStokes2DUpwindFdFRoe
int NavierStokes2DUpwindFdFRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:950
mathfunctions.h
Contains function definitions for common mathematical functions.
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
NavierStokes2DUpwindRusanovModified
int NavierStokes2DUpwindRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:499
NavierStokes2DUpwinddFRoe
int NavierStokes2DUpwinddFRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:589
NavierStokes2DUpwindRoe
int NavierStokes2DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:34
hypar.h
Contains structure definition for hypar.
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_NavierStokes2DGetFlowVar_
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Definition: navierstokes2d.h:81
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
basic.h
Some basic definitions and macros.
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
HyPar::ndims
int ndims
Definition: hypar.h:26
NavierStokes2DUpwinddFRF
int NavierStokes2DUpwinddFRF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:670
arrayfunctions.h
Contains macros and function definitions for common array operations.
absolute
#define absolute(a)
Definition: math_ops.h:32
max
#define max(a, b)
Definition: math_ops.h:18
NavierStokes2DUpwindSWFS
int NavierStokes2DUpwindSWFS(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:316
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
NavierStokes2DUpwindFdFRusanovModified
int NavierStokes2DUpwindFdFRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:1235
_ArrayAdd1D_
#define _ArrayAdd1D_(x, a, b, size)
Definition: arrayfunctions.h:212
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
matmult_native.h
Contains macros and function definitions for common matrix multiplication.
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
NavierStokes2DUpwindLLF
int NavierStokes2DUpwindLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DUpwind.c:223