HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Functions
NavierStokes2DUpwind.c File Reference

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

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

Go to the source code of this file.

Functions

int NavierStokes2DUpwindRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindRF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindLLF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindSWFS (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindRusanov (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindRusanovModified (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwinddFRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwinddFRF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwinddFLLF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwinddFRusanovModified (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindFdFRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindFdFRF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindFdFLLF (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int NavierStokes2DUpwindFdFRusanovModified (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 

Detailed Description

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

Author
Debojyoti Ghosh

Definition in file NavierStokes2DUpwind.c.

Function Documentation

◆ NavierStokes2DUpwindRoe()

int NavierStokes2DUpwindRoe ( 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 2D 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 34 of file NavierStokes2DUpwind.c.

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 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwindRF()

int NavierStokes2DUpwindRF ( 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 120 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58

◆ NavierStokes2DUpwindLLF()

int NavierStokes2DUpwindLLF ( 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 2D 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, AIAA Journal, 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 223 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwindSWFS()

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

Steger-Warming Flux-Splitting scheme

  • Steger, J.L., Warming, R.F., "Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods", J. Comput. Phys., 40(2), 1981, pp. 263-293, http://dx.doi.org/10.1016/0021-9991(81)90210-2.

Note that this method cannot be used for 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 316 of file NavierStokes2DUpwind.c.

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 }
_NavierStokes2DGetFlowVar_
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Definition: navierstokes2d.h:81
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_ArrayAdd1D_
#define _ArrayAdd1D_(x, a, b, size)
Definition: arrayfunctions.h:212
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ NavierStokes2DUpwindRusanov()

int NavierStokes2DUpwindRusanov ( 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 2D 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, AIAA Journal, 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 428 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DGetFlowVar_
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Definition: navierstokes2d.h:81
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwindRusanovModified()

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

Modified Rusanov's upwinding scheme: NavierStokes2DUpwindRusanov() 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 499 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_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
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwinddFRoe()

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

The Roe upwinding scheme (NavierStokes2DUpwindRoe) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes2DStiffFlux, _NavierStokes2DSetStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 589 of file NavierStokes2DUpwind.c.

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 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwinddFRF()

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

The characteristic-based Roe-fixed upwinding scheme (NavierStokes2DUpwindRF) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes2DStiffFlux, _NavierStokes2DSetStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 670 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ NavierStokes2DUpwinddFLLF()

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

The characteristic-based local Lax-Friedrich upwinding scheme (NavierStokes2DUpwindLLF) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes2DStiffFlux, _NavierStokes2DSetStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 764 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwinddFRusanovModified()

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

The modified Rusanov upwinding scheme (NavierStokes2DUpwindRusanovModified()) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes2DStiffFlux, _NavierStokes2DSetStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 859 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_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
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwindFdFRoe()

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

The Roe upwinding scheme (NavierStokes2DUpwindRoe) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes2DNonStiffFlux, _NavierStokes2DSetStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 950 of file NavierStokes2DUpwind.c.

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 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwindFdFRF()

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

The characteristic-based Roe-fixed upwinding scheme (NavierStokes2DUpwindRF) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes2DNonStiffFlux, _NavierStokes2DSetNonStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 1046 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ NavierStokes2DUpwindFdFLLF()

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

The characteristic-based local Lax-Friedrich upwinding scheme (NavierStokes2DUpwindLLF) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes2DNonStiffFlux, _NavierStokes2DSetNonStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 1140 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_NavierStokes2DLeftEigenvectors_
#define _NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
Definition: navierstokes2d.h:247
_NavierStokes2DEigenvalues_
#define _NavierStokes2DEigenvalues_(u, D, gamma, dir)
Definition: navierstokes2d.h:213
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
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::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes2DUpwindFdFRusanovModified()

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

The modified Rusanov upwinding scheme (NavierStokes2DUpwindRusanovModified()) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes2DNonStiffFlux, _NavierStokes2DSetNonStiffFlux_). 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 or y)
sSolver object of type HyPar
tCurrent solution time

Definition at line 1235 of file NavierStokes2DUpwind.c.

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_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
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_NavierStokes2DRoeAverage_
#define _NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
Definition: navierstokes2d.h:171
_NavierStokes2DRightEigenvectors_
#define _NavierStokes2DRightEigenvectors_(u, R, ga, dir)
Definition: navierstokes2d.h:309
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
MatVecMult4
#define MatVecMult4(N, y, A, x)
Definition: matmult_native.h:147
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
MatMult4
#define MatMult4(N, A, X, Y)
Definition: matmult_native.h:59
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
_ArrayIndex1D2_
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:56
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18