Initialization of the physics-related variables and function pointers for the 3D Navier-Stokes system. More...

#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <boundaryconditions.h>
#include <physicalmodels/navierstokes3d.h>
#include <mpivars.h>
#include <hypar.h>

Functions
double	NavierStokes3DComputeCFL (void , void , double, double)

int	NavierStokes3DFlux (double , double , int, void *, double)

int	NavierStokes3DStiffFlux (double , double , int, void *, double)

int	NavierStokes3DNonStiffFlux (double , double , int, void *, double)

int	NavierStokes3DRoeAverage (double , double , double , void )

int	NavierStokes3DParabolicFunction (double , double , void , void , double)

int	NavierStokes3DSource (double , double , void , void , double)

int	NavierStokes3DJacobian (double , double , void *, int, int, int)

int	NavierStokes3DStiffJacobian (double , double , void *, int, int, int)

int	NavierStokes3DLeftEigenvectors (double , double , void *, int)

int	NavierStokes3DRightEigenvectors (double , double , void *, int)

int	NavierStokes3DUpwindRoe (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwindRF (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwindLLF (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwindRusanov (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwinddFRoe (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwindFdFRoe (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwindRusanovModified (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwinddFRusanovModified (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DUpwindFdFRusanovModified (double , double , double , double , double , double , int, void *, double)

int	NavierStokes3DGravityField (void , void )

int	NavierStokes3DModifiedSolution (double , double , int, void , void , double)

int	NavierStokes3DIBAdiabatic (void , void , double *, double)

int	NavierStokes3DIBIsothermal (void , void , double *, double)

int	NavierStokes3DPreStep (double , void , void *, double)

int	NavierStokes3DIBForces (void , void , double)

int	gpuNavierStokes3DInitialize (void , void )

int	gpuNavierStokes3DFlux (double , double , int, void *, double)

int	gpuNavierStokes3DParabolicFunction (double , double , void , void , double)

int	gpuNavierStokes3DSource (double , double , void , void , double)

int	gpuNavierStokes3DModifiedSolution (double , double , int, void , void , double)

int	gpuNavierStokes3DUpwindRusanov (double , double , double , double , double , double , int, void *, double)

int	gpuNavierStokes3DUpwindRoe (double , double , double , double , double , double , int, void *, double)

int	gpuNavierStokes3DPreStep (double , void , void *, double)

int	NavierStokes3DInitialize (void s, void m)

Detailed Description

Initialization of the physics-related variables and function pointers for the 3D Navier-Stokes system.

Author: Debojyoti Ghosh

Definition in file NavierStokes3DInitialize.c.

Function Documentation

double NavierStokes3DComputeCFL	(	void *	s,
		void *	m,
		double	dt,
		double	t
	)

Computes the maximum CFL number over the domain. Note that the CFL is computed over the local domain on this processor only.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
dt	Time step size for which to compute the CFL
t	Time

Definition at line 16 of file NavierStokes3DComputeCFL.c.

 {
   HyPar             *solver = (HyPar*)   s;
   NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
 
   int *dim    = solver->dim_local;
   int ghosts  = solver->ghosts;
   int ndims   = solver->ndims;
   int index[ndims];
   double *u   = solver->u;
 
   double max_cfl = 0;
   int done = 0; _ArraySetValue_(index,ndims,0);
   while (!done) {
     int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
     double rho, vx, vy, vz, e, P, c, dxinv, dyinv, dzinv, local_cfl[3];
     _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,vx,vy,vz,e,P,param->gamma);
 
     c = sqrt(param->gamma*P/rho); /* speed of sound */
     _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
     _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv); /* 1/dy */
     _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv); /* 1/dz */
 
     local_cfl[_XDIR_] = (absolute(vx)+c)*dt*dxinv; /* local cfl for this grid point (x) */
     local_cfl[_YDIR_] = (absolute(vy)+c)*dt*dyinv; /* local cfl for this grid point (y) */
     local_cfl[_ZDIR_] = (absolute(vz)+c)*dt*dzinv; /* local cfl for this grid point (z) */
     if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
     if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];
     if (local_cfl[_ZDIR_] > max_cfl) max_cfl = local_cfl[_ZDIR_];
 
     _ArrayIncrementIndex_(ndims,dim,index,done);
   }
 
   return(max_cfl);
 }

int NavierStokes3DFlux	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the hyperbolic flux function for the 3D Navier-Stokes equations:

\begin{eqnarray} dir = x, & {\bf f}\left({\bf u}\right) = \left[\begin{array}{c} \rho u \\ \rho u^2 + p \\ \rho u v \\ \rho u w \\ (e+p)u \end{array}\right], \\ dir = y, & {\bf f}\left({\bf u}\right) = \left[\begin{array}{c} \rho v \\ \rho u v \\ \rho v^2 + p \\ \rho v w \\ (e+p)v \end{array}\right], \\ dir = z, & {\bf f}\left({\bf u}\right) = \left[\begin{array}{c} \rho w \\ \rho u w \\ \rho v w \\ \rho w^2 + p \\ (e+p)w \end{array}\right] \end{eqnarray}

Note: the flux function needs to be computed at the ghost points as well.

Parameters

f	Array to hold the computed flux vector (same layout as u)
u	Array with the solution vector
dir	Spatial dimension (x, y, or z) for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 23 of file NavierStokes3DFlux.c.

 {
   HyPar             *solver = (HyPar*)   s;
   NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
   int               i;
 
   int *dim    = solver->dim_local;
   int ghosts  = solver->ghosts;
   int ndims   = solver->ndims;
   int index[ndims], bounds[ndims], offset[ndims];
 
   /* set bounds for array index to include ghost points */
   _ArrayCopy1D_(dim,bounds,ndims);
   for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
 
   /* set offset such that index is compatible with ghost point arrangement */
   _ArraySetValue_(offset,ndims,-ghosts);
 
   int done = 0; _ArraySetValue_(index,ndims,0);
   while (!done) {
     int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
     double rho, vx, vy, vz, e, P;
     _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,vx,vy,vz,e,P,param->gamma);
     _NavierStokes3DSetFlux_((f+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,vx,vy,vz,e,P,dir);
     _ArrayIncrementIndex_(ndims,bounds,index,done);
   }
 
   return(0);
 }

int NavierStokes3DStiffFlux	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the stiff flux, given the solution vector. The stiff flux is the component of the total flux that represents the acoustic modes (see _NavierStokes3DSetStiffFlux_). Here, the linearized approximation to the stiff flux is computed as:

\begin{equation} {\bf f}\left({\bf u}\right) = A_f{\bf u} \end{equation}

where \(A_f = A_f\left({\bf u}_{ref}\right)\) is the fast Jacobian (NavierStokes3D::fast_jac) evaluated for the solution at the beginning of each time step ( \({\bf u}_{ref}\) is NavierStokes3D::solution). This is done in NavierStokes3DPreStep().

Note: the flux function needs to be computed at the ghost points as well.

Parameters

f	Array to hold the computed flux vector (same layout as u)
u	Array with the solution vector
dir	Spatial dimension (x,y, or z) for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 70 of file NavierStokes3DFlux.c.

 {
   HyPar             *solver = (HyPar*)   s;
   NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
   int               *dim    = solver->dim_local;
   int               ghosts  = solver->ghosts;
   static const int  JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
   static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
 
   /* set bounds for array index to include ghost points */
   _ArrayAddCopy1D_(dim,(2*ghosts),bounds,_MODEL_NDIMS_);
   /* set offset such that index is compatible with ghost point arrangement */
   _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
 
   int done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
   while (!done) {
     int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
     double *Af = param->fast_jac+(_MODEL_NDIMS_*p+dir)*JacSize;
     MatVecMult5(_MODEL_NVARS_,(f+_MODEL_NVARS_*p),Af,(u+_MODEL_NVARS_*p));
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
   }
 
   return(0);
 }

int NavierStokes3DNonStiffFlux	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the non-stiff flux, given the solution vector. The non-stiff flux is the component of the total flux that represents the entropy modes (see _NavierStokes3DSetNonStiffFlux_). Here, the linearized approximation to the non-stiff flux is computed as:

\begin{equation} {\bf f}\left({\bf u}\right) - A_f{\bf u} \end{equation}

where \({\bf f}\left({\bf u}\right)\) is the total flux computed in NavierStokes3DFlux(), and \(A_f{\bf u}\) is the linearized stiff flux computed in NavierStokes3DStiffFlux().

Note: the flux function needs to be computed at the ghost points as well.

Parameters

f	Array to hold the computed flux vector (same layout as u)
u	Array with the solution vector
dir	Spatial dimension (x,y, or z) for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 111 of file NavierStokes3DFlux.c.

 {
   HyPar             *solver = (HyPar*)   s;
   NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
   int               *dim    = solver->dim_local;
   int               ghosts  = solver->ghosts;
   static const int  JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
   static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
   static double     ftot[_MODEL_NVARS_], fstiff[_MODEL_NVARS_];
 
   /* set bounds for array index to include ghost points */
   _ArrayAddCopy1D_(dim,(2*ghosts),bounds,_MODEL_NDIMS_);
   /* set offset such that index is compatible with ghost point arrangement */
   _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
 
   int done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
   while (!done) {
     int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
     /* compute total flux */
     double rho, vx, vy, vz, e, P;
     _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,vx,vy,vz,e,P,param->gamma);
     _NavierStokes3DSetFlux_(ftot,_NavierStokes3D_stride_,rho,vx,vy,vz,e,P,dir);
     /* compute stiff stuff */
     double *Af = param->fast_jac+(_MODEL_NDIMS_*p+dir)*JacSize;
     MatVecMult5(_MODEL_NVARS_,fstiff,Af,(u+_MODEL_NVARS_*p));
     /* subtract stiff flux from total flux */
     _ArraySubtract1D_((f+_MODEL_NVARS_*p),ftot,fstiff,_MODEL_NVARS_);
     /* Done */
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
   }
 
   return(0);
 }

int NavierStokes3DRoeAverage	(	double *	uavg,
		double *	uL,
		double *	uR,
		void *	p
	)

Compute the Roe-averaged state for the 3D Navier Stokes equations. This function just calls the macro _NavierStokes3DRoeAverage_ and is not used by any functions within the 3D Navier Stokes module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

uavg	The computed Roe-averaged state
uL	Left state (conserved variables)
uR	Right state (conserved variables)
p	Object of type NavierStokes3D with physics-related variables

Definition at line 18 of file NavierStokes3DFunctions.c.

 {
   NavierStokes3D *param  = (NavierStokes3D*) p;
   _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,uL,uR,param->gamma);
   return(0);
 }

int NavierStokes3DParabolicFunction	(	double *	par,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Compute the viscous terms in the 3D Navier Stokes equations: this function computes the following:

\begin{equation} \frac {\partial} {\partial x} \left[\begin{array}{c} 0 \\ \tau_{xx} \\ \tau_{yx} \\ \tau_{zx} \\ u \tau_{xx} + v \tau_{yx} + w \tau_{zx} - q_x \end{array}\right] + \frac {\partial} {\partial y} \left[\begin{array}{c} 0 \\ \tau_{xy} \\ \tau_{yy} \\ \tau_{zy} \\ u \tau_{xy} + v \tau_{yy} + w \tau_{zy} - q_y \end{array}\right] + \frac {\partial} {\partial z} \left[\begin{array}{c} 0 \\ \tau_{xz} \\ \tau_{yz} \\ \tau_{zz} \\ u \tau_{xz} + v \tau_{yz} + w \tau_{zz} - q_z \end{array}\right] \end{equation}

where

\begin{align} \tau_{xx} &= \frac{2}{3}\left(\frac{\mu}{Re}\right)\left(2\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} - \frac{\partial w}{\partial z}\right),\\ \tau_{xy} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right),\\ \tau_{xz} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}\right),\\ \tau_{yx} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right),\\ \tau_{yy} &= \frac{2}{3}\left(\frac{\mu}{Re}\right)\left(-\frac{\partial u}{\partial x} +2\frac{\partial v}{\partial y} - \frac{\partial w}{\partial z}\right),\\ \tau_{yz} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial w}{\partial y} + \frac{\partial v}{\partial z}\right),\\ \tau_{zx} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}\right),\\ \tau_{zy} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}\right),\\ \tau_{zz} &= \frac{2}{3}\left(\frac{\mu}{Re}\right)\left(-\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} + 2\frac{\partial v}{\partial y}\right),\\ q_x &= -\frac{mu}{\left(\gamma-1\right)Re Pr}\frac{\partial T}{\partial x}, \\ q_y &= -\frac{mu}{\left(\gamma-1\right)Re Pr}\frac{\partial T}{\partial y}, \\ q_z &= -\frac{mu}{\left(\gamma-1\right)Re Pr}\frac{\partial T}{\partial z} \end{align}

and the temperature is \(T = \gamma p/\rho\). \(Re\) and \(Pr\) are the Reynolds and Prandtl numbers, respectively. Note that this function computes the entire parabolic term, and thus bypasses HyPar's parabolic function calculation interfaces. NavierStokes3DInitialize() assigns this function to HyPar::ParabolicFunction.

Reference:

Tannehill, Anderson and Pletcher, Computational Fluid Mechanics and Heat Transfer, Chapter 5, Section 5.1.7 (However, the non-dimensional velocity and the Reynolds number is based on speed of sound, instead of the freestream velocity).

Parameters

par	Array to hold the computed viscous terms
u	Solution vector array
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 50 of file NavierStokes3DParabolicFunction.c.

 {
   HyPar           *solver   = (HyPar*) s;
   MPIVariables    *mpi      = (MPIVariables*) m;
   NavierStokes3D  *physics  = (NavierStokes3D*) solver->physics;
   int             i,j,k,v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int imax   = solver->dim_local[0];
   int jmax   = solver->dim_local[1];
   int kmax   = solver->dim_local[2];
   int *dim   = solver->dim_local;
   int size   = solver->npoints_local_wghosts;
 
   _ArraySetValue_(par,size*_MODEL_NVARS_,0.0);
   if (physics->Re <= 0) return(0); /* inviscid flow */
   solver->count_par++;
 
   static double two_third    = 2.0/3.0;
   double        inv_gamma_m1 = 1.0 / (physics->gamma-1.0);
   double        inv_Re       = 1.0 / physics->Re;
   double        inv_Pr       = 1.0 / physics->Pr;
 
 #if defined(CPU_STAT)
   clock_t startEvent, stopEvent;
 #endif
 
   double *Q; /* primitive variables */
   Q = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
 
 #if defined(CPU_STAT)
   startEvent = clock();
 #endif
 
   for (i=-ghosts; i<(imax+ghosts); i++) {
     for (j=-ghosts; j<(jmax+ghosts); j++) {
       for (k=-ghosts; k<(kmax+ghosts); k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         double energy,pressure;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
         _NavierStokes3DGetFlowVar_( (u+p),
                                     _NavierStokes3D_stride_,
                                     Q[p+0],
                                     Q[p+1],
                                     Q[p+2],
                                     Q[p+3],
                                     energy,
                                     pressure,
                                     physics->gamma);
         Q[p+4] = physics->gamma*pressure/Q[p+0]; /* temperature */
       }
     }
   }
 
 #if defined(CPU_STAT)
   stopEvent = clock();
   printf("%-50s CPU time (secs) = %.6f\n",
          "NaverStokes3DParabolicFunction_Q",(double)(stopEvent-startEvent)/CLOCKS_PER_SEC);
 #endif
 
   double *QDerivX = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
   double *QDerivY = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
   double *QDerivZ = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
 
   IERR solver->FirstDerivativePar(QDerivX,Q,_XDIR_,1,solver,mpi); CHECKERR(ierr);
   IERR solver->FirstDerivativePar(QDerivY,Q,_YDIR_,1,solver,mpi); CHECKERR(ierr);
   IERR solver->FirstDerivativePar(QDerivZ,Q,_ZDIR_,1,solver,mpi); CHECKERR(ierr);
 
   IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
                                  solver->ghosts,mpi,QDerivX); CHECKERR(ierr);
   IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
                                  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
   IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,
                                  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
 
   for (i=-ghosts; i<(imax+ghosts); i++) {
     for (j=-ghosts; j<(jmax+ghosts); j++) {
       for (k=-ghosts; k<(kmax+ghosts); k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         double dxinv, dyinv, dzinv;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
         _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
         _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
         _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
         _ArrayScale1D_((QDerivX+p),dxinv,_MODEL_NVARS_);
         _ArrayScale1D_((QDerivY+p),dyinv,_MODEL_NVARS_);
         _ArrayScale1D_((QDerivZ+p),dzinv,_MODEL_NVARS_);
       }
     }
   }
 
   double *FViscous = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
   double *FDeriv   = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
 
   /* Along X */
   for (i=-ghosts; i<(imax+ghosts); i++) {
     for (j=-ghosts; j<(jmax+ghosts); j++) {
       for (k=-ghosts; k<(kmax+ghosts); k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
 
         double uvel, vvel, wvel, T, Tx,
                ux, uy, uz, vx, vy, wx, wz;
         uvel = (Q+p)[1];
         vvel = (Q+p)[2];
         wvel = (Q+p)[3];
         T    = (Q+p)[4];
         Tx   = (QDerivX+p)[4];
         ux   = (QDerivX+p)[1];
         vx   = (QDerivX+p)[2];
         wx   = (QDerivX+p)[3];
         uy   = (QDerivY+p)[1];
         vy   = (QDerivY+p)[2];
         uz   = (QDerivZ+p)[1];
         wz   = (QDerivZ+p)[3];
 
         /* calculate viscosity coeff based on Sutherland's law */
         double mu = raiseto(T, 0.76);
         /*
         if (p/_MODEL_NVARS_ == 49841) {
           printf("[CPU] 49841 mu = %.20e T = %.20e log(0.76) = %.20e T*log(0.76) = %.20e\n",
                  mu, T, log(0.76), T*log(0.76));
         }
         */
 
         double tau_xx, tau_xy, tau_xz, qx;
         tau_xx = two_third * (mu*inv_Re) * (2*ux - vy - wz);
         tau_xy = (mu*inv_Re) * (uy + vx);
         tau_xz = (mu*inv_Re) * (uz + wx);
         qx     = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tx;
 
         (FViscous+p)[0] = 0.0;
         (FViscous+p)[1] = tau_xx;
         (FViscous+p)[2] = tau_xy;
         (FViscous+p)[3] = tau_xz;
         (FViscous+p)[4] = uvel*tau_xx + vvel*tau_xy + wvel*tau_xz + qx;
       }
     }
   }
 
   IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
   for (i=0; i<imax; i++) {
     for (j=0; j<jmax; j++) {
       for (k=0; k<kmax; k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         double dxinv;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
         _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
         for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dxinv * (FDeriv+p)[v] );
       }
     }
   }
 
   /* Along Y */
   for (i=-ghosts; i<(imax+ghosts); i++) {
     for (j=-ghosts; j<(jmax+ghosts); j++) {
       for (k=-ghosts; k<(kmax+ghosts); k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
 
         double uvel, vvel, wvel, T, Ty,
                ux, uy, vx, vy, vz, wy, wz;
         uvel = (Q+p)[1];
         vvel = (Q+p)[2];
         wvel = (Q+p)[3];
         T    = (Q+p)[4];
         Ty   = (QDerivY+p)[4];
         ux   = (QDerivX+p)[1];
         vx   = (QDerivX+p)[2];
         uy   = (QDerivY+p)[1];
         vy   = (QDerivY+p)[2];
         wy   = (QDerivY+p)[3];
         vz   = (QDerivZ+p)[2];
         wz   = (QDerivZ+p)[3];
 
         /* calculate viscosity coeff based on Sutherland's law */
         double mu = raiseto(T, 0.76);
 
         double tau_yx, tau_yy, tau_yz, qy;
         tau_yx = (mu*inv_Re) * (uy + vx);
         tau_yy = two_third * (mu*inv_Re) * (-ux + 2*vy - wz);
         tau_yz = (mu*inv_Re) * (vz + wy);
         qy     = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Ty;
 
         (FViscous+p)[0] = 0.0;
         (FViscous+p)[1] = tau_yx;
         (FViscous+p)[2] = tau_yy;
         (FViscous+p)[3] = tau_yz;
         (FViscous+p)[4] = uvel*tau_yx + vvel*tau_yy + wvel*tau_yz + qy;
       }
     }
   }
 
   IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
   for (i=0; i<imax; i++) {
     for (j=0; j<jmax; j++) {
       for (k=0; k<kmax; k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         double dyinv;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
         _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
         for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
       }
     }
   }
 
   /* Along Z */
   for (i=-ghosts; i<(imax+ghosts); i++) {
     for (j=-ghosts; j<(jmax+ghosts); j++) {
       for (k=-ghosts; k<(kmax+ghosts); k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
 
         double uvel, vvel, wvel, T, Tz,
                ux, uz, vy, vz, wx, wy, wz;
         uvel = (Q+p)[1];
         vvel = (Q+p)[2];
         wvel = (Q+p)[3];
         T    = (Q+p)[4];
         Tz   = (QDerivZ+p)[4];
         ux   = (QDerivX+p)[1];
         wx   = (QDerivX+p)[3];
         vy   = (QDerivY+p)[2];
         wy   = (QDerivY+p)[3];
         uz   = (QDerivZ+p)[1];
         vz   = (QDerivZ+p)[2];
         wz   = (QDerivZ+p)[3];
 
         /* calculate viscosity coeff based on Sutherland's law */
         double mu = raiseto(T,0.76);
 
         double tau_zx, tau_zy, tau_zz, qz;
         tau_zx = (mu*inv_Re) * (uz + wx);
         tau_zy = (mu*inv_Re) * (vz + wy);
         tau_zz = two_third * (mu*inv_Re) * (-ux - vy + 2*wz);
         qz     = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tz;
 
         (FViscous+p)[0] = 0.0;
         (FViscous+p)[1] = tau_zx;
         (FViscous+p)[2] = tau_zy;
         (FViscous+p)[3] = tau_zz;
         (FViscous+p)[4] = uvel*tau_zx + vvel*tau_zy + wvel*tau_zz + qz;
       }
     }
   }
 
   IERR solver->FirstDerivativePar(FDeriv,FViscous,_ZDIR_,-1,solver,mpi); CHECKERR(ierr);
   for (i=0; i<imax; i++) {
     for (j=0; j<jmax; j++) {
       for (k=0; k<kmax; k++) {
         int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
         double dzinv;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
         _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
         for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dzinv * (FDeriv+p)[v] );
       }
     }
   }
 
   free(Q);
   free(QDerivX);
   free(QDerivY);
   free(QDerivZ);
   free(FViscous);
   free(FDeriv);
 
   if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,_MODEL_NVARS_);
   return(0);
 }

int NavierStokes3DSource	(	double *	source,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Computes the gravitational source term using a well-balanced formulation: the source term is rewritten as follows:

\begin{equation} \left[\begin{array}{c} 0 \\ -\rho {\bf g}\cdot{\bf \hat{i}} \\ -\rho {\bf g}\cdot{\bf \hat{j}} \\ -\rho {\bf g}\cdot{\bf \hat{k}} \\ -\rho u {\bf g}\cdot{\bf \hat{i}} - \rho v {\bf g}\cdot{\bf \hat{j}} - \rho w {\bf g}\cdot{\bf \hat{k}} \end{array}\right] = \left[\begin{array}{ccccc} 0 & p_0 \varrho^{-1} & 0 & 0 & p_0 u \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial x}\left[\begin{array}{c} 0 \\ \varphi \\ 0 \\ 0 \\ \varphi \end{array}\right] + \left[\begin{array}{ccccc} 0 & 0 & p_0 \varrho^{-1} & 0 & p_0 v \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial y}\left[\begin{array}{c} 0 \\ 0 \\ \varphi \\ 0 \\ \varphi \end{array}\right] + \left[\begin{array}{ccccc} 0 & 0 & 0 & p_0 \varrho^{-1} & p_0 w \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial y}\left[\begin{array}{c} 0 \\ 0 \\ 0 \\ \varphi \\ \varphi \end{array}\right] \end{equation}

where \(\varphi = \varphi\left(x,y\right)\) and \(\varrho = \varrho\left(x,y\right)\) are computed in NavierStokes3DGravityField(). The derivatives are computed in an exactly identical manner as the hyperbolic flux (i.e., using a conservative finite-difference formulation) (see HyperbolicFunction()).

References:

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

source	Array to hold the computed source
u	Solution vector array
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 38 of file NavierStokes3DSource.c.

 {
   HyPar           *solver = (HyPar* )         s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
 
   if ((param->grav_x == 0.0) && (param->grav_y == 0.0) && (param->grav_z == 0.0))
     return(0); /* no gravitational forces */
 
   int     v, done, p, p1, p2, dir;
   double  *SourceI = solver->fluxI; /* interace source term       */
   double  *SourceC = solver->fluxC; /* cell-centered source term  */
   double  *SourceL = solver->fL;
   double  *SourceR = solver->fR;
 
   int     ghosts  = solver->ghosts;
   int     *dim    = solver->dim_local;
   double  *x      = solver->x;
   double  *dxinv  = solver->dxinv;
   double  RT      =  param->p0 / param->rho0;
   static int index[_MODEL_NDIMS_],index1[_MODEL_NDIMS_],index2[_MODEL_NDIMS_],dim_interface[_MODEL_NDIMS_];
   static double grav[_MODEL_NDIMS_];
 
   grav[_XDIR_] = param->grav_x;
   grav[_YDIR_] = param->grav_y;
   grav[_ZDIR_] = param->grav_z;
   for (dir = 0; dir < _MODEL_NDIMS_; dir++) {
     if (grav[dir] != 0.0) {
       /* set interface dimensions */
       _ArrayCopy1D_(dim,dim_interface,_MODEL_NDIMS_); dim_interface[dir]++;
       /* calculate the split source function exp(-phi/RT) */
       IERR NavierStokes3DSourceFunction(SourceC,u,x,solver,mpi,t,dir); CHECKERR(ierr);
       /* calculate the left and right interface source terms */
       IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,dir,solver,mpi,0); CHECKERR(ierr);
       IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,dir,solver,mpi,0); CHECKERR(ierr);
       /* calculate the final interface source term */
       IERR NavierStokes3DSourceUpwind(SourceI,SourceL,SourceR,u,dir,solver,t);
       /* calculate the final cell-centered source term */
       done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
       while (!done) {
         _ArrayCopy1D_(index,index1,_MODEL_NDIMS_);
         _ArrayCopy1D_(index,index2,_MODEL_NDIMS_); index2[dir]++;
         _ArrayIndex1D_(_MODEL_NDIMS_,dim          ,index ,ghosts,p );
         _ArrayIndex1D_(_MODEL_NDIMS_,dim_interface,index1,0     ,p1);
         _ArrayIndex1D_(_MODEL_NDIMS_,dim_interface,index2,0     ,p2);
 
         double dx_inverse; _GetCoordinate_(dir,index[dir],dim,ghosts,dxinv,dx_inverse);
         double rho, vel[_MODEL_NDIMS_], e, P; 
         _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],e,P,param->gamma);
         double term[_MODEL_NVARS_] = {0.0, rho*RT*(dir==_XDIR_), rho*RT*(dir==_YDIR_), rho*RT*(dir==_ZDIR_), rho*RT*vel[dir]};
         for (v=0; v<_MODEL_NVARS_; v++) {
           source[_MODEL_NVARS_*p+v] += (  (term[v]*param->grav_field_f[p])
                                         * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse );
         }
         vel[0] = P; /* useless statement to avoid compiler warnings */
         _ArrayIncrementIndex_(_MODEL_NDIMS_,dim,index,done);
       }
     }
   }
 
   return(0);
 }

int NavierStokes3DJacobian	(	double *	Jac,
		double *	u,
		void *	p,
		int	dir,
		int	nvars,
		int	upw
	)

Function to compute the flux Jacobian of the 3D Navier-Stokes equations, given the solution at a grid point. The Jacobian is square matrix of size nvar=5, and is returned as a 1D array (double) of 25 elements in row-major format.

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 25
u	solution at a grid point (array of size nvar = 5)
p	object containing the physics-related parameters
dir	dimension (0 -> x, 1 -> y, 2 -> z)
nvars	number of vector components
upw	0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 15 of file NavierStokes3DJacobian.c.

 {
   NavierStokes3D *param = (NavierStokes3D*) p;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];
 
   /* get the eigenvalues and left,right eigenvectors */
   _NavierStokes3DEigenvalues_      (u,_NavierStokes3D_stride_,D,param->gamma,dir);
   _NavierStokes3DLeftEigenvectors_ (u,_NavierStokes3D_stride_,L,param->gamma,dir);
   _NavierStokes3DRightEigenvectors_(u,_NavierStokes3D_stride_,R,param->gamma,dir);
 
   int aupw = absolute(upw), k;
   k = 0;  D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   k = 6;  D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   k = 12; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   k = 18; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
 
   MatMult5(_MODEL_NVARS_,DL,D,L);
   MatMult5(_MODEL_NVARS_,Jac,R,DL);
 
   return(0);
 }

int NavierStokes3DStiffJacobian	(	double *	Jac,
		double *	u,
		void *	p,
		int	dir,
		int	nvars,
		int	upw
	)

Function to compute the Jacobian of the fast flux (representing the acoustic waves) of the 3D Navier-Stokes equations, given the solution at a grid point. The Jacobian is square matrix of size nvar=5, and is returned as a 1D array (double) of 25 elements in row-major format.

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 25
u	solution at a grid point (array of size nvar = 5)
p	object containing the physics-related parameters
dir	dimension (0 -> x, 1 -> y, 2 -> z)
nvars	number of vector components
upw	0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 53 of file NavierStokes3DJacobian.c.

 {
   NavierStokes3D *param = (NavierStokes3D*) p;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];
 
   /* get the eigenvalues and left,right eigenvectors */
   _NavierStokes3DEigenvalues_      (u,_NavierStokes3D_stride_,D,param->gamma,dir);
   _NavierStokes3DLeftEigenvectors_ (u,_NavierStokes3D_stride_,L,param->gamma,dir);
   _NavierStokes3DRightEigenvectors_(u,_NavierStokes3D_stride_,R,param->gamma,dir);
 
   int aupw = absolute(upw), k;
   if (dir == _XDIR_) {
     k = 0;  D[k] = 0.0;
     k = 6;  D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
     k = 12; D[k] = 0.0;
     k = 18; D[k] = 0.0;
     k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   } else if (dir == _YDIR_) {
     k = 0;  D[k] = 0.0;
     k = 6;  D[k] = 0.0;
     k = 12; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
     k = 18; D[k] = 0.0;
     k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   } else if (dir == _ZDIR_) {
     k = 0;  D[k] = 0.0;
     k = 6;  D[k] = 0.0;
     k = 12; D[k] = 0.0;
     k = 18; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
     k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
   }
 
   MatMult5(_MODEL_NVARS_,DL,D,L);
   MatMult5(_MODEL_NVARS_,Jac,R,DL);
 
   return(0);
 }

int NavierStokes3DLeftEigenvectors	(	double *	u,
		double *	L,
		void *	p,
		int	dir
	)

Compute the left eigenvections for the 3D Navier Stokes equations. This function just calls the macro _NavierStokes3DLeftEigenvectors_ and is not used by any functions within the 3D Navier Stokes module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

u	Conserved solution at a grid point
L	Array of size nvar^2 = 5^2 to save the matrix of left eigenvectors in (row-major format).
p	Object of type NavierStokes3D with physics-related variables
dir	Spatial dimension (x, y, or z)

Definition at line 20 of file NavierStokes3DEigen.c.

 {
   NavierStokes3D *param  = (NavierStokes3D*)  p;
   _NavierStokes3DLeftEigenvectors_(u,_NavierStokes3D_stride_,L,param->gamma,dir);
   return(0);
 }

int NavierStokes3DRightEigenvectors	(	double *	u,
		double *	R,
		void *	p,
		int	dir
	)

Compute the right eigenvections for the 3D Navier Stokes equations. This function just calls the macro _NavierStokes3DRightEigenvectors_ and is not used by any functions within the 3D Navier Stokes module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

u	Conserved solution at a grid point
R	Array of size nvar^2 = 5^2 to save the matrix of right eigenvectors in (row-major format).
p	Object of type NavierStokes3D with physics-related variables
dir	Spatial dimension (x, y, or z)

Definition at line 39 of file NavierStokes3DEigen.c.

 {
   NavierStokes3D *param  = (NavierStokes3D*)  p;
   _NavierStokes3DRightEigenvectors_(u,_NavierStokes3D_stride_,R,param->gamma,dir);
   return(0);
 }

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

Roe's upwinding scheme.

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

Roe, P. L., “Approximate Riemann solvers, parameter vectors, and difference schemes,” Journal of Computational Physics, Vol. 43, No. 2, 1981, pp. 357–372, http://dx.doi.org/10.1016/0021-9991(81)90128-5.

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

Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x, y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 40 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)    s;
   NavierStokes3D  *param  = (NavierStokes3D*)  solver->physics;
   int             done;
 
   int *dim  = solver->dim_local;
 
   int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
   _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
                 modA[_MODEL_NVARS_*_MODEL_NVARS_];
 
   done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
 
 #if defined(CPU_STAT)
   clock_t startEvent, stopEvent;
   startEvent = clock();
 #endif
 
   while (!done) {
     _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
       int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
 
       /* Roe's upwinding scheme */
 
       udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
       udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
       udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
       udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
       udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
 
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DEigenvalues_       (uavg,dummy,D,param->gamma,dir);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
 
       /* Harten's Entropy Fix - Page 362 of Leveque */
       int k;
       double delta = 0.000001, delta2 = delta*delta;
       k=0;  D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=6;  D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
 
       MatMult5(_MODEL_NVARS_,DL,D,L);
       MatMult5(_MODEL_NVARS_,modA,R,DL);
       MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
 
       fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
       fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
       fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
       fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
       fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
 #if defined(CPU_STAT)
   stopEvent = clock();
   printf("%-50s CPU time (secs) = %.6f dir = %d\n",
           "NavierStokes3DUpwindRoe", (double)(stopEvent-startEvent)/CLOCKS_PER_SEC, dir);
 #endif
 
   return(0);
 }

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

Characteristic-based Roe-fixed upwinding scheme.

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

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

C.-W. Shu, and S. Osher, "Efficient implementation of essentially non-oscillatory schemes, II", J. Comput. Phys., 83 (1989), pp. 32–78, http://dx.doi.org/10.1016/0021-9991(89)90222-2.

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

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x, y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 140 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)    s;
   NavierStokes3D  *param  = (NavierStokes3D*)  solver->physics;
   int             done,k;
 
   int *dim  = solver->dim_local;
 
   int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
   _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
 
   done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
   while (!done) {
     _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
       int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
              ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
 
       /* Roe-Fixed upwinding scheme */
 
       _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
 
       _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
 
       /* calculate characteristic fluxes and variables */
       MatVecMult5(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
       MatVecMult5(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
       MatVecMult5(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
       MatVecMult5(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
 
       double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
       _NavierStokes3DEigenvalues_((uL+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
       eigL[0] = D[0];
       eigL[1] = D[6];
       eigL[2] = D[12];
       eigL[3] = D[18];
       eigL[4] = D[24];
       _NavierStokes3DEigenvalues_((uR+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
       eigR[0] = D[0];
       eigR[1] = D[6];
       eigR[2] = D[12];
       eigR[3] = D[18];
       eigR[4] = D[24];
       _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
       eigC[0] = D[0];
       eigC[1] = D[6];
       eigC[2] = D[12];
       eigC[3] = D[18];
       eigC[4] = D[24];
 
       for (k = 0; k < _MODEL_NVARS_; k++) {
         if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0))      fc[k] = fcL[k];
         else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
         else {
           double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
           fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
         }
       }
 
       /* calculate the interface flux from the characteristic flux */
       MatVecMult5(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

Characteristic-based local Lax-Friedrich upwinding scheme.

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

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

C.-W. Shu, and S. Osher, "Efficient implementation of essentially non-oscillatory schemes, II", J. Comput. Phys., 83 (1989), pp. 32–78, http://dx.doi.org/10.1016/0021-9991(89)90222-2.

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

Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x, y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 244 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)    s;
   NavierStokes3D  *param  = (NavierStokes3D*)  solver->physics;
   int             done;
 
   int *dim  = solver->dim_local;
 
   int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
   _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
 
   done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
   while (!done) {
     _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
       int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
              ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
 
       /* Local Lax-Friedrich upwinding scheme */
 
       _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
 
       _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
 
       /* calculate characteristic fluxes and variables */
       MatVecMult5(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
       MatVecMult5(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
       MatVecMult5(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
       MatVecMult5(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
 
       double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
       _NavierStokes3DEigenvalues_((uL+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
       eigL[0] = D[0];
       eigL[1] = D[6];
       eigL[2] = D[12];
       eigL[3] = D[18];
       eigL[4] = D[24];
       _NavierStokes3DEigenvalues_((uR+_MODEL_NVARS_*p),_NavierStokes3D_stride_,D,param->gamma,dir);
       eigR[0] = D[0];
       eigR[1] = D[6];
       eigR[2] = D[12];
       eigR[3] = D[18];
       eigR[4] = D[24];
       _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
       eigC[0] = D[0];
       eigC[1] = D[6];
       eigC[2] = D[12];
       eigC[3] = D[18];
       eigC[4] = D[24];
 
       double alpha;
       alpha = max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
       fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
       alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
       fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
       alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
       fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
       alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
       fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
       alpha = max3(absolute(eigL[4]),absolute(eigC[4]),absolute(eigR[4]));
       fc[4] = 0.5 * (fcL[4] + fcR[4] + alpha * (ucL[4]-ucR[4]));
 
       /* calculate the interface flux from the characteristic flux */
       MatVecMult5(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

Rusanov's upwinding scheme.

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

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

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

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

Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x,y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 349 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)          s;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
   int             *dim    = solver->dim_local, done;
 
   int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
   bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
 
   static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
 
   static int index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
              indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_];
 
   done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
   while (!done) {
     _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
 
       int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       int q = p*_MODEL_NVARS_;
 
       /* Modified Rusanov's upwinding scheme */
 
       udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
       udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
       udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
       udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
       udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
 
       _NavierStokes3DRoeAverage_(uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
 
       double c, vel[_MODEL_NDIMS_], rho,E,P;
       _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaL = c + absolute(vel[dir]);
       _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaR = c + absolute(vel[dir]);
       _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaavg = c + absolute(vel[dir]);
 
       double kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
       double alpha  = kappa*max3(alphaL,alphaR,alphaavg);
 
       fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-alpha*udiff[0];
       fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-alpha*udiff[1];
       fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-alpha*udiff[2];
       fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-alpha*udiff[3];
       fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-alpha*udiff[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

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

Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x, y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 430 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)    s;
   NavierStokes3D  *param  = (NavierStokes3D*)  solver->physics;
   int             done;
 
   int     *dim  = solver->dim_local;
   double  *uref = param->solution;
 
   int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
   _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
                 modA[_MODEL_NVARS_*_MODEL_NVARS_];
 
   done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
   while (!done) {
     _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
       int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
 
       /* Roe's upwinding scheme */
 
       udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
       udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
       udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
       udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
       udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
 
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DEigenvalues_       (uavg,dummy,D,param->gamma,dir);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
 
       /* Harten's Entropy Fix - Page 362 of Leveque */
       int k;
       double delta = 0.000001, delta2 = delta*delta;
       if (dir == _XDIR_) {
         k=0;  D[k] = 0.0;
         k=6;  D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
         k=12; D[k] = 0.0;
         k=18; D[k] = 0.0;
         k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       } else if (dir == _YDIR_) {
         k=0;  D[k] = 0.0;
         k=6;  D[k] = 0.0;
         k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
         k=18; D[k] = 0.0;
         k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       } else if (dir == _ZDIR_) {
         k=0;  D[k] = 0.0;
         k=6;  D[k] = 0.0;
         k=12; D[k] = 0.0;
         k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
         k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       }
 
       MatMult5(_MODEL_NVARS_,DL,D,L);
       MatMult5(_MODEL_NVARS_,modA,R,DL);
       MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
 
       fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
       fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
       fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
       fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
       fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

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

Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x, y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 527 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)    s;
   NavierStokes3D  *param  = (NavierStokes3D*)  solver->physics;
   int             done;
 
   int     *dim  = solver->dim_local;
   double  *uref = param->solution;
 
   int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] =  1;
   _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
                 modA[_MODEL_NVARS_*_MODEL_NVARS_];
 
   done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
   while (!done) {
     _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
       int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       int k;
       double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_],
              udiss_total[_MODEL_NVARS_],udiss_stiff[_MODEL_NVARS_];
       double delta = 0.000001, delta2 = delta*delta;
 
       /* Roe's upwinding scheme */
 
       udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
       udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
       udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
       udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
       udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
 
       /* Compute total dissipation */
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DEigenvalues_       (uavg,dummy,D,param->gamma,dir);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
       k=0;  D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=6;  D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       MatMult5(_MODEL_NVARS_,DL,D,L);
       MatMult5(_MODEL_NVARS_,modA,R,DL);
       MatVecMult5(_MODEL_NVARS_,udiss_total,modA,udiff);
 
       /* Compute dissipation corresponding to acoustic modes */
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DEigenvalues_       (uavg,dummy,D,param->gamma,dir);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
       if (dir == _XDIR_) {
         k=0;  D[k] = 0.0;
         k=6;  D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
         k=12; D[k] = 0.0;
         k=18; D[k] = 0.0;
         k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       } else if (dir == _YDIR_) {
         k=0;  D[k] = 0.0;
         k=6;  D[k] = 0.0;
         k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
         k=18; D[k] = 0.0;
         k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       } else if (dir == _ZDIR_) {
         k=0;  D[k] = 0.0;
         k=6;  D[k] = 0.0;
         k=12; D[k] = 0.0;
         k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
         k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
       }
       MatMult5(_MODEL_NVARS_,DL,D,L);
       MatMult5(_MODEL_NVARS_,modA,R,DL);
       MatVecMult5(_MODEL_NVARS_,udiss_stiff,modA,udiff);
 
      /* Compute the dissipation term for the entropy modes */
       _ArraySubtract1D_(udiss,udiss_total,udiss_stiff,_MODEL_NVARS_);
 
       fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
       fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
       fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
       fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
       fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

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

Ghosh, D., Constantinescu, E. M., "Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning", Submitted (http://arxiv.org/abs/1510.05751).

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x,y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 638 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)          s;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
   int             *dim    = solver->dim_local, done;
 
   static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
   bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
 
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
                 modA[_MODEL_NVARS_*_MODEL_NVARS_];
 
   static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
              index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
 
   static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
 
   done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
   while (!done) {
     _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
 
       int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       int q = p*_MODEL_NVARS_;
 
       /* Modified Rusanov's upwinding scheme */
 
       udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
       udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
       udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
       udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
       udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
 
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
 
       double c, vel[_MODEL_NDIMS_], rho,E,P;
 
       _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaL = c + absolute(vel[dir]);
       double betaL = absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaR = c + absolute(vel[dir]);
       double betaR = absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaavg = c + absolute(vel[dir]);
       double betaavg = absolute(vel[dir]);
 
       double kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
       double alpha  = kappa*max3(alphaL,alphaR,alphaavg);
       double beta   = kappa*max3(betaL,betaR,betaavg);
 
       _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
       if (dir == _XDIR_) {
         D[0]  = beta;
         D[6]  = alpha;
         D[12] = beta;
         D[18] = beta;
         D[24] = alpha;
       } else if (dir == _YDIR_) {
         D[0]  = beta;
         D[6]  = beta;
         D[12] = alpha;
         D[18] = beta;
         D[24] = alpha;
       } else if (dir == _ZDIR_) {
         D[0]  = beta;
         D[6]  = beta;
         D[12] = beta;
         D[18] = alpha;
         D[24] = alpha;
       }
       MatMult5    (_MODEL_NVARS_,DL,D,L);
       MatMult5    (_MODEL_NVARS_,modA,R,DL);
       MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
 
       fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
       fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
       fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
       fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
       fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

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

Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x,y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 756 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)          s;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
   int             *dim    = solver->dim_local, done;
   double          *uref   = param->solution;
 
   static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
   bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
 
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
                 modA[_MODEL_NVARS_*_MODEL_NVARS_];
 
   static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
              index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
 
   static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
 
   done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
   while (!done) {
     _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
 
       int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       int q = p*_MODEL_NVARS_;
 
       /* Modified Rusanov's upwinding scheme */
 
       udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
       udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
       udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
       udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
       udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
 
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
 
       double c, vel[_MODEL_NDIMS_], rho,E,P;
 
       _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaL = c + absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaR = c + absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       double alphaavg = c + absolute(vel[dir]);
 
       double kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
       double alpha  = kappa*max3(alphaL,alphaR,alphaavg);
 
       _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
       if (dir == _XDIR_) {
         D[6]  = alpha;
         D[24] = alpha;
       } else if (dir == _YDIR_) {
         D[12] = alpha;
         D[24] = alpha;
       } else if (dir == _ZDIR_) {
         D[18] = alpha;
         D[24] = alpha;
       }
       MatMult5    (_MODEL_NVARS_,DL,D,L);
       MatMult5    (_MODEL_NVARS_,modA,R,DL);
       MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
 
       fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
       fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
       fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
       fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
       fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

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

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

Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x,y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 862 of file NavierStokes3DUpwind.c.

 {
   HyPar           *solver = (HyPar*)          s;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
   int             *dim    = solver->dim_local, done;
   double          *uref   = param->solution;
 
   static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
   bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
   bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
 
   static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                 L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_],
                 modA[_MODEL_NVARS_*_MODEL_NVARS_];
 
   static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
              index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
 
   static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
 
   done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
   while (!done) {
     _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
     for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
 
       int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
       _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
       _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
       int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
       int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
       int q = p*_MODEL_NVARS_;
       double udiff[_MODEL_NVARS_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_],
              udiss_total[_MODEL_NVARS_],udiss_acoustic[_MODEL_NVARS_];
 
       /* Modified Rusanov's upwinding scheme */
       /* Modified Rusanov's upwinding scheme */
       double c, vel[_MODEL_NDIMS_], rho,E,P, alphaL, alphaR, alphaavg,
              betaL, betaR, betaavg, alpha, beta,
              kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
 
       udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
       udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
       udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
       udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
       udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
 
       /* Compute total dissipation */
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
 
       _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       alphaL = c + absolute(vel[dir]);
       betaL = absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       alphaR = c + absolute(vel[dir]);
       betaR = absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       alphaavg = c + absolute(vel[dir]);
       betaavg = absolute(vel[dir]);
 
       kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
       alpha  = kappa*max3(alphaL,alphaR,alphaavg);
       beta   = kappa*max3(betaL,betaR,betaavg);
 
       _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
       if (dir == _XDIR_) {
         D[0]  = beta;
         D[6]  = alpha;
         D[12] = beta;
         D[18] = beta;
         D[24] = alpha;
       } else if (dir == _YDIR_) {
         D[0]  = beta;
         D[6]  = beta;
         D[12] = alpha;
         D[18] = beta;
         D[24] = alpha;
       } else if (dir == _ZDIR_) {
         D[0]  = beta;
         D[6]  = beta;
         D[12] = beta;
         D[18] = alpha;
         D[24] = alpha;
       }
       MatMult5    (_MODEL_NVARS_,DL,D,L);
       MatMult5    (_MODEL_NVARS_,modA,R,DL);
       MatVecMult5 (_MODEL_NVARS_,udiss_total,modA,udiff);
 
       /* Compute dissipation for the linearized acoustic modes */
       _NavierStokes3DRoeAverage_        (uavg,_NavierStokes3D_stride_,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
       _NavierStokes3DLeftEigenvectors_  (uavg,dummy,L,param->gamma,dir);
       _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
 
       _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       alphaL = c + absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       alphaR = c + absolute(vel[dir]);
 
       _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
       c = sqrt(param->gamma*P/rho);
       alphaavg = c + absolute(vel[dir]);
 
       kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
       alpha  = kappa*max3(alphaL,alphaR,alphaavg);
 
       _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
       if (dir == _XDIR_) {
         D[6]  = alpha;
         D[24] = alpha;
       } else if (dir == _YDIR_) {
         D[12] = alpha;
         D[24] = alpha;
       } else if (dir == _ZDIR_) {
         D[18] = alpha;
         D[24] = alpha;
       }
       MatMult5    (_MODEL_NVARS_,DL,D,L);
       MatMult5    (_MODEL_NVARS_,modA,R,DL);
       MatVecMult5 (_MODEL_NVARS_,udiss_acoustic,modA,udiff);
 
       /* Compute dissipation for the entropy modes */
       _ArraySubtract1D_(udiss,udiss_total,udiss_acoustic,_MODEL_NVARS_);
 
       fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
       fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
       fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
       fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
       fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
     }
     _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
   }
 
   return(0);
 }

int NavierStokes3DGravityField	(	void *	s,
		void *	m
	)

This function computes the pressure and density variation functions for the hydrostatic balance of the type specified by the user (NavierStokes3D:HB). The pressure and density in hydrostatic balance are given by

\begin{equation} \rho = \rho_0\varrho\left(x,y\right),\ p = p_0\varphi\left(x,y\right) \end{equation}

where \(\rho_0\) and \(p_0\) are the reference density (NavierStokes3D::rho0) and pressure (NavierStokes3D::p0). This function computes \(\varrho\) (NavierStokes3D::grav_field_f) and \(\varphi\) (NavierStokes3D::grav_field_g). For flows without gravity, \(\varrho = \varphi = 1\).

References:

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

s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 33 of file NavierStokes3DGravityField.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
 
   double  *f      = param->grav_field_f;
   double  *g      = param->grav_field_g;
   int     *dim    = solver->dim_local;
   int     ghosts  = solver->ghosts;
   int     index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_],
           offset[_MODEL_NDIMS_], d, done;
 
   /* set bounds for array index to include ghost points */
   _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_);
   for (d=0; d<_MODEL_NDIMS_; d++) bounds[d] += 2*ghosts;
   /* set offset such that index is compatible with ghost point arrangement */
   _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
 
   double p0   = param->p0;
   double rho0 = param->rho0;
   double RT   = p0 / rho0;
   double gamma= param->gamma;
   double R    = param->R;
   double Cp   = gamma*R / (gamma-1.0);
   double T0   = p0 / (rho0 * R);
   double gx   = param->grav_x;
   double gy   = param->grav_y;
   double gz   = param->grav_z;
   double Nbv  = param->N_bv;
   
   /* set the value of the gravity field */
   done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
   if (param->HB == 1) {
     while (!done) {
       int p;         _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
       double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
       double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
       double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
       f[p] = exp( (gx*xcoord+gy*ycoord+gz*zcoord)/RT);
       g[p] = exp(-(gx*xcoord+gy*ycoord+gz*zcoord)/RT);
       _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
     }
   } else if (param->HB == 2) {
     while (!done) {
       int p;         _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
       double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
       double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
       double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
       f[p] = raiseto((1.0-(gx*xcoord+gy*ycoord+gz*zcoord)/(Cp*T0)), (-1.0 /(gamma-1.0)));
       g[p] = raiseto((1.0-(gx*xcoord+gy*ycoord+gz*zcoord)/(Cp*T0)), (gamma/(gamma-1.0)));
       _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
     }
   } else if (param->HB == 3) {
     while (!done) {
       int p;         _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
       double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
       if ((gx != 0) || (gy != 0)) {
         if (!mpi->rank) {
           fprintf(stderr,"Error in NavierStokes3DGravityField(): HB = 3 is implemented only for ");
           fprintf(stderr,"gravity force along the z-coordinate.\n");
         }
         return(1);
       }
       double Pexner = 1 + ((gz*gz)/(Cp*T0*Nbv*Nbv)) * (exp(-(Nbv*Nbv/gz)*zcoord)-1.0);
       f[p] = raiseto(Pexner, (-1.0 /(gamma-1.0))) * exp(Nbv*Nbv*zcoord/gz);
       g[p] = raiseto(Pexner, (gamma/(gamma-1.0)));
       _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
     }
   } else {
     while (!done) {
       int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
       f[p] = g[p] = 1.0;
       _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
     }
   }
 
   /* A sensible simulation will not specify peridic boundary conditions 
    * along a direction in which gravity acts.
    * Gravity will be zero along the dimension periodic BCs are specified,
    * so the value of the gravity potential will be the same along those
    * grid lines.
    * Thus, for both these cases, extrapolate the gravity field at the 
    * boundaries */
   int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_];
   for (d = 0; d < _MODEL_NDIMS_; d++) {
     /* left boundary */
     if (!mpi->ip[d]) {
       _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
       _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
       done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
       while (!done) {
         _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = ghosts-1-indexb[d];
         int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
         int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
         f[p1] = f[p2];
         g[p1] = g[p2];
         _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
       }
     }
     /* right boundary */
     if (mpi->ip[d] == mpi->iproc[d]-1) {
       _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
       _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
       done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
       while (!done) {
         _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = dim[d]-1-indexb[d];
         int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
         int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
         f[p1] = f[p2];
         g[p1] = g[p2];
         _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
       }
     }
   }
 
   return(0);
 }

int NavierStokes3DModifiedSolution	(	double *	uC,
		double *	u,
		int	d,
		void *	s,
		void *	m,
		double	waqt
	)

This function computes the modified solution for the well-balanced treatment of the gravitational source terms. The modified solution vector is given by

\begin{equation} {\bf u}^* = \left[\begin{array}{c} \rho \varrho^{-1}\left(x,y\right) \\ \rho u \varrho^{-1}\left(x,y\right) \\ \rho v \varrho^{-1}\left(x,y\right) \\ \rho w \varrho^{-1}\left(x,y\right) \\ e^* \end{array}\right] \end{equation}

where

\begin{equation} e^* = \frac{p \varphi^{-1}\left(x,y\right)}{\gamma-1} + \frac{1}{2}\rho \varrho^{-1}\left(x,y\right) \left(u^2+v^2+w^2\right) \end{equation}

\(\varrho\) and \(\varphi\) are computed in NavierStokes3DGravityField(). For flows without gravity, \({\bf u}^* = {\bf u}\).

References:

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

uC	Array to hold the computed modified solution
u	Solution vector array
d	spatial dimension (not used)
s	Solver object of type HyPar
m	MPI object of time MPIVariables
waqt	Current simulation time

Definition at line 31 of file NavierStokes3DModifiedSolution.c.

 {
   HyPar           *solver = (HyPar*)          s;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
 
   int     ghosts  = solver->ghosts;
   int     *dim    = solver->dim_local;
   int     ndims   = solver->ndims;
   int     index[ndims], bounds[ndims], offset[ndims];
 
   /* set bounds for array index to include ghost points */
   _ArrayCopy1D_(dim,bounds,ndims);
   int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
 
   /* set offset such that index is compatible with ghost point arrangement */
   _ArraySetValue_(offset,ndims,-ghosts);
 
   int done = 0; _ArraySetValue_(index,ndims,0);
   double inv_gamma_m1 = 1.0 / (param->gamma-1.0);
   while (!done) {
     int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
     double rho, uvel, vvel, wvel, E, P;
     _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,E,P,param->gamma);
     uC[_MODEL_NVARS_*p+0] = u[_MODEL_NVARS_*p+0] * param->grav_field_f[p];
     uC[_MODEL_NVARS_*p+1] = u[_MODEL_NVARS_*p+1] * param->grav_field_f[p];
     uC[_MODEL_NVARS_*p+2] = u[_MODEL_NVARS_*p+2] * param->grav_field_f[p];
     uC[_MODEL_NVARS_*p+3] = u[_MODEL_NVARS_*p+3] * param->grav_field_f[p];
     uC[_MODEL_NVARS_*p+4] = (P*inv_gamma_m1)*(1.0/param->grav_field_g[p]) + (0.5*rho*(uvel*uvel+vvel*vvel+wvel*wvel))*param->grav_field_f[p];
     _ArrayIncrementIndex_(ndims,bounds,index,done);
   }
 
   return(0);
 }

int NavierStokes3DIBAdiabatic	(	void *	s,
		void *	m,
		double *	u,
		double	t
	)

Apply no-slip adiabatic wall boundary conditions on the immersed boundary points (grid points within the immersed body that are within stencil-width distance of interior points, i.e., points in the interior of the computational domain).

Parameters

s	Solver object of type HyPar
m	Solver object of type HyPar
u	Array with the solution vector
t	Current simulation time

Definition at line 19 of file NavierStokes3DImmersedBoundary.c.

 {
   HyPar             *solver   = (HyPar*)   s;
   MPIVariables      *mpi      = (MPIVariables*) m;
   ImmersedBoundary  *IB       = (ImmersedBoundary*) solver->ib;
   IBNode            *boundary = IB->boundary;
   NavierStokes3D    *param    = (NavierStokes3D*) solver->physics;
   static double     v[_MODEL_NVARS_];
   int               n, j, k, nb = IB->n_boundary_nodes;
 
   if (!solver->flag_ib) return(0);
 
   /* Ideally, this shouldn't be here - But this function is called everywhere
      (through ApplyIBConditions()) *before* MPIExchangeBoundariesnD is called! */
   MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,solver->ghosts,mpi,u);
 
   double inv_gamma_m1 = 1.0 / (param->gamma - 1.0);
 
   double ramp_fac = 1.0;
   if (param->t_ib_ramp > 0) {
     double x = t/param->t_ib_ramp;
     if (!strcmp(param->ib_ramp_type,_IB_RAMP_LINEAR_)) {
       ramp_fac = x;
     } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_SMOOTHEDSLAB_)) {
       double a = 0.0;
       double b = 1.0;
       double c = 0.5;
       double r = param->t_ib_ramp/param->t_ib_width;
       ramp_fac = (a*exp(c*r)+b*exp(r*x))/(exp(c*r)+exp(r*x));
     } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_DISABLE_)) {
       ramp_fac = 0.0;
     } else {
       fprintf(stderr,"Error in NavierStokes3DImmersedBoundary():\n");
       fprintf(stderr,"  Ramp type %s not recognized.\n", param->ib_ramp_type);
       return 1;
     }
   }
 
   for (n=0; n<nb; n++) {
 
     int     node_index = boundary[n].p;
     double  *alpha = &(boundary[n].interp_coeffs[0]);
     int     *nodes = &(boundary[n].interp_nodes[0]);
     double  factor = boundary[n].surface_distance / boundary[n].interp_node_distance;
 
     _ArraySetValue_(v,_MODEL_NVARS_,0.0);
     for (j=0; j<_IB_NNODES_; j++) {
       for (k=0; k<_MODEL_NVARS_; k++) {
         v[k] += ( alpha[j] * u[_MODEL_NVARS_*nodes[j]+k] );
       }
     }
 
     double rho, uvel, vvel, wvel, energy, pressure;
     _NavierStokes3DGetFlowVar_(v,_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,param->gamma);
 
     double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
     _NavierStokes3DGetFlowVar_( (u+_MODEL_NVARS_*node_index),
                                 _NavierStokes3D_stride_,
                                 rho_gpt,
                                 uvel_gpt,
                                 vvel_gpt,
                                 wvel_gpt,
                                 energy_gpt,
                                 pressure_gpt,
                                 param->gamma );
 
     double rho_ib_target, uvel_ib_target, vvel_ib_target, wvel_ib_target, pressure_ib_target;
     rho_ib_target = rho;
     pressure_ib_target = pressure;
     uvel_ib_target = -uvel * factor;
     vvel_ib_target = -vvel * factor;
     wvel_ib_target = -wvel * factor;
 
     double rho_ib, uvel_ib, vvel_ib, wvel_ib, energy_ib, pressure_ib;
     rho_ib      = ramp_fac * rho_ib_target      + (1.0-ramp_fac) * rho_gpt;
     pressure_ib = ramp_fac * pressure_ib_target + (1.0-ramp_fac) * pressure_gpt;
     uvel_ib     = ramp_fac * uvel_ib_target     + (1.0-ramp_fac) * uvel_gpt;
     vvel_ib     = ramp_fac * vvel_ib_target     + (1.0-ramp_fac) * vvel_gpt;
     wvel_ib     = ramp_fac * wvel_ib_target     + (1.0-ramp_fac) * wvel_gpt;
     energy_ib   = inv_gamma_m1*pressure_ib 
                   + 0.5*rho_ib*(uvel_ib*uvel_ib+vvel_ib*vvel_ib+wvel_ib*wvel_ib);
 
     u[_MODEL_NVARS_*node_index+0] = rho_ib;
     u[_MODEL_NVARS_*node_index+1] = rho_ib * uvel_ib;
     u[_MODEL_NVARS_*node_index+2] = rho_ib * vvel_ib;
     u[_MODEL_NVARS_*node_index+3] = rho_ib * wvel_ib;
     u[_MODEL_NVARS_*node_index+4] = energy_ib;
   }
 
   return 0;
 }

int NavierStokes3DIBIsothermal	(	void *	s,
		void *	m,
		double *	u,
		double	t
	)

Apply no-slip isothermal wall boundary conditions on the immersed boundary points (grid points within the immersed body that are within stencil-width distance of interior points, i.e., points in the interior of the computational domain).

Parameters

s	Solver object of type HyPar
m	Solver object of type HyPar
u	Array with the solution vector
t	Current simulation time

Definition at line 119 of file NavierStokes3DImmersedBoundary.c.

 {
   HyPar             *solver   = (HyPar*)   s;
   MPIVariables      *mpi      = (MPIVariables*) m;
   ImmersedBoundary  *IB       = (ImmersedBoundary*) solver->ib;
   IBNode            *boundary = IB->boundary;
   NavierStokes3D    *param    = (NavierStokes3D*) solver->physics;
   static double     v[_MODEL_NVARS_];
   int               n, j, k, nb = IB->n_boundary_nodes;
 
   if (!solver->flag_ib) return(0);
 
   /* Ideally, this shouldn't be here - But this function is called everywhere
      (through ApplyIBConditions()) *before* MPIExchangeBoundariesnD is called! */
   MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,solver->ghosts,mpi,u);
 
   double inv_gamma_m1 = 1.0 / (param->gamma - 1.0);
 
   double ramp_fac = 1.0;
   if (param->t_ib_ramp > 0) {
     double x = t/param->t_ib_ramp;
     if (!strcmp(param->ib_ramp_type,_IB_RAMP_LINEAR_)) {
       ramp_fac = x;
     } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_SMOOTHEDSLAB_)) {
       double a = 0.0;
       double b = 1.0;
       double c = 0.5;
       double r = param->t_ib_ramp/param->t_ib_width;
       ramp_fac = (a*exp(c*r)+b*exp(r*x))/(exp(c*r)+exp(r*x));
     } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_DISABLE_)) {
       ramp_fac = 0.0;
     } else {
       fprintf(stderr,"Error in NavierStokes3DImmersedBoundary():\n");
       fprintf(stderr,"  Ramp type %s not recognized.\n", param->ib_ramp_type);
       return 1;
     }
   }
 
   for (n=0; n<nb; n++) {
 
     int     node_index = boundary[n].p;
     double  *alpha = &(boundary[n].interp_coeffs[0]);
     int     *nodes = &(boundary[n].interp_nodes[0]);
     double  factor = boundary[n].surface_distance / boundary[n].interp_node_distance;
 
     _ArraySetValue_(v,_MODEL_NVARS_,0.0);
     for (j=0; j<_IB_NNODES_; j++) {
       for (k=0; k<_MODEL_NVARS_; k++) {
         v[k] += ( alpha[j] * u[_MODEL_NVARS_*nodes[j]+k] );
       }
     }
 
     double rho, uvel, vvel, wvel, energy, pressure, temperature;
     _NavierStokes3DGetFlowVar_(v,_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,param->gamma);
     temperature = pressure / rho;
 
     double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt, temperature_gpt;
     _NavierStokes3DGetFlowVar_( (u+_MODEL_NVARS_*node_index),
                                 _NavierStokes3D_stride_,
                                 rho_gpt,
                                 uvel_gpt,
                                 vvel_gpt,
                                 wvel_gpt,
                                 energy_gpt,
                                 pressure_gpt,
                                 param->gamma );
     temperature_gpt = pressure_gpt / rho_gpt;
 
     double  rho_ib_target, 
             uvel_ib_target, vvel_ib_target, wvel_ib_target, 
             pressure_ib_target, 
             temperature_ib_target;
     temperature_ib_target = (1.0+factor)*param->T_ib_wall - factor * temperature;
     if (    (temperature_ib_target < param->T_ib_wall/param->ib_T_tol)
          || (temperature_ib_target > param->T_ib_wall*param->ib_T_tol) ) {
       temperature_ib_target = param->T_ib_wall;
     }
     pressure_ib_target = pressure;
     rho_ib_target = pressure_ib_target / temperature_ib_target;
     uvel_ib_target = - factor * uvel;
     vvel_ib_target = - factor * vvel;
     wvel_ib_target = - factor * wvel;
 
     double rho_ib, uvel_ib, vvel_ib, wvel_ib, energy_ib, pressure_ib;
     rho_ib      = ramp_fac * rho_ib_target      + (1.0-ramp_fac) * rho_gpt;
     pressure_ib = ramp_fac * pressure_ib_target + (1.0-ramp_fac) * pressure_gpt;
     uvel_ib     = ramp_fac * uvel_ib_target     + (1.0-ramp_fac) * uvel_gpt;
     vvel_ib     = ramp_fac * vvel_ib_target     + (1.0-ramp_fac) * vvel_gpt;
     wvel_ib     = ramp_fac * wvel_ib_target     + (1.0-ramp_fac) * wvel_gpt;
     energy_ib   = inv_gamma_m1*pressure_ib 
                   + 0.5*rho_ib*(uvel_ib*uvel_ib+vvel_ib*vvel_ib+wvel_ib*wvel_ib);
 
     u[_MODEL_NVARS_*node_index+0] = rho_ib;
     u[_MODEL_NVARS_*node_index+1] = rho_ib * uvel_ib;
     u[_MODEL_NVARS_*node_index+2] = rho_ib * vvel_ib;
     u[_MODEL_NVARS_*node_index+3] = rho_ib * wvel_ib;
     u[_MODEL_NVARS_*node_index+4] = energy_ib;
   }
 
   return(0);
 }

int NavierStokes3DPreStep	(	double *	u,
		void *	s,
		void *	m,
		double	waqt
	)

Pre-step function for the 3D Navier Stokes equations: This function is called at the beginning of each time step.

The solution at the beginning of the step is copied into NavierStokes3D::solution for the linearized flux partitioning.
The linearized "fast" Jacobian (representing the acoustic modes) NavierStokes3D::fast_jac is computed as:
\begin{equation} A_f\left({\bf u}\right) = R\left({\bf u}\right)\Lambda_f\left({\bf u}\right)L\left({\bf u}\right) \end{equation}
where \(R\) and \(L\) are the matrices of right and left eigenvectors, and,
\begin{equation} \Lambda_f = diag\left[0,0,0,u+c,u-c \right] \end{equation}
with \(c\) being the speed of sound.

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

u	Solution vector
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current simulation time

Definition at line 32 of file NavierStokes3DPreStep.c.

 {
   HyPar             *solver = (HyPar*)   s;
   NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
   int               *dim    = solver->dim_local;
   int               ghosts  = solver->ghosts, dir, p;
   double            *A;
 
   static const int  ndims   = _MODEL_NDIMS_;
   static const int  JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
   static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
   static double     D[_MODEL_NVARS_*_MODEL_NVARS_],L[_MODEL_NVARS_*_MODEL_NVARS_],
                     R[_MODEL_NVARS_*_MODEL_NVARS_],DL[_MODEL_NVARS_*_MODEL_NVARS_];
 
   /* set bounds for array index to include ghost points */
   _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
   /* set offset such that index is compatible with ghost point arrangement */
   _ArraySetValue_(offset,ndims,-ghosts);
   /* copy the solution to act as a reference for linearization */
   _ArrayCopy1D_(u,param->solution,(solver->npoints_local_wghosts*_MODEL_NVARS_));
 
   int done = 0; _ArraySetValue_(index,ndims,0);
   while (!done) {
     _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
     int q = _MODEL_NVARS_*p;
 
     dir = _XDIR_;
     A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
     /* get the eigenvalues, and left & right eigenvectors */
     _NavierStokes3DEigenvalues_      ((u+q),_NavierStokes3D_stride_,D,param->gamma,dir);
     _NavierStokes3DLeftEigenvectors_ ((u+q),_NavierStokes3D_stride_,L,param->gamma,dir);
     _NavierStokes3DRightEigenvectors_((u+q),_NavierStokes3D_stride_,R,param->gamma,dir);
     /* remove the entropy modes (corresponding to eigenvalues u) */
     D[0] = D[12] = D[18] = 0.0;
     /* assemble the Jacobian */
     MatMult5(_MODEL_NVARS_,DL,D,L );
     MatMult5(_MODEL_NVARS_,A ,R,DL);
 
     dir = _YDIR_;
     A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
     /* get the eigenvalues, and left & right eigenvectors */
     _NavierStokes3DEigenvalues_      ((u+q),_NavierStokes3D_stride_,D,param->gamma,dir);
     _NavierStokes3DLeftEigenvectors_ ((u+q),_NavierStokes3D_stride_,L,param->gamma,dir);
     _NavierStokes3DRightEigenvectors_((u+q),_NavierStokes3D_stride_,R,param->gamma,dir);
     /* remove the entropy modes (corresponding to eigenvalues v) */
     D[0] = D[6] = D[18] = 0.0;
     /* assemble the Jacobian */
     MatMult5(_MODEL_NVARS_,DL,D,L );
     MatMult5(_MODEL_NVARS_,A ,R,DL);
 
     dir = _ZDIR_;
     A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
     /* get the eigenvalues, and left & right eigenvectors */
     _NavierStokes3DEigenvalues_      ((u+q),_NavierStokes3D_stride_,D,param->gamma,dir);
     _NavierStokes3DLeftEigenvectors_ ((u+q),_NavierStokes3D_stride_,L,param->gamma,dir);
     _NavierStokes3DRightEigenvectors_((u+q),_NavierStokes3D_stride_,R,param->gamma,dir);
     /* remove the entropy modes (corresponding to eigenvalues v) */
     D[0] = D[6] = D[12] = 0.0;
     /* assemble the Jacobian */
     MatMult5(_MODEL_NVARS_,DL,D,L );
     MatMult5(_MODEL_NVARS_,A ,R,DL);
 
     _ArrayIncrementIndex_(ndims,bounds,index,done);
   }
   return(0);
 }

int NavierStokes3DIBForces	(	void *	s,
		void *	m,
		double	a_t
	)

Calculate the aerodynamic forces on the immersed body surface and write them to file

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
a_t	Current simulation time

Definition at line 264 of file NavierStokes3DIBForces.c.

 {
   HyPar             *solver  = (HyPar*)          s;
   MPIVariables      *mpi     = (MPIVariables*)   m; 
   NavierStokes3D    *physics = (NavierStokes3D*) solver->physics;
   ImmersedBoundary  *IB      = (ImmersedBoundary*) solver->ib;
   int ierr;
 
   if (!solver->flag_ib) return(0);
 
   int npts = solver->npoints_local_wghosts;
 
   double* pressure = (double*) calloc (npts, sizeof(double));
   ierr = NavierStokes3DComputePressure(pressure, solver->u, solver);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  NavierStokes3DComputePressure() returned with error.\n");
     return 1;
   }
   double* temperature = (double*) calloc(npts, sizeof(double));
   ierr = NavierStokes3DComputeTemperature(temperature, solver->u, solver);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  NavierStokes3DComputeTemperature() returned with error.\n");
     return 1;
   }
 
   /* Compute surface pressure */
   double* p_surface = NULL;
   ierr = IBComputeFacetVar(solver, mpi, pressure, 1, &p_surface);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  IBComputeFacetVar() returned with error.\n");
     return 1;
   }
   /* Compute surface temperature */
   double* T_surface = NULL;
   ierr = IBComputeFacetVar(solver, mpi, temperature, 1, &T_surface);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  IBComputeFacetVar() returned with error.\n");
     return 1;
   }
   /* Compute normal surface pressure gradient */
   double *ngrad_p_surface = NULL;
   ierr = IBComputeNormalGradient(solver, mpi, pressure, 1, &ngrad_p_surface);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  IBComputeNormalGradient() returned with error.\n");
     return 1;
   }
   /* Compute normal temperature gradient */
   double *ngrad_T_surface = NULL;
   ierr = IBComputeNormalGradient(solver, mpi, temperature, 1, &ngrad_T_surface);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  IBComputeNormalGradient() returned with error.\n");
     return 1;
   }
   /* Compute shear forces */
   double *shear = NULL;
   ierr = ComputeShear(solver, mpi, solver->u, &shear);
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  ComputeShear() returned with error.\n");
     return 1;
   }
 
   char surface_filename[_MAX_STRING_SIZE_] = "surface";
   if (solver->nsims == 1) {
     if (!strcmp(solver->op_overwrite,"no")) {
       strcat(surface_filename,solver->filename_index);
     }
   } else {
     char index[_MAX_STRING_SIZE_];
     GetStringFromInteger(solver->my_idx, index, (int)log10(solver->nsims)+1);
     strcat(surface_filename, "_");
     strcat(surface_filename, index);
     strcat(surface_filename, "_");
   }
   strcat(surface_filename,".dat");
   if (!mpi->rank) {
     printf("Writing immersed body surface data file %s.\n",surface_filename);
   }
   ierr = WriteSurfaceData(  mpi,
                             IB,
                             p_surface,
                             T_surface,
                             ngrad_p_surface,
                             ngrad_T_surface,
                             shear,
                             surface_filename );
   if (ierr) {
     fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
     fprintf(stderr,"  WriteSurfaceData() returned with error\n");
     return 1;
   }
 
   free(pressure);
   free(temperature);
   if (p_surface) free(p_surface);
   if (T_surface) free(T_surface);
   if (ngrad_p_surface) free(ngrad_p_surface);
   if (ngrad_T_surface) free(ngrad_T_surface);
   if (shear) free(shear);
 
   return 0;
 }

int gpuNavierStokes3DInitialize	(	void *	s,
		void *	m
	)

Initialize GPU-related arrays.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 10 of file NavierStokes3DInitialize_GPU.c.

 {
     HyPar          *solver   = (HyPar*)          s;
     NavierStokes3D *physics  = (NavierStokes3D*) solver->physics;
 
     int *dim = solver->dim_local;
     int ghosts = solver->ghosts;
     int d, size = 1; for (d = 0; d <_MODEL_NDIMS_; d++) size *= (dim[d] + 2*ghosts);
 
     gpuMalloc((void**)&physics->gpu_Q, size*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_QDerivX, size*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_QDerivY, size*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_QDerivZ, size*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_FViscous, size*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_FDeriv, size*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_grav_field_f, size*sizeof(double));
     gpuMalloc((void**)&physics->gpu_grav_field_g, size*sizeof(double));
     gpuMalloc((void**)&physics->gpu_fast_jac, _MODEL_NDIMS_*size*_MODEL_NVARS_*_MODEL_NVARS_*sizeof(double));
     gpuMalloc((void**)&physics->gpu_solution, size*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_Q, 0, size*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_QDerivX, 0, size*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_QDerivY, 0, size*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_QDerivZ, 0, size*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_FViscous, 0, size*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_FDeriv, 0, size*_MODEL_NVARS_*sizeof(double));
     gpuMemcpy(physics->gpu_grav_field_f, physics->grav_field_f, size*sizeof(double), gpuMemcpyHostToDevice);
     gpuMemcpy(physics->gpu_grav_field_g, physics->grav_field_g, size*sizeof(double), gpuMemcpyHostToDevice);
     gpuMemset(physics->gpu_fast_jac, 0, _MODEL_NDIMS_*size*_MODEL_NVARS_*_MODEL_NVARS_*sizeof(double));
     gpuMemset(physics->gpu_solution, 0, size*_MODEL_NVARS_*sizeof(double));
 
     return (0);
 }

int gpuNavierStokes3DFlux	(	double *	,
		double *	,
		int	,
		void *	,
		double
	)

int gpuNavierStokes3DParabolicFunction	(	double *	,
		double *	,
		void *	,
		void *	,
		double
	)

int gpuNavierStokes3DSource	(	double *	,
		double *	,
		void *	,
		void *	,
		double
	)

int gpuNavierStokes3DModifiedSolution	(	double *	,
		double *	,
		int	,
		void *	,
		void *	,
		double
	)

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

Rusanov's upwinding scheme.

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

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

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

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

Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension (x,y, or z)
s	Solver object of type HyPar
t	Current solution time

Definition at line 340 of file NavierStokes3DUpwind_GPU.cu.

 {
   HyPar           *solver = (HyPar*)          s;
   NavierStokes3D  *param  = (NavierStokes3D*) solver->physics;
   int             *dim    = solver->dim_local;
 
   int bounds_inter[_MODEL_NDIMS_];
   _ArrayCopy1D_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
   int npoints_grid; _ArrayProduct1D_(bounds_inter,_MODEL_NDIMS_,npoints_grid);
   int nblocks = (npoints_grid-1)/GPU_THREADS_PER_BLOCK + 1;
 
   gpuNavierStokes3DUpwindRusanov_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
     npoints_grid, solver->npoints_local_wghosts, dir, solver->ghosts, param->gamma, solver->gpu_dim_local,
     param->gpu_grav_field_g, fL, fR, uL, uR, u, fI
   );
   cudaDeviceSynchronize();
 
   return 0;
 }

int gpuNavierStokes3DUpwindRoe	(	double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		int	,
		void *	,
		double
	)

int gpuNavierStokes3DPreStep	(	double *	u,
		void *	s,
		void *	m,
		double	waqt
	)

Pre-step function for the 3D Navier Stokes equations: This function is called at the beginning of each time step.

The solution at the beginning of the step is copied into NavierStokes3D::solution for the linearized flux partitioning.
The linearized "fast" Jacobian (representing the acoustic modes) NavierStokes3D::fast_jac is computed as:
\begin{equation} A_f\left({\bf u}\right) = R\left({\bf u}\right)\Lambda_f\left({\bf u}\right)L\left({\bf u}\right) \end{equation}
where \(R\) and \(L\) are the matrices of right and left eigenvectors, and,
\begin{equation} \Lambda_f = diag\left[0,0,0,u+c,u-c \right] \end{equation}
with \(c\) being the speed of sound.

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

u	Solution vector
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current simulation time

Definition at line 192 of file NavierStokes3DPreStep_GPU.cu.

 {
   HyPar             *solver = (HyPar*)   s;
   NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
 
   int nblocks = (solver->npoints_local_wghosts-1)/GPU_THREADS_PER_BLOCK + 1;
 
   gpuArrayCopy1D(u,param->gpu_solution,(solver->npoints_local_wghosts)*_MODEL_NVARS_);
   gpuNavierStokes3DPreStep_kernel<<<nblocks, GPU_THREADS_PER_BLOCK>>>(
     solver->npoints_local_wghosts, param->gamma, u, param->gpu_fast_jac
   );
   cudaDeviceSynchronize();
 
   return 0;
 }

int NavierStokes3DInitialize	(	void *	s,
		void *	m
	)

Initialize the 3D Navier-Stokes (NavierStokes3D) module: Sets the default parameters, read in and set physics-related parameters, and set the physics-related function pointers in HyPar.

This file reads the file "physics.inp" that must have the following format:

begin
    <keyword>   <value>
    <keyword>   <value>
    <keyword>   <value>
    ...
    <keyword>   <value>
end

where the list of keywords are:

Keyword name	Type	Variable	Default value
gamma	double	NavierStokes3D::gamma	1.4
Pr	double	NavierStokes3D::Pr	0.72
Re	double	NavierStokes3D::Re	-1
Minf	double	NavierStokes3D::Minf	1.0
gravity	double,double,double	NavierStokes3D::grav_x,NavierStokes3D::grav_y,NavierStokes3D::grav_z	0.0,0.0,0.0
rho_ref	double	NavierStokes3D::rho0	1.0
p_ref	double	NavierStokes3D::p0	1.0
HB	int	NavierStokes3D::HB	1
R	double	NavierStokes3D::R	1.0
upwinding	char[]	NavierStokes3D::upw_choice	"roe" (_ROE_)
ib_wall_type	char[]	NavierStokes3D::ib_wall_type	"adiabatic" (_IB_ADIABATIC_)
ib_ramp_time	double	NavierStokes3D::t_ib_ramp	-1
ib_ramp_width	double	NavierStokes3D::t_ib_width	0.05
ib_ramp_type	char[]	NavierStokes3D::ib_ramp_type	"linear" (_IB_RAMP_LINEAR_)

if "ib_wall_type" (NavierStokes3D::ib_wall_type) is specified as "isothermal", it should be followed by the wall temperature (#NavierStokes3D::T_ib_wall), i.e,

begin ... ib_wall_type isothermal 1.0 ... end
If "HB" (NavierStokes3D::HB) is specified as 3, it should be followed by the the Brunt-Vaisala frequency (NavierStokes3D::N_bv), i.e.

begin ... HB 3 0.01 ... end

Note: "physics.inp" is optional; if absent, default values will be used.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 116 of file NavierStokes3DInitialize.c.

 {
   HyPar           *solver  = (HyPar*)          s;
   MPIVariables    *mpi     = (MPIVariables*)   m;
   NavierStokes3D  *physics = (NavierStokes3D*) solver->physics;
   int             ferr     = 0;
 
   static int count = 0;
 
   if (solver->nvars != _MODEL_NVARS_) {
     fprintf(stderr,"Error in NavierStokes3DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
     return(1);
   }
   if (solver->ndims != _MODEL_NDIMS_) {
     fprintf(stderr,"Error in NavierStokes3DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
     return(1);
   }
 
   /* default values */
   physics->gamma     = 1.4;
   physics->Pr        = 0.72;
   physics->Re        = -1;
   physics->Minf      = 1.0;
   physics->C1        = 1.458e-6;
   physics->C2        = 110.4;
   physics->grav_x    = 0.0;
   physics->grav_y    = 0.0;
   physics->grav_z    = 0.0;
   physics->rho0      = 1.0;
   physics->p0        = 1.0;
   physics->HB        = 1;
   physics->R         = 1.0;
   physics->N_bv      = 0.0;
   physics->T_ib_wall = -DBL_MAX;
   physics->t_ib_ramp = -1.0;
   physics->t_ib_width= 0.05;
   physics->ib_T_tol  = 5;
   strcpy(physics->upw_choice,"roe");
   strcpy(physics->ib_write_surface_data,"yes");
   strcpy(physics->ib_wall_type,"adiabatic");
   strcpy(physics->ib_ramp_type,"linear");
 
   /* reading physical model specific inputs - all processes */
   if (!mpi->rank) {
     FILE *in;
     if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
     in = fopen("physics.inp","r");
     if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
     else {
       char word[_MAX_STRING_SIZE_];
       ferr = fscanf(in,"%s",word);
       if (ferr != 1) {
         fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
         return 1;
       }
       if (!strcmp(word, "begin")){
           while (strcmp(word, "end")){
               ferr = fscanf(in,"%s",word);
           if (ferr != 1) {
             fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
             return 1;
           }
           if (!strcmp(word, "gamma")) {
             ferr = fscanf(in,"%lf",&physics->gamma);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"upwinding")) {
             ferr = fscanf(in,"%s",physics->upw_choice);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"Pr")) {
             ferr = fscanf(in,"%lf",&physics->Pr);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"Re")) {
             ferr = fscanf(in,"%lf",&physics->Re);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"Minf")) {
             ferr = fscanf(in,"%lf",&physics->Minf);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"gravity")) {
             ferr = fscanf(in,"%lf",&physics->grav_x);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             ferr = fscanf(in,"%lf",&physics->grav_y);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             ferr = fscanf(in,"%lf",&physics->grav_z);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"rho_ref")) {
             ferr = fscanf(in,"%lf",&physics->rho0);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"p_ref")) {
             ferr = fscanf(in,"%lf",&physics->p0);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"HB")) {
             ferr = fscanf(in,"%d",&physics->HB);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             if (physics->HB==3) {
               ferr = fscanf(in,"%lf",&physics->N_bv);
               if (ferr != 1) {
                 fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
                 return 1;
               }
             }
           } else if (!strcmp(word,"R")) {
             ferr = fscanf(in,"%lf",&physics->R);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"ib_surface_data")) {
             ferr = fscanf(in,"%s",physics->ib_write_surface_data);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
           } else if (!strcmp(word,"ib_wall_type")) {
             ferr = fscanf(in,"%s",physics->ib_wall_type);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             if (!strcmp(physics->ib_wall_type,_IB_ISOTHERMAL_)) {
               ferr = fscanf(in,"%lf",&physics->T_ib_wall);
               if (ferr != 1) {
                 fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
                 return 1;
               }
             }
             if (!solver->flag_ib) {
               printf("Warning: in NavierStokes3DInitialize().\n");
               printf("Warning: no immersed body present; specification of ib_wall_type unnecessary.\n");
             }
           } else if (!strcmp(word,"ib_ramp_time")) {
             ferr = fscanf(in,"%lf",&physics->t_ib_ramp);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             if (!solver->flag_ib) {
               printf("Warning: in NavierStokes3DInitialize().\n");
               printf("Warning: no immersed body present; specification of ib_ramp_time unnecessary.\n");
             }
           } else if (!strcmp(word,"ib_ramp_width")) {
             ferr = fscanf(in,"%lf",&physics->t_ib_width);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             if (!solver->flag_ib) {
               printf("Warning: in NavierStokes3DInitialize().\n");
               printf("Warning: no immersed body present; specification of ib_ramp_width unnecessary.\n");
             }
           } else if (!strcmp(word,"ib_ramp_type")) {
             ferr = fscanf(in,"%s",physics->ib_ramp_type);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             if (!solver->flag_ib) {
               printf("Warning: in NavierStokes3DInitialize().\n");
               printf("Warning: no immersed body present; specification of ib_ramp_type unnecessary.\n");
             }
           } else if (!strcmp(word,"ib_T_tolerance")) {
             ferr = fscanf(in,"%lf",&physics->ib_T_tol);
             if (ferr != 1) {
               fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
               return 1;
             }
             if (!solver->flag_ib) {
               printf("Warning: in NavierStokes3DInitialize().\n");
               printf("Warning: no immersed body present; specification of ib_T_tolerance unnecessary.\n");
             }
           } else if (strcmp(word,"end")) {
             char useless[_MAX_STRING_SIZE_];
             ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
             printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
             printf("recognized or extraneous. Ignoring.\n");
           }
         }
         } else {
           fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
         return(1);
         }
     }
     fclose(in);
   }
 
   IERR MPIBroadcast_character (physics->upw_choice            ,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_character (physics->ib_write_surface_data ,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_character (physics->ib_wall_type          ,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_character (physics->ib_ramp_type          ,_MAX_STRING_SIZE_,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->gamma                ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->Pr                   ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->Re                   ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->Minf                 ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->grav_x               ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->grav_y               ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->grav_z               ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->rho0                 ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->p0                   ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->R                    ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->N_bv                 ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->T_ib_wall            ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->t_ib_ramp            ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->t_ib_width           ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_double    (&physics->ib_T_tol             ,1                ,0,&mpi->world); CHECKERR(ierr);
   IERR MPIBroadcast_integer   (&physics->HB                   ,1                ,0,&mpi->world); CHECKERR(ierr);
 
   /* if file output is disabled in HyPar, respect that */
   if (!strcmp(solver->op_file_format,"none")) {
     if (!strcmp(physics->ib_write_surface_data,"yes")) {
       if (!mpi->rank) {
         printf("Warning from NavierStokes3DInitialize(): solver->op_file_format is set to \"none\", thus ");
         printf("setting physics->ib_write_surface_data to \"no\" (no solution files will be written).\n");
       }
     }
     strcpy(physics->ib_write_surface_data,"no");
   }
 
   /* Scaling Re by M_inf */
   physics->Re /= physics->Minf;
 
   /* check that a well-balanced upwinding scheme is being used for cases with gravity */
   if (   ((physics->grav_x != 0.0) || (physics->grav_y != 0.0) || (physics->grav_z != 0.0) )
       && (strcmp(physics->upw_choice,_RUSANOV_)) ) {
     if (!mpi->rank) {
       fprintf(stderr,"Error in NavierStokes3DInitialize: %s upwinding is needed for flows ",_RUSANOV_);
       fprintf(stderr,"with gravitational forces.\n");
     }
     return(1);
   }
   /* check that solver has the correct choice of diffusion formulation, if viscous flow */
   if (strcmp(solver->spatial_type_par,_NC_2STAGE_) && (physics->Re > 0)) {
     if (!mpi->rank) {
       fprintf(stderr,"Error in NavierStokes3DInitialize(): Parabolic term spatial discretization must be \"%s\"\n",_NC_2STAGE_);
     }
     return(1);
   }
 
   /* initializing physical model-specific functions */
 #if defined(HAVE_CUDA)
   if (solver->use_gpu) {
     solver->PreStep     = gpuNavierStokes3DPreStep;
     solver->FFunction   = gpuNavierStokes3DFlux;
     solver->SFunction   = gpuNavierStokes3DSource;
     solver->UFunction   = gpuNavierStokes3DModifiedSolution;
   } else {
 #endif
     solver->PreStep               = NavierStokes3DPreStep;
     solver->ComputeCFL            = NavierStokes3DComputeCFL;
     solver->FFunction             = NavierStokes3DFlux;
     solver->SFunction             = NavierStokes3DSource;
     solver->UFunction             = NavierStokes3DModifiedSolution;
     solver->AveragingFunction     = NavierStokes3DRoeAverage;
     solver->GetLeftEigenvectors   = NavierStokes3DLeftEigenvectors;
     solver->GetRightEigenvectors  = NavierStokes3DRightEigenvectors;
 #if defined(HAVE_CUDA)
   }
 #endif
 
   if (solver->flag_ib) {
     if (!strcmp(physics->ib_wall_type,_IB_ADIABATIC_)) {
       solver->IBFunction = NavierStokes3DIBAdiabatic;
     } else if (!strcmp(physics->ib_wall_type,_IB_ISOTHERMAL_)) {
       solver->IBFunction = NavierStokes3DIBIsothermal;
     } else {
       fprintf(stderr, "Error in NavierStokes3DInitialize()\n");
       fprintf(stderr, "  invalid value for IB wall type (%s).\n",
               physics->ib_wall_type );
     }
     if (!strcmp(physics->ib_write_surface_data,"yes")) {
       solver->PhysicsOutput = NavierStokes3DIBForces;
     }
   }
 
 #if defined(HAVE_CUDA)
   if (solver->use_gpu) {
 
     if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
       if (!mpi->rank) {
         fprintf(stderr,"Error in NavierStokes3DInitialize(): Not available on GPU!");
       }
       return 1;
     } else {
       solver->JFunction = NavierStokes3DJacobian;
       if (!strcmp(physics->upw_choice,_ROE_    )) {
         solver->Upwind = gpuNavierStokes3DUpwindRoe;
       } else if (!strcmp(physics->upw_choice,_RUSANOV_)) {
         solver->Upwind = gpuNavierStokes3DUpwindRusanov;
       } else {
         if (!mpi->rank) {
           fprintf(stderr,"Error in NavierStokes3DInitialize(): %s is not a valid upwinding scheme on GPU. ",
                   physics->upw_choice);
           fprintf(stderr,"Choices are %s and %s.\n",_ROE_,_RUSANOV_);
         }
         return(1);
       }
     }
 
   } else {
 #endif
 
     if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
       solver->FdFFunction = NavierStokes3DNonStiffFlux;
       solver->dFFunction  = NavierStokes3DStiffFlux;
       solver->JFunction   = NavierStokes3DStiffJacobian;
       if (!strcmp(physics->upw_choice,_ROE_)) {
         solver->Upwind    = NavierStokes3DUpwindRoe;
         solver->UpwinddF  = NavierStokes3DUpwinddFRoe;
         solver->UpwindFdF = NavierStokes3DUpwindFdFRoe;
       } else if (!strcmp(physics->upw_choice,_RUSANOV_)) {
         solver->Upwind    = NavierStokes3DUpwindRusanovModified;
         solver->UpwinddF  = NavierStokes3DUpwinddFRusanovModified;
         solver->UpwindFdF = NavierStokes3DUpwindFdFRusanovModified;
       } else {
         if (!mpi->rank) {
           fprintf(stderr,"Error in NavierStokes3DInitialize(): %s is not a valid upwinding scheme ",
                   physics->upw_choice);
           fprintf(stderr,"for use with split hyperbolic flux form. Use %s or %s.\n",
                   _ROE_,_RUSANOV_);
         }
         return(1);
       }
     } else {
       solver->JFunction      = NavierStokes3DJacobian;
       if      (!strcmp(physics->upw_choice,_ROE_    )) solver->Upwind = NavierStokes3DUpwindRoe;
       else if (!strcmp(physics->upw_choice,_RF_     )) solver->Upwind = NavierStokes3DUpwindRF;
       else if (!strcmp(physics->upw_choice,_LLF_    )) solver->Upwind = NavierStokes3DUpwindLLF;
       else if (!strcmp(physics->upw_choice,_RUSANOV_)) solver->Upwind = NavierStokes3DUpwindRusanov;
       else {
         if (!mpi->rank) {
           fprintf(stderr,"Error in NavierStokes3DInitialize(): %s is not a valid upwinding scheme. ",
                   physics->upw_choice);
           fprintf(stderr,"Choices are %s, %s, %s, and %s.\n",_ROE_,_RF_,_LLF_,_RUSANOV_);
         }
         return(1);
       }
     }
 
 #if defined(HAVE_CUDA)
   }
 #endif
 
   /* set the value of gamma in all the boundary objects */
   int n;
   DomainBoundary  *boundary = (DomainBoundary*) solver->boundary;
   for (n = 0; n < solver->nBoundaryZones; n++)  boundary[n].gamma = physics->gamma;
 
   /* finally, hijack the main solver's dissipation function pointer
    * to this model's own function, since it's difficult to express
    * the dissipation terms in the general form                      */
 #if defined(HAVE_CUDA)
   if (solver->use_gpu) {
     solver->ParabolicFunction = gpuNavierStokes3DParabolicFunction;
   } else {
 #endif
     solver->ParabolicFunction = NavierStokes3DParabolicFunction;
 #if defined(HAVE_CUDA)
   }
 #endif
 
   /* allocate array to hold the gravity field */
   int *dim    = solver->dim_local;
   int ghosts  = solver->ghosts;
   int d, size = 1; for (d=0; d<_MODEL_NDIMS_; d++) size *= (dim[d] + 2*ghosts);
   physics->grav_field_f = (double*) calloc (size, sizeof(double));
   physics->grav_field_g = (double*) calloc (size, sizeof(double));
   /* allocate arrays to hold the fast Jacobian for split form of the hyperbolic flux */
   physics->fast_jac     = (double*) calloc (_MODEL_NDIMS_*size*_MODEL_NVARS_*_MODEL_NVARS_,sizeof(double));
   physics->solution     = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
   /* initialize the gravity fields */
   IERR NavierStokes3DGravityField(solver,mpi); CHECKERR(ierr);
 
 #if defined(HAVE_CUDA)
   if (solver->use_gpu) gpuNavierStokes3DInitialize(s,m);
 #endif
 
   count++;
   return(0);
 }