HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NavierStokes3DUpwind.c File Reference

Contains functions to compute the upwind flux at grid interfaces for the 3D Navier Stokes equations. More...

#include <stdlib.h>
#include <math.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <matmult_native.h>
#include <physicalmodels/navierstokes3d.h>
#include <hypar.h>

Go to the source code of this file.

Functions

int NavierStokes3DUpwindRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwindRF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwindLLF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwindRusanov (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwinddFRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwindFdFRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwindRusanovModified (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwinddFRusanovModified (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes3DUpwindFdFRusanovModified (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 

Variables

static const int dummy = 1
 

Detailed Description

Contains functions to compute the upwind flux at grid interfaces for the 3D Navier Stokes equations.

Author
Debojyoti Ghosh

Definition in file NavierStokes3DUpwind.c.

Function Documentation

int NavierStokes3DUpwindRoe ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Roe's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \left| A\left({\bf u}_{j+1/2}^L,{\bf u}_{j+1/2}^R\right) \right| \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

This upwinding scheme is modified for the balanced discretization of the 3D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 40 of file NavierStokes3DUpwind.c.

51 {
52  HyPar *solver = (HyPar*) s;
53  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
54  int done;
55 
56  int *dim = solver->dim_local;
57 
58  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
59  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
60  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
64 
65  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
66 
67 #if defined(CPU_STAT)
68  clock_t startEvent, stopEvent;
69  startEvent = clock();
70 #endif
71 
72  while (!done) {
73  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
74  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
75  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
76  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
77  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
78  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
79  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
80  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
81 
82  /* Roe's upwinding scheme */
83 
84  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
85  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
86  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
87  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
88  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
89 
91  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
92  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
93  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
94 
95  /* Harten's Entropy Fix - Page 362 of Leveque */
96  int k;
97  double delta = 0.000001, delta2 = delta*delta;
98  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
99  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
100  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
101  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
102  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
103 
104  MatMult5(_MODEL_NVARS_,DL,D,L);
105  MatMult5(_MODEL_NVARS_,modA,R,DL);
106  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
107 
108  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
109  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
110  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
111  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
112  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
113  }
114  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
115  }
116 
117 #if defined(CPU_STAT)
118  stopEvent = clock();
119  printf("%-50s CPU time (secs) = %.6f dir = %d\n",
120  "NavierStokes3DUpwindRoe", (double)(stopEvent-startEvent)/CLOCKS_PER_SEC, dir);
121 #endif
122 
123  return(0);
124 }
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _ArrayCopy1D3_(x, y, size)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define MatVecMult5(N, y, A, x)
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D_(x, y, size)
#define MatMult5(N, A, X, Y)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwindRF ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Characteristic-based Roe-fixed upwinding scheme.

\begin{align} \alpha_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,L}, \\ \alpha_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,R}, \\ v_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,L}, \\ v_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,R}, \\ \alpha_{j+1/2}^k &= \left\{ \begin{array}{cc} \alpha_{j+1/2}^{k,L} & {\rm if}\ \lambda_{j,j+1/2,j+1} > 0 \\ \alpha_{j+1/2}^{k,R} & {\rm if}\ \lambda_{j,j+1/2,j+1} < 0 \\ \frac{1}{2}\left[ \alpha_{j+1/2}^{k,L} + \alpha_{j+1/2}^{k,R} - \left(\max_{\left[j,j+1\right]} \lambda\right) \left( v_{j+1/2}^{k,R} - v_{j+1/2}^{k,L} \right) \right] & {\rm otherwise} \end{array}\right., \\ {\bf f}_{j+1/2} &= \sum_{k=1}^3 \alpha_{j+1/2}^k {\bf r}_{j+1/2}^k \end{align}

where \({\bf l}\), \({\bf r}\), and \(\lambda\) are the left-eigenvectors, right-eigenvectors and eigenvalues. The subscripts denote the grid locations.

Note that this upwinding scheme cannot be used for solving flows with non-zero gravitational forces.

Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 140 of file NavierStokes3DUpwind.c.

151 {
152  HyPar *solver = (HyPar*) s;
153  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
154  int done,k;
155 
156  int *dim = solver->dim_local;
157 
158  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
159  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
160  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
162 
163  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
164  while (!done) {
165  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
166  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
167  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
168  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
169  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
170 
171  /* Roe-Fixed upwinding scheme */
172 
174 
175  _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
176  _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
177 
178  /* calculate characteristic fluxes and variables */
183 
184  double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
186  eigL[0] = D[0];
187  eigL[1] = D[6];
188  eigL[2] = D[12];
189  eigL[3] = D[18];
190  eigL[4] = D[24];
192  eigR[0] = D[0];
193  eigR[1] = D[6];
194  eigR[2] = D[12];
195  eigR[3] = D[18];
196  eigR[4] = D[24];
197  _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
198  eigC[0] = D[0];
199  eigC[1] = D[6];
200  eigC[2] = D[12];
201  eigC[3] = D[18];
202  eigC[4] = D[24];
203 
204  for (k = 0; k < _MODEL_NVARS_; k++) {
205  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
206  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
207  else {
208  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
209  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
210  }
211  }
212 
213  /* calculate the interface flux from the characteristic flux */
214  MatVecMult5(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
215  }
216  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
217  }
218 
219  return(0);
220 }
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
#define max3(a, b, c)
Definition: math_ops.h:27
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _ArrayCopy1D3_(x, y, size)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define MatVecMult5(N, y, A, x)
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwindLLF ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Characteristic-based local Lax-Friedrich upwinding scheme.

\begin{align} \alpha_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,L}, \\ \alpha_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,R}, \\ v_{j+1/2}^{k,L} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,L}, \\ v_{j+1/2}^{k,R} &= \sum_{k=1}^3 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,R}, \\ \alpha_{j+1/2}^k &= \frac{1}{2}\left[ \alpha_{j+1/2}^{k,L} + \alpha_{j+1/2}^{k,R} - \left(\max_{\left[j,j+1\right]} \lambda\right) \left( v_{j+1/2}^{k,R} - v_{j+1/2}^{k,L} \right) \right], \\ {\bf f}_{j+1/2} &= \sum_{k=1}^3 \alpha_{j+1/2}^k {\bf r}_{j+1/2}^k \end{align}

where \({\bf l}\), \({\bf r}\), and \(\lambda\) are the left-eigenvectors, right-eigenvectors and eigenvalues. The subscripts denote the grid locations.

This upwinding scheme is modified for the balanced discretization of the 3D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 244 of file NavierStokes3DUpwind.c.

255 {
256  HyPar *solver = (HyPar*) s;
257  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
258  int done;
259 
260  int *dim = solver->dim_local;
261 
262  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
263  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
264  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
266 
267  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
268  while (!done) {
269  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
270  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
271  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
272  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
273  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
274 
275  /* Local Lax-Friedrich upwinding scheme */
276 
278 
279  _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
280  _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
281 
282  /* calculate characteristic fluxes and variables */
287 
288  double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
290  eigL[0] = D[0];
291  eigL[1] = D[6];
292  eigL[2] = D[12];
293  eigL[3] = D[18];
294  eigL[4] = D[24];
296  eigR[0] = D[0];
297  eigR[1] = D[6];
298  eigR[2] = D[12];
299  eigR[3] = D[18];
300  eigR[4] = D[24];
301  _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
302  eigC[0] = D[0];
303  eigC[1] = D[6];
304  eigC[2] = D[12];
305  eigC[3] = D[18];
306  eigC[4] = D[24];
307 
308  double alpha;
309  alpha = max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
310  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
311  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
312  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
313  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
314  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
315  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
316  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
317  alpha = max3(absolute(eigL[4]),absolute(eigC[4]),absolute(eigR[4]));
318  fc[4] = 0.5 * (fcL[4] + fcR[4] + alpha * (ucL[4]-ucR[4]));
319 
320  /* calculate the interface flux from the characteristic flux */
322  }
323  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
324  }
325 
326  return(0);
327 }
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
#define max3(a, b, c)
Definition: math_ops.h:27
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _ArrayCopy1D3_(x, y, size)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define MatVecMult5(N, y, A, x)
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwindRusanov ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Rusanov's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \max_{j,j+1} \nu_j \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

where \(\nu = c + \left|u\right|\).

  • Rusanov, V. V., "The calculation of the interaction of non-stationary shock waves and obstacles," USSR Computational Mathematics and Mathematical Physics, Vol. 1, No. 2, 1962, pp. 304–320

This upwinding scheme is modified for the balanced discretization of the 3D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 349 of file NavierStokes3DUpwind.c.

360 {
361  HyPar *solver = (HyPar*) s;
362  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
363  int *dim = solver->dim_local, done;
364 
365  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
366  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
367  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
368 
369  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
370 
371  static int index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
372  indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_];
373 
374  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
375  while (!done) {
376  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
377  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
378 
379  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
380  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
381  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
382  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
383  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
384  int q = p*_MODEL_NVARS_;
385 
386  /* Modified Rusanov's upwinding scheme */
387 
388  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
389  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
390  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
391  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
392  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
393 
394  _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
395 
396  double c, vel[_MODEL_NDIMS_], rho,E,P;
397  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
398  c = sqrt(param->gamma*P/rho);
399  double alphaL = c + absolute(vel[dir]);
400  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
401  c = sqrt(param->gamma*P/rho);
402  double alphaR = c + absolute(vel[dir]);
403  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
404  c = sqrt(param->gamma*P/rho);
405  double alphaavg = c + absolute(vel[dir]);
406 
407  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
408  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
409 
410  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-alpha*udiff[0];
411  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-alpha*udiff[1];
412  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-alpha*udiff[2];
413  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-alpha*udiff[3];
414  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-alpha*udiff[4];
415  }
416  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
417  }
418 
419  return(0);
420 }
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
#define _ArraySetValue_(x, size, value)
#define max3(a, b, c)
Definition: math_ops.h:27
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
double * grav_field_g
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
#define max(a, b)
Definition: math_ops.h:18
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwinddFRoe ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

The Roe upwinding scheme (NavierStokes3DUpwindRoe) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes3DStiffFlux, _NavierStokes3DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \( u\pm a\) are used. Reference:

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 430 of file NavierStokes3DUpwind.c.

441 {
442  HyPar *solver = (HyPar*) s;
443  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
444  int done;
445 
446  int *dim = solver->dim_local;
447  double *uref = param->solution;
448 
449  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
450  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
451  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
455 
456  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
457  while (!done) {
458  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
459  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
460  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
461  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
462  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
463  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
464  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
465  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
466 
467  /* Roe's upwinding scheme */
468 
469  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
470  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
471  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
472  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
473  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
474 
476  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
477  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
478  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
479 
480  /* Harten's Entropy Fix - Page 362 of Leveque */
481  int k;
482  double delta = 0.000001, delta2 = delta*delta;
483  if (dir == _XDIR_) {
484  k=0; D[k] = 0.0;
485  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
486  k=12; D[k] = 0.0;
487  k=18; D[k] = 0.0;
488  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
489  } else if (dir == _YDIR_) {
490  k=0; D[k] = 0.0;
491  k=6; D[k] = 0.0;
492  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
493  k=18; D[k] = 0.0;
494  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
495  } else if (dir == _ZDIR_) {
496  k=0; D[k] = 0.0;
497  k=6; D[k] = 0.0;
498  k=12; D[k] = 0.0;
499  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
500  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
501  }
502 
503  MatMult5(_MODEL_NVARS_,DL,D,L);
504  MatMult5(_MODEL_NVARS_,modA,R,DL);
505  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
506 
507  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
508  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
509  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
510  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
511  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
512  }
513  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
514  }
515 
516  return(0);
517 }
#define _YDIR_
Definition: euler2d.h:41
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _ArrayCopy1D3_(x, y, size)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define MatVecMult5(N, y, A, x)
#define _ZDIR_
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D_(x, y, size)
#define MatMult5(N, A, X, Y)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwindFdFRoe ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

The Roe upwinding scheme (NavierStokes3DUpwindRoe) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes3DNonStiffFlux, _NavierStokes3DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \(u\) are used. Reference:

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 527 of file NavierStokes3DUpwind.c.

538 {
539  HyPar *solver = (HyPar*) s;
540  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
541  int done;
542 
543  int *dim = solver->dim_local;
544  double *uref = param->solution;
545 
546  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
547  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
548  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
552 
553  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
554  while (!done) {
555  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
556  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
557  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
558  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
559  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
560  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
561  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
562  int k;
563  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_],
564  udiss_total[_MODEL_NVARS_],udiss_stiff[_MODEL_NVARS_];
565  double delta = 0.000001, delta2 = delta*delta;
566 
567  /* Roe's upwinding scheme */
568 
569  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
570  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
571  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
572  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
573  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
574 
575  /* Compute total dissipation */
577  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
578  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
579  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
580  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
581  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
582  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
583  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
584  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
585  MatMult5(_MODEL_NVARS_,DL,D,L);
586  MatMult5(_MODEL_NVARS_,modA,R,DL);
587  MatVecMult5(_MODEL_NVARS_,udiss_total,modA,udiff);
588 
589  /* Compute dissipation corresponding to acoustic modes */
591  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
592  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
593  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
594  if (dir == _XDIR_) {
595  k=0; D[k] = 0.0;
596  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
597  k=12; D[k] = 0.0;
598  k=18; D[k] = 0.0;
599  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
600  } else if (dir == _YDIR_) {
601  k=0; D[k] = 0.0;
602  k=6; D[k] = 0.0;
603  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
604  k=18; D[k] = 0.0;
605  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
606  } else if (dir == _ZDIR_) {
607  k=0; D[k] = 0.0;
608  k=6; D[k] = 0.0;
609  k=12; D[k] = 0.0;
610  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
611  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
612  }
613  MatMult5(_MODEL_NVARS_,DL,D,L);
614  MatMult5(_MODEL_NVARS_,modA,R,DL);
615  MatVecMult5(_MODEL_NVARS_,udiss_stiff,modA,udiff);
616 
617  /* Compute the dissipation term for the entropy modes */
618  _ArraySubtract1D_(udiss,udiss_total,udiss_stiff,_MODEL_NVARS_);
619 
620  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
621  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
622  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
623  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
624  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
625  }
626  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
627  }
628 
629  return(0);
630 }
#define _YDIR_
Definition: euler2d.h:41
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _ArrayCopy1D3_(x, y, size)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArraySubtract1D_(x, a, b, size)
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define MatVecMult5(N, y, A, x)
#define _ZDIR_
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D_(x, y, size)
#define MatMult5(N, A, X, Y)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwindRusanovModified ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Modified Rusanov's upwinding scheme: NavierStokes3DUpwindRusanov() modified as described in the following paper (for consistent characteristic-based splitting):

  • Ghosh, D., Constantinescu, E. M., "Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning", Submitted (http://arxiv.org/abs/1510.05751).
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 638 of file NavierStokes3DUpwind.c.

649 {
650  HyPar *solver = (HyPar*) s;
651  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
652  int *dim = solver->dim_local, done;
653 
654  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
655  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
656  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
657 
658  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
659  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
660  modA[_MODEL_NVARS_*_MODEL_NVARS_];
661 
662  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
663  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
664 
665  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
666 
667  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
668  while (!done) {
669  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
670  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
671 
672  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
673  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
674  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
675  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
676  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
677  int q = p*_MODEL_NVARS_;
678 
679  /* Modified Rusanov's upwinding scheme */
680 
681  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
682  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
683  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
684  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
685  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
686 
687  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
688  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
689  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
690 
691  double c, vel[_MODEL_NDIMS_], rho,E,P;
692 
693  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
694  c = sqrt(param->gamma*P/rho);
695  double alphaL = c + absolute(vel[dir]);
696  double betaL = absolute(vel[dir]);
697 
698  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
699  c = sqrt(param->gamma*P/rho);
700  double alphaR = c + absolute(vel[dir]);
701  double betaR = absolute(vel[dir]);
702 
703  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
704  c = sqrt(param->gamma*P/rho);
705  double alphaavg = c + absolute(vel[dir]);
706  double betaavg = absolute(vel[dir]);
707 
708  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
709  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
710  double beta = kappa*max3(betaL,betaR,betaavg);
711 
712  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
713  if (dir == _XDIR_) {
714  D[0] = beta;
715  D[6] = alpha;
716  D[12] = beta;
717  D[18] = beta;
718  D[24] = alpha;
719  } else if (dir == _YDIR_) {
720  D[0] = beta;
721  D[6] = beta;
722  D[12] = alpha;
723  D[18] = beta;
724  D[24] = alpha;
725  } else if (dir == _ZDIR_) {
726  D[0] = beta;
727  D[6] = beta;
728  D[12] = beta;
729  D[18] = alpha;
730  D[24] = alpha;
731  }
732  MatMult5 (_MODEL_NVARS_,DL,D,L);
733  MatMult5 (_MODEL_NVARS_,modA,R,DL);
734  MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
735 
736  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
737  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
738  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
739  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
740  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
741  }
742  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
743  }
744 
745  return(0);
746 }
#define _YDIR_
Definition: euler2d.h:41
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
#define _ArraySetValue_(x, size, value)
#define max3(a, b, c)
Definition: math_ops.h:27
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define MatVecMult5(N, y, A, x)
#define _ZDIR_
#define _ArrayCopy1D_(x, y, size)
#define MatMult5(N, A, X, Y)
double * grav_field_g
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
#define max(a, b)
Definition: math_ops.h:18
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwinddFRusanovModified ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

The modified Rusanov upwinding scheme (NavierStokes3DUpwindRusanovModified()) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes3DStiffFlux, _NavierStokes3DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \( u\pm a\) are used. Reference:

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 756 of file NavierStokes3DUpwind.c.

767 {
768  HyPar *solver = (HyPar*) s;
769  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
770  int *dim = solver->dim_local, done;
771  double *uref = param->solution;
772 
773  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
774  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
775  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
776 
777  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
778  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
779  modA[_MODEL_NVARS_*_MODEL_NVARS_];
780 
781  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
782  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
783 
784  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
785 
786  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
787  while (!done) {
788  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
789  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
790 
791  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
792  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
793  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
794  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
795  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
796  int q = p*_MODEL_NVARS_;
797 
798  /* Modified Rusanov's upwinding scheme */
799 
800  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
801  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
802  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
803  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
804  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
805 
806  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
807  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
808  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
809 
810  double c, vel[_MODEL_NDIMS_], rho,E,P;
811 
812  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
813  c = sqrt(param->gamma*P/rho);
814  double alphaL = c + absolute(vel[dir]);
815 
816  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
817  c = sqrt(param->gamma*P/rho);
818  double alphaR = c + absolute(vel[dir]);
819 
820  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
821  c = sqrt(param->gamma*P/rho);
822  double alphaavg = c + absolute(vel[dir]);
823 
824  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
825  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
826 
827  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
828  if (dir == _XDIR_) {
829  D[6] = alpha;
830  D[24] = alpha;
831  } else if (dir == _YDIR_) {
832  D[12] = alpha;
833  D[24] = alpha;
834  } else if (dir == _ZDIR_) {
835  D[18] = alpha;
836  D[24] = alpha;
837  }
838  MatMult5 (_MODEL_NVARS_,DL,D,L);
839  MatMult5 (_MODEL_NVARS_,modA,R,DL);
840  MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
841 
842  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
843  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
844  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
845  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
846  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
847  }
848  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
849  }
850 
851  return(0);
852 }
#define _YDIR_
Definition: euler2d.h:41
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
#define _ArraySetValue_(x, size, value)
#define max3(a, b, c)
Definition: math_ops.h:27
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define MatVecMult5(N, y, A, x)
#define _ZDIR_
#define _ArrayCopy1D_(x, y, size)
#define MatMult5(N, A, X, Y)
double * grav_field_g
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
#define max(a, b)
Definition: math_ops.h:18
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75
static const int _NavierStokes3D_stride_
int NavierStokes3DUpwindFdFRusanovModified ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

The modified Rusanov upwinding scheme (NavierStokes3DUpwindRusanovModified()) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes3DNonStiffFlux, _NavierStokes3DSetNonStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \(u\) are used. Reference:

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 862 of file NavierStokes3DUpwind.c.

873 {
874  HyPar *solver = (HyPar*) s;
875  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
876  int *dim = solver->dim_local, done;
877  double *uref = param->solution;
878 
879  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
880  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
881  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
882 
883  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
884  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
885  modA[_MODEL_NVARS_*_MODEL_NVARS_];
886 
887  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
888  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
889 
890  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
891 
892  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
893  while (!done) {
894  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
895  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
896 
897  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
898  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
899  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
900  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
901  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
902  int q = p*_MODEL_NVARS_;
903  double udiff[_MODEL_NVARS_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_],
904  udiss_total[_MODEL_NVARS_],udiss_acoustic[_MODEL_NVARS_];
905 
906  /* Modified Rusanov's upwinding scheme */
907  /* Modified Rusanov's upwinding scheme */
908  double c, vel[_MODEL_NDIMS_], rho,E,P, alphaL, alphaR, alphaavg,
909  betaL, betaR, betaavg, alpha, beta,
910  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
911 
912  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
913  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
914  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
915  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
916  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
917 
918  /* Compute total dissipation */
919  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
920  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
921  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
922 
923  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
924  c = sqrt(param->gamma*P/rho);
925  alphaL = c + absolute(vel[dir]);
926  betaL = absolute(vel[dir]);
927 
928  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
929  c = sqrt(param->gamma*P/rho);
930  alphaR = c + absolute(vel[dir]);
931  betaR = absolute(vel[dir]);
932 
933  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
934  c = sqrt(param->gamma*P/rho);
935  alphaavg = c + absolute(vel[dir]);
936  betaavg = absolute(vel[dir]);
937 
938  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
939  alpha = kappa*max3(alphaL,alphaR,alphaavg);
940  beta = kappa*max3(betaL,betaR,betaavg);
941 
942  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
943  if (dir == _XDIR_) {
944  D[0] = beta;
945  D[6] = alpha;
946  D[12] = beta;
947  D[18] = beta;
948  D[24] = alpha;
949  } else if (dir == _YDIR_) {
950  D[0] = beta;
951  D[6] = beta;
952  D[12] = alpha;
953  D[18] = beta;
954  D[24] = alpha;
955  } else if (dir == _ZDIR_) {
956  D[0] = beta;
957  D[6] = beta;
958  D[12] = beta;
959  D[18] = alpha;
960  D[24] = alpha;
961  }
962  MatMult5 (_MODEL_NVARS_,DL,D,L);
963  MatMult5 (_MODEL_NVARS_,modA,R,DL);
964  MatVecMult5 (_MODEL_NVARS_,udiss_total,modA,udiff);
965 
966  /* Compute dissipation for the linearized acoustic modes */
967  _NavierStokes3DRoeAverage_ (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
968  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
969  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
970 
971  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
972  c = sqrt(param->gamma*P/rho);
973  alphaL = c + absolute(vel[dir]);
974 
975  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
976  c = sqrt(param->gamma*P/rho);
977  alphaR = c + absolute(vel[dir]);
978 
979  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
980  c = sqrt(param->gamma*P/rho);
981  alphaavg = c + absolute(vel[dir]);
982 
983  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
984  alpha = kappa*max3(alphaL,alphaR,alphaavg);
985 
986  _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
987  if (dir == _XDIR_) {
988  D[6] = alpha;
989  D[24] = alpha;
990  } else if (dir == _YDIR_) {
991  D[12] = alpha;
992  D[24] = alpha;
993  } else if (dir == _ZDIR_) {
994  D[18] = alpha;
995  D[24] = alpha;
996  }
997  MatMult5 (_MODEL_NVARS_,DL,D,L);
998  MatMult5 (_MODEL_NVARS_,modA,R,DL);
999  MatVecMult5 (_MODEL_NVARS_,udiss_acoustic,modA,udiff);
1000 
1001  /* Compute dissipation for the entropy modes */
1002  _ArraySubtract1D_(udiss,udiss_total,udiss_acoustic,_MODEL_NVARS_);
1003 
1004  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
1005  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
1006  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
1007  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
1008  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
1009  }
1010  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1011  }
1012 
1013  return(0);
1014 }
#define _YDIR_
Definition: euler2d.h:41
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
#define _ArraySetValue_(x, size, value)
#define max3(a, b, c)
Definition: math_ops.h:27
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _MODEL_NVARS_
Definition: euler1d.h:58
static const int dummy
void * physics
Definition: hypar.h:266
int * dim_local
Definition: hypar.h:37
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArraySubtract1D_(x, a, b, size)
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define MatVecMult5(N, y, A, x)
#define _ZDIR_
#define _ArrayCopy1D_(x, y, size)
#define MatMult5(N, A, X, Y)
double * grav_field_g
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
#define max(a, b)
Definition: math_ops.h:18
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _XDIR_
Definition: euler1d.h:75
static const int _NavierStokes3D_stride_

Variable Documentation

const int dummy = 1
static

Definition at line 18 of file NavierStokes3DUpwind.c.