NavierStokes2DInitialize.c File Reference

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

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

Go to the source code of this file.

Functions
double	NavierStokes2DComputeCFL (void *, void *, double, double)
int	NavierStokes2DFlux (double *, double *, int, void *, double)
int	NavierStokes2DStiffFlux (double *, double *, int, void *, double)
int	NavierStokes2DNonStiffFlux (double *, double *, int, void *, double)
int	NavierStokes2DRoeAverage (double *, double *, double *, void *)
int	NavierStokes2DParabolicFunction (double *, double *, void *, void *, double)
int	NavierStokes2DSource (double *, double *, void *, void *, double)
int	NavierStokes2DJacobian (double *, double *, void *, int, int, int)
int	NavierStokes2DStiffJacobian (double *, double *, void *, int, int, int)
int	NavierStokes2DLeftEigenvectors (double *, double *, void *, int)
int	NavierStokes2DRightEigenvectors (double *, double *, void *, int)
int	NavierStokes2DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindRF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindSWFS (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindRusanov (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwinddFRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwinddFRF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwinddFLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwinddFRusanov (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindFdFRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindFdFRF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindFdFLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindFdFRusanov (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindRusanovModified (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwinddFRusanovModified (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DUpwindFdFRusanovModified (double *, double *, double *, double *, double *, double *, int, void *, double)
int	NavierStokes2DGravityField (void *, void *)
int	NavierStokes2DModifiedSolution (double *, double *, int, void *, void *, double)
int	NavierStokes2DPreStep (double *, void *, void *, double)
int	NavierStokes2DInitialize (void *s, void *m)

Detailed Description

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

Author

Debojyoti Ghosh

Definition in file NavierStokes2DInitialize.c.

Function Documentation

double NavierStokes2DComputeCFL	(	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 NavierStokes2DComputeCFL.c.

{
  HyPar   *solver = (HyPar*)   s;
  NavierStokes2D *param  = (NavierStokes2D*) solver->physics;
  _DECLARE_IERR_;

  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,e,P,c,dxinv,dyinv,local_cfl[2];
    _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,vx,vy,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 */

    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) */
    if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
    if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];

    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(max_cfl);
}

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

Compute the hyperbolic flux function for the 2D 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 \\ (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 \\ (e+p)v \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 or y) for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 21 of file NavierStokes2DFlux.c.

{
  HyPar           *solver = (HyPar*)   s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             *dim    = solver->dim_local;
  int             ghosts  = solver->ghosts;
  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 rho, vx, vy, e, P;
    _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,vx,vy,e,P,param->gamma);
    _NavierStokes2DSetFlux_((f+_MODEL_NVARS_*p),rho,vx,vy,e,P,dir);
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
  }

  return(0);
}

int NavierStokes2DStiffFlux	(	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 _NavierStokes2DSetStiffFlux_). 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 (NavierStokes2D::fast_jac) evaluated for the solution at the beginning of each time step ( \({\bf u}_{ref}\) is NavierStokes2D::solution). This is done in NavierStokes2DPreStep().

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 or y) for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 63 of file NavierStokes2DFlux.c.

{
  HyPar             *solver = (HyPar*)   s;
  NavierStokes2D    *param  = (NavierStokes2D*) 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+(2*p+dir)*JacSize;
    MatVecMult4(_MODEL_NVARS_,(f+_MODEL_NVARS_*p),Af,(u+_MODEL_NVARS_*p));
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
  }

  return(0);
}

int NavierStokes2DNonStiffFlux	(	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 _NavierStokes2DSetNonStiffFlux_). 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 NavierStokes2DFlux(), and \(A_f{\bf u}\) is the linearized stiff flux computed in NavierStokes2DStiffFlux().

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 or y) for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 104 of file NavierStokes2DFlux.c.

{
  HyPar             *solver = (HyPar*)   s;
  NavierStokes2D    *param  = (NavierStokes2D*) 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, e, P;
    _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,vx,vy,e,P,param->gamma);
    _NavierStokes2DSetFlux_(ftot,rho,vx,vy,e,P,dir);
    /* compute stiff stuff */
    double *Af = param->fast_jac+(2*p+dir)*JacSize;
    MatVecMult4(_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 NavierStokes2DRoeAverage	(	double *	uavg,
		double *	uL,
		double *	uR,
		void *	p
	)

Compute the Roe-averaged state for the 2D Navier Stokes equations. This function just calls the macro _NavierStokes2DRoeAverage_ and is not used by any functions within the 2D 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 NavierStokes2D with physics-related variables

Definition at line 15 of file NavierStokes2DFunctions.c.

{
  NavierStokes2D *param  = (NavierStokes2D*) p;
  _NavierStokes2DRoeAverage_(uavg,uL,uR,param->gamma);
  return(0);
}

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

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

\begin{equation} \frac {\partial} {\partial x} \left[\begin{array}{c} 0 \\ \tau_{xx} \\ \tau_{yx} \\ u \tau_{xx} + v \tau_{yx} - q_x \end{array}\right] + \frac {\partial} {\partial y} \left[\begin{array}{c} 0 \\ \tau_{xy} \\ \tau_{yy} \\ u \tau_{xy} + v \tau_{yy} - q_y \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}\right),\\ \tau_{xy} &= \left(\frac{\mu}{Re}\right)\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\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}\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} \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. NavierStokes2DInitialize() 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 38 of file NavierStokes2DParabolicFunction.c.

{
  HyPar           *solver   = (HyPar*) s;
  MPIVariables    *mpi      = (MPIVariables*) m;
  NavierStokes2D  *physics  = (NavierStokes2D*) solver->physics;
  int             i,j,v;
  _DECLARE_IERR_;

  int ghosts = solver->ghosts;
  int imax   = solver->dim_local[0];
  int jmax   = solver->dim_local[1];
  int *dim   = solver->dim_local;
  int nvars  = solver->nvars;
  int ndims  = solver->ndims;
  int size   = (imax+2*ghosts)*(jmax+2*ghosts)*nvars;

  _ArraySetValue_(par,size,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;

  double *Q; /* primitive variables */
  Q = (double*) calloc (size,sizeof(double));
  for (i=-ghosts; i<(imax+ghosts); i++) {
    for (j=-ghosts; j<(jmax+ghosts); j++) {
      int p,index[2]; index[0]=i; index[1]=j;
      double energy,pressure;
      _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
      _NavierStokes2DGetFlowVar_( (u+p),Q[p+0],Q[p+1],Q[p+2],energy,
                                  pressure,physics->gamma);
      Q[p+3] = physics->gamma*pressure/Q[p+0]; /* temperature */
    }
  }

  double *QDerivX = (double*) calloc (size,sizeof(double));
  double *QDerivY = (double*) calloc (size,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 MPIExchangeBoundariesnD(solver->ndims,solver->nvars,dim,
                               solver->ghosts,mpi,QDerivX);     CHECKERR(ierr);
  IERR MPIExchangeBoundariesnD(solver->ndims,solver->nvars,dim,
                               solver->ghosts,mpi,QDerivY);     CHECKERR(ierr);

  for (i=-ghosts; i<(imax+ghosts); i++) {
    for (j=-ghosts; j<(jmax+ghosts); j++) {
      int p,index[2]; index[0]=i; index[1]=j;
      double dxinv, dyinv;
      _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
      _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
      _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
      _ArrayScale1D_((QDerivX+p),dxinv,nvars);
      _ArrayScale1D_((QDerivY+p),dyinv,nvars);
    }
  }

  double *FViscous = (double*) calloc (size,sizeof(double));
  double *FDeriv   = (double*) calloc (size,sizeof(double));

  /* Along X */
  for (i=-ghosts; i<(imax+ghosts); i++) {
    for (j=-ghosts; j<(jmax+ghosts); j++) {
      int p,index[2]; index[0]=i; index[1]=j;;
      _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;

      double uvel, vvel, T, Tx, ux, uy, vx, vy;
      uvel = (Q+p)[1];
      vvel = (Q+p)[2];
      T    = (Q+p)[3];
      Tx   = (QDerivX+p)[3];
      ux   = (QDerivX+p)[1];
      vx   = (QDerivX+p)[2];
      uy   = (QDerivY+p)[1];
      vy   = (QDerivY+p)[2];

      /* calculate viscosity coeff based on Sutherland's law */
      double mu = raiseto(T, 0.76);

      double tau_xx, tau_xy, qx;
      tau_xx = two_third * (mu*inv_Re) * (2*ux - vy);
      tau_xy = (mu*inv_Re) * (uy + vx);
      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] = uvel*tau_xx + vvel*tau_xy + qx;
    }
  }
  IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
  for (i=0; i<imax; i++) {
    for (j=0; j<jmax; j++) {
      int p,index[2]; index[0]=i; index[1]=j;
      double dxinv;
      _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
      _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
      for (v=0; v<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++) {
      int p,index[2]; index[0]=i; index[1]=j;
      _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;

      double uvel, vvel, T, Ty, ux, uy, vx, vy;
      uvel = (Q+p)[1];
      vvel = (Q+p)[2];
      T    = (Q+p)[3];
      Ty   = (QDerivY+p)[3];
      ux   = (QDerivX+p)[1];
      vx   = (QDerivX+p)[2];
      uy   = (QDerivY+p)[1];
      vy   = (QDerivY+p)[2];

      /* calculate viscosity coeff based on Sutherland's law */
      double mu = raiseto(T, 0.76);

      double tau_yx, tau_yy, qy;
      tau_yx = (mu*inv_Re) * (uy + vx);
      tau_yy = two_third * (mu*inv_Re) * (-ux + 2*vy);
      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] = uvel*tau_yx + vvel*tau_yy + qy;
    }
  }
  IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
  for (i=0; i<imax; i++) {
    for (j=0; j<jmax; j++) {
      int p,index[2]; index[0]=i; index[1]=j;
      double dyinv;
      _ArrayIndex1D_(ndims,dim,index,ghosts,p); p *= nvars;
      _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
      for (v=0; v<nvars; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
    }
  }

  free(Q);
  free(QDerivX);
  free(QDerivY);
  free(FViscous);
  free(FDeriv);

  return(0);
}

int NavierStokes2DSource	(	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 u {\bf g}\cdot{\bf \hat{i}} - \rho v {\bf g}\cdot{\bf \hat{j}} \end{array}\right] = \left[\begin{array}{cccc} 0 & p_0 \varrho^{-1} & 0 & p_0 u \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial x}\left[\begin{array}{c} 0 \\ \varphi \\ 0 \\ \varphi \end{array}\right] + \left[\begin{array}{cccc} 0 & 0 & p_0 \varrho^{-1} & p_0 v \varrho^{-1} \end{array}\right] \cdot \frac{\partial}{\partial y}\left[\begin{array}{c} 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 NavierStokes2DGravityField(). 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 37 of file NavierStokes2DSource.c.

{
  HyPar           *solver = (HyPar* )         s;
  MPIVariables    *mpi    = (MPIVariables*)   m;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;

  if ((param->grav_x == 0.0) && (param->grav_y == 0.0))
    return(0); /* no gravitational forces */

  int     v, done, p, p1, p2;
  double  *SourceI = solver->fluxI; /* interace source term       */
  double  *SourceC = solver->fluxC; /* cell-centered source term  */
  double  *SourceL = solver->fL;
  double  *SourceR = solver->fR;

  int     ndims   = solver->ndims;
  int     ghosts  = solver->ghosts;
  int     *dim    = solver->dim_local;
  double  *x      = solver->x;
  double  *dxinv  = solver->dxinv;
  double  RT      =  param->p0 / param->rho0;
  int     index[ndims],index1[ndims],index2[ndims],dim_interface[ndims];

  /* Along X-direction */
  if (param->grav_x != 0.0) {
    /* set interface dimensions */
    _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;
    /* calculate the split source function exp(-phi/RT) */
    IERR NavierStokes2DSourceFunction(SourceC,u,x,solver,mpi,t,_XDIR_); CHECKERR(ierr);
    /* calculate the left and right interface source terms */
    IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
    IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
    /* calculate the final interface source term */
    IERR NavierStokes2DSourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
    /* calculate the final cell-centered source term */
    done = 0; _ArraySetValue_(index,ndims,0);
    while (!done) {
      _ArrayCopy1D_(index,index1,ndims);
      _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
      _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
      _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
      _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);

      double dx_inverse; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
      double rho, uvel, vvel, e, P; 
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,e,P,param->gamma);
      double term[_MODEL_NVARS_] = {0.0, rho*RT, 0.0, rho*RT*uvel};
      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 );
      }
      uvel = P; /* useless statement to avoid compiler warnings */
      _ArrayIncrementIndex_(ndims,dim,index,done);
    }
  }

  /* Along Y-direction */
  if (param->grav_y != 0.0) {
    /* set interface dimensions */
    _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_YDIR_]++;
    /* calculate the split source function exp(-phi/RT) */
    IERR NavierStokes2DSourceFunction(SourceC,u,x,solver,mpi,t,_YDIR_); CHECKERR(ierr);
    /* calculate the left and right interface source terms */
    IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
    IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
    /* calculate the final interface source term */
    IERR NavierStokes2DSourceUpwind(SourceI,SourceL,SourceR,u,_YDIR_,solver,t);

    /* calculate the final cell-centered source term */
    done = 0; _ArraySetValue_(index,ndims,0);
    while (!done) {
      _ArrayCopy1D_(index,index1,ndims);
      _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
      _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
      _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
      _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);

      double dy_inverse; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
      double rho, uvel, vvel, e, P; 
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,e,P,param->gamma);
      double term[_MODEL_NVARS_] = {0.0, 0.0, rho*RT, rho*RT*vvel};
      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])*dy_inverse );
      }
      uvel = P; /* useless statement to avoid compiler warnings */
      _ArrayIncrementIndex_(ndims,dim,index,done);
    }
  }

  return(0);
}

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

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

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 16
u	solution at a grid point (array of size nvar = 4)
p	object containing the physics-related parameters
dir	dimension (0 -> x, 1 -> y)
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 14 of file NavierStokes2DJacobian.c.

{
  NavierStokes2D *param = (NavierStokes2D*) 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 */
  _NavierStokes2DEigenvalues_      (u,D,param->gamma,dir);
  _NavierStokes2DLeftEigenvectors_ (u,L,param->gamma,dir);
  _NavierStokes2DRightEigenvectors_(u,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 = 5;  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 = 10; 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 = 15; 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]) );

  MatMult4(_MODEL_NVARS_,DL,D,L);
  MatMult4(_MODEL_NVARS_,Jac,R,DL);

  return(0);
}

int NavierStokes2DStiffJacobian	(	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 2D Navier-Stokes equations, given the solution at a grid point. The Jacobian is square matrix of size nvar=4, and is returned as a 1D array (double) of 16 elements in row-major format.

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 16
u	solution at a grid point (array of size nvar = 4)
p	object containing the physics-related parameters
dir	dimension (0 -> x, 1 -> y)
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 51 of file NavierStokes2DJacobian.c.

{
  NavierStokes2D *param = (NavierStokes2D*) 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 */
  _NavierStokes2DEigenvalues_      (u,D,param->gamma,dir);
  _NavierStokes2DLeftEigenvectors_ (u,L,param->gamma,dir);
  _NavierStokes2DRightEigenvectors_(u,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 = 5;  D[k] = ( dir == _YDIR_ ? 0.0 : absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) ) );
  k = 10; D[k] = ( dir == _XDIR_ ? 0.0 : absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) ) );
  k = 15; D[k] = 0.0;

  MatMult4(_MODEL_NVARS_,DL,D,L);
  MatMult4(_MODEL_NVARS_,Jac,R,DL);

  return(0);
}

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

Compute the left eigenvections for the 2D Navier Stokes equations. This function just calls the macro _NavierStokes2DLeftEigenvectors_ and is not used by any functions within the 2D 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 = 4^2 to save the matrix of left eigenvectors in (row-major format).
p	Object of type NavierStokes2D with physics-related variables
dir	Spatial dimension (x or y)

Definition at line 19 of file NavierStokes2DEigen.c.

{
  NavierStokes2D *param  = (NavierStokes2D*)  p;
  _NavierStokes2DLeftEigenvectors_(u,L,param->gamma,dir);
  return(0);
}

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

Compute the right eigenvections for the 2D Navier Stokes equations. This function just calls the macro _NavierStokes2DRightEigenvectors_ and is not used by any functions within the 2D 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 = 4^2 to save the matrix of right eigenvectors in (row-major format).
p	Object of type NavierStokes2D with physics-related variables
dir	Spatial dimension (x or y)

Definition at line 38 of file NavierStokes2DEigen.c.

{
  NavierStokes2D *param  = (NavierStokes2D*)  p;
  _NavierStokes2DRightEigenvectors_(u,R,param->gamma,dir);
  return(0);
}

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

Roe's upwinding scheme.

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

• 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 2D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 34 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             done;

  int *dim  = solver->dim_local;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; 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_];

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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]);

      _NavierStokes2DRoeAverage_        (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

       /* Harten's Entropy Fix - Page 362 of Leveque */
      int k;
      double delta = 0.000001, delta2 = delta*delta;
      double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
      k=0;  D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=5;  D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=10; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=15; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );

      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_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];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

Characteristic-based Roe-fixed upwinding scheme.

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

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

• 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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 120 of file NavierStokes2DUpwind.c.

{
  HyPar    *solver = (HyPar*)    s;
  NavierStokes2D  *param  = (NavierStokes2D*)  solver->physics;
  int      done,k;

  int *dim  = solver->dim_local;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; int index_outer[2] = {0,0}, index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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 */

      _NavierStokes2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);

      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_(uavg,R,param->gamma,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      double eigL[4],eigC[4],eigR[4];
      _NavierStokes2DEigenvalues_((uL+_MODEL_NVARS_*p),D,param->gamma,dir);
      eigL[0] = D[0];
      eigL[1] = D[5];
      eigL[2] = D[10];
      eigL[3] = D[15];
      _NavierStokes2DEigenvalues_((uR+_MODEL_NVARS_*p),D,param->gamma,dir);
      eigR[0] = D[0];
      eigR[1] = D[5];
      eigR[2] = D[10];
      eigR[3] = D[15];
      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      eigC[0] = D[0];
      eigC[1] = D[5];
      eigC[2] = D[10];
      eigC[3] = D[15];

      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 */
      MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

Characteristic-based local Lax-Friedrich upwinding scheme.

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

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

• 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 2D Navier Stokes equations when there is a non-zero gravitational force. See the reference below. For flows without any gravitational forces, it reduces to its original form.

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 223 of file NavierStokes2DUpwind.c.

{
  HyPar    *solver = (HyPar*)    s;
  NavierStokes2D  *param  = (NavierStokes2D*)  solver->physics;
  int      done;

  int *dim  = solver->dim_local;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; int index_outer[2] = {0,0}, index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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 uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);

      /* Local Lax-Friedrich upwinding scheme */

      _NavierStokes2DRoeAverage_        (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      double eigL[4],eigC[4],eigR[4];
      _NavierStokes2DEigenvalues_((u+_MODEL_NVARS_*pL),D,param->gamma,dir);
      eigL[0] = D[0];
      eigL[1] = D[5];
      eigL[2] = D[10];
      eigL[3] = D[15];
      _NavierStokes2DEigenvalues_((u+_MODEL_NVARS_*pR),D,param->gamma,dir);
      eigR[0] = D[0];
      eigR[1] = D[5];
      eigR[2] = D[10];
      eigR[3] = D[15];
      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      eigC[0] = D[0];
      eigC[1] = D[5];
      eigC[2] = D[10];
      eigC[3] = D[15];

      double alpha;
      alpha = kappa * 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]));

      /* calculate the interface flux from the characteristic flux */
      MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

Steger-Warming Flux-Splitting scheme

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

Note that this method cannot be used for flows with non-zero gravitational forces.

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 316 of file NavierStokes2DUpwind.c.

{
  HyPar             *solver = (HyPar*)    s;
  NavierStokes2D    *param  = (NavierStokes2D*)  solver->physics;
  int               done,k;
  _DECLARE_IERR_;

  int ndims = solver->ndims;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double fp[_MODEL_NVARS_], fm[_MODEL_NVARS_],uavg[_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      double rho,vx,vy,e,P,c,gamma=param->gamma,term,Mach,lp[_MODEL_NVARS_],lm[_MODEL_NVARS_];

      /* Steger Warming flux splitting */
      _NavierStokes2DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param->gamma);
      _NavierStokes2DGetFlowVar_(uavg,rho,vx,vy,e,P,param->gamma);
      Mach = (dir==_XDIR_ ? vx : vy) / sqrt(gamma*P/rho);

      if (Mach < -1.0) {

        _ArrayCopy1D_((fR+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);

      } else if (Mach < 1.0) {

        double kx = 0, ky = 0;
        kx = (dir==_XDIR_ ? 1.0 : 0.0);
        ky = (dir==_YDIR_ ? 1.0 : 0.0);

        _NavierStokes2DGetFlowVar_((uL+_MODEL_NVARS_*p),rho,vx,vy,e,P,param->gamma);
        c = sqrt(gamma*P/rho);
        term = rho/(2.0*gamma);
        lp[0] = lp[1] = kx*vx + ky*vy;
        lp[2] = lp[0] + c;
        lp[3] = lp[0] - c;
        for (k=0; k<_MODEL_NVARS_; k++) if (lp[k] < 0.0) lp[k] = 0.0;

        fp[0] = term * (2.0*(gamma-1.0)*lp[0] + lp[2] + lp[3]);
        fp[1] = term * (2.0*(gamma-1.0)*lp[0]*vx + lp[2]*(vx+c*kx) + lp[3]*(vx-c*kx));
        fp[2] = term * (2.0*(gamma-1.0)*lp[0]*vy + lp[2]*(vy+c*ky) + lp[3]*(vy-c*ky));
        fp[3] = term * ((gamma-1.0)*lp[0]*(vx*vx+vy*vy) + 0.5*lp[2]*((vx+c*kx)*(vx+c*kx) + (vy+c*ky)*(vy+c*ky))
                        + 0.5*lp[3]*((vx-c*kx)*(vx-c*kx) + (vy-c*ky)*(vy-c*ky))
                        + ((3.0-gamma)*(lp[2]+lp[3])*c*c)/(2.0*(gamma-1.0)) );

        _NavierStokes2DGetFlowVar_((uR+_MODEL_NVARS_*p),rho,vx,vy,e,P,param->gamma);
        c = sqrt(gamma*P/rho);
        term = rho/(2.0*gamma);
        lm[0] = lm[1] = kx*vx + ky*vy;
        lm[2] = lm[0] + c;
        lm[3] = lm[0] - c;
        for (k=0; k<_MODEL_NVARS_; k++) if (lm[k] > 0.0) lm[k] = 0.0;

        fm[0] = term * (2.0*(gamma-1.0)*lm[0] + lm[2] + lm[3]);
        fm[1] = term * (2.0*(gamma-1.0)*lm[0]*vx + lm[2]*(vx+c*kx) + lm[3]*(vx-c*kx));
        fm[2] = term * (2.0*(gamma-1.0)*lm[0]*vy + lm[2]*(vy+c*ky) + lm[3]*(vy-c*ky));
        fm[3] = term * ((gamma-1.0)*lm[0]*(vx*vx+vy*vy) + 0.5*lm[2]*((vx+c*kx)*(vx+c*kx) + (vy+c*ky)*(vy+c*ky))
                        + 0.5*lm[3]*((vx-c*kx)*(vx-c*kx) + (vy-c*ky)*(vy-c*ky))
                        + ((3.0-gamma)*(lm[2]+lm[3])*c*c)/(2.0*(gamma-1.0)) );

        _ArrayAdd1D_((fI+_MODEL_NVARS_*p),fp,fm,_MODEL_NVARS_);

      } else {

        _ArrayCopy1D_((fL+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);

      }

    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

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

Rusanov's upwinding scheme.

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

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

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 428 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             *dim    = solver->dim_local, done;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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_];

      /* Modified Rusanov'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]);

      _NavierStokes2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);

      double c, vel[_MODEL_NDIMS_], rho,E,P;
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaL = c + absolute(vel[dir]), betaL = absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaR = c + absolute(vel[dir]), betaR = absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaavg = c + absolute(vel[dir]), betaavg = absolute(vel[dir]);

      double kappa  = max(param->grav_field_g[pL],param->grav_field_g[pR]);
      double alpha  = kappa*max3(alphaL,alphaR,alphaavg);

      fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
      fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
      fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
      fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 589 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             done;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; 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_];

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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]);

      _NavierStokes2DRoeAverage_        (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

       /* Harten's Entropy Fix - Page 362 of Leveque */
      int k;
      double delta = 0.000001, delta2 = delta*delta;
      double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
      k=0;  D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=5;  D[k] = (dir == _YDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
      k=10; D[k] = (dir == _XDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
      k=15; D[k] = 0.0;

      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_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];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 670 of file NavierStokes2DUpwind.c.

{
  HyPar    *solver = (HyPar*)    s;
  NavierStokes2D  *param  = (NavierStokes2D*)  solver->physics;
  int      done,k;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; int index_outer[2] = {0,0}, index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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 uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];

      /* Roe-Fixed upwinding scheme */

      _NavierStokes2DRoeAverage_(uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);

      _NavierStokes2DEigenvalues_      (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_(uavg,R,param->gamma,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      double eigL[4],eigC[4],eigR[4];
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
      eigL[0] = D[0];
      eigL[1] = (dir == _YDIR_ ? 0.0 : D[5]);
      eigL[2] = (dir == _XDIR_ ? 0.0 : D[10]);
      eigL[3] = 0.0;
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
      eigR[0] = D[0];
      eigR[1] = (dir == _YDIR_ ? 0.0 : D[5]);
      eigR[2] = (dir == _XDIR_ ? 0.0 : D[10]);
      eigR[3] = 0.0;
      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      eigC[0] = D[0];
      eigC[1] = (dir == _YDIR_ ? 0.0 : D[5]);
      eigC[2] = (dir == _XDIR_ ? 0.0 : D[10]);
      eigC[3] = 0.0;

      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 */
      MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

The characteristic-based local Lax-Friedrich upwinding scheme (NavierStokes2DUpwindLLF) for the partitioned hyperbolic flux that comprises of the acoustic waves only (see NavierStokes2DStiffFlux, _NavierStokes2DSetStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \( u\pm a\) are used. Reference:

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 764 of file NavierStokes2DUpwind.c.

{
  HyPar    *solver = (HyPar*)    s;
  NavierStokes2D  *param  = (NavierStokes2D*)  solver->physics;
  int      done;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; int index_outer[2] = {0,0}, index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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 uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);

      /* Local Lax-Friedrich upwinding scheme */

      _NavierStokes2DRoeAverage_        (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      double eigL[4],eigC[4],eigR[4];
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
      eigL[0] = D[0];
      eigL[1] = (dir == _YDIR_ ? 0.0 : D[5]);
      eigL[2] = (dir == _XDIR_ ? 0.0 : D[10]);
      eigL[3] = 0.0;
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
      eigR[0] = D[0];
      eigR[1] = (dir == _YDIR_ ? 0.0 : D[5]);
      eigR[2] = (dir == _XDIR_ ? 0.0 : D[10]);
      eigR[3] = 0.0;
      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      eigC[0] = D[0];
      eigC[1] = (dir == _YDIR_ ? 0.0 : D[5]);
      eigC[2] = (dir == _XDIR_ ? 0.0 : D[10]);
      eigC[3] = 0.0;

      double alpha;
      alpha = kappa * 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]));

      /* calculate the interface flux from the characteristic flux */
      MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

int NavierStokes2DUpwinddFRusanov	(	double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		int	,
		void *	,
		double
	)

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 950 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             done;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; 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_];

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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;
      double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);

      /* 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]);

      /* Compute total dissipation */
      _NavierStokes2DRoeAverage_        (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
      k=0;  D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=5;  D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=10; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=15; D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_MODEL_NVARS_,udiss_total,modA,udiff);

      /* Compute dissipation corresponding to acoustic modes */
      _NavierStokes2DRoeAverage_        (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
      k=0;  D[k] = kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
      k=5;  D[k] = (dir == _YDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
      k=10; D[k] = (dir == _XDIR_ ? 0.0 : kappa * (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) ) );
      k=15; D[k] = 0.0;
      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_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];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 1046 of file NavierStokes2DUpwind.c.

{
  HyPar    *solver = (HyPar*)    s;
  NavierStokes2D  *param  = (NavierStokes2D*)  solver->physics;
  int      done,k;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; int index_outer[2] = {0,0}, index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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 uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];

      /* Roe-Fixed upwinding scheme */

      _NavierStokes2DRoeAverage_(uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);

      _NavierStokes2DEigenvalues_      (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_ (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_(uavg,R,param->gamma,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      double eigL[4],eigC[4],eigR[4];
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
      eigL[0] = 0.0;
      eigL[1] = (dir == _XDIR_ ? 0.0 : D[5]);
      eigL[2] = (dir == _YDIR_ ? 0.0 : D[10]);
      eigL[3] = D[15];
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
      eigR[0] = 0.0;
      eigR[1] = (dir == _XDIR_ ? 0.0 : D[5]);
      eigR[2] = (dir == _YDIR_ ? 0.0 : D[10]);
      eigR[3] = D[15];
      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      eigC[0] = 0.0;
      eigC[1] = (dir == _XDIR_ ? 0.0 : D[5]);
      eigC[2] = (dir == _YDIR_ ? 0.0 : D[10]);
      eigC[3] = D[15];

      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 */
      MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

The characteristic-based local Lax-Friedrich upwinding scheme (NavierStokes2DUpwindLLF) for the partitioned hyperbolic flux that comprises of the entropy waves only (see NavierStokes2DNonStiffFlux, _NavierStokes2DSetNonStiffFlux_). Thus, only the characteristic fields / eigen-modes corresponding to \(u\) are used. Reference:

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 1140 of file NavierStokes2DUpwind.c.

{
  HyPar    *solver = (HyPar*)    s;
  NavierStokes2D  *param  = (NavierStokes2D*)  solver->physics;
  int      done;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[dir]++;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_],
                L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; int index_outer[2] = {0,0}, index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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 uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
      double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);

      /* Local Lax-Friedrich upwinding scheme */

      _NavierStokes2DRoeAverage_        (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DEigenvalues_       (uavg,D,param->gamma,dir);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult4(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult4(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      double eigL[4],eigC[4],eigR[4];
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pL),D,param->gamma,dir);
      eigL[0] = 0.0;
      eigL[1] = (dir == _XDIR_ ? 0.0 : D[5]);
      eigL[2] = (dir == _YDIR_ ? 0.0 : D[10]);
      eigL[3] = D[15];
      _NavierStokes2DEigenvalues_((uref+_MODEL_NVARS_*pR),D,param->gamma,dir);
      eigR[0] = 0.0;
      eigR[1] = (dir == _XDIR_ ? 0.0 : D[5]);
      eigR[2] = (dir == _YDIR_ ? 0.0 : D[10]);
      eigR[3] = D[15];
      _NavierStokes2DEigenvalues_(uavg,D,param->gamma,dir);
      eigC[0] = 0.0;
      eigC[1] = (dir == _XDIR_ ? 0.0 : D[5]);
      eigC[2] = (dir == _YDIR_ ? 0.0 : D[10]);
      eigC[3] = D[15];

      double alpha;
      alpha = kappa * 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]));

      /* calculate the interface flux from the characteristic flux */
      MatVecMult4(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

int NavierStokes2DUpwindFdFRusanov	(	double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		int	,
		void *	,
		double
	)

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 499 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             done;

  int *dim  = solver->dim_local;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; 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_];

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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_];

      /* Modified Rusanov'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]);

      _NavierStokes2DRoeAverage_        (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

      double c, vel[_MODEL_NDIMS_], rho,E,P;
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaL = c + absolute(vel[dir]), betaL = absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaR = c + absolute(vel[dir]), betaR = absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaavg = c + absolute(vel[dir]), 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);
      D[0]  = alpha;
      D[5]  = (dir == _XDIR_ ? alpha : beta);
      D[10] = (dir == _YDIR_ ? alpha : beta);
      D[15] = beta;
      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_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];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 859 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             done;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; 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_];

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_];

      /* Rusanov'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]);

      _NavierStokes2DRoeAverage_        (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);

      double c, vel[_MODEL_NDIMS_], rho,E,P;
      _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaL = c + absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      double alphaR = c + absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],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);
      D[0]  = alpha;
      D[5]  = (dir == _YDIR_ ? 0.0 : alpha);
      D[10] = (dir == _XDIR_ ? 0.0 : alpha);
      D[15] = 0.0;

      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_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];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

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

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

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

Parameters

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 or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 1235 of file NavierStokes2DUpwind.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) solver->physics;
  int             done;

  int     *dim  = solver->dim_local;
  double  *uref = param->solution;

  int bounds_outer[2], bounds_inter[2];
  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[dir] = 1;
  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; 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_];

  done = 0; int index_outer[2] = {0,0}; int index_inter[2];
  while (!done) {
    index_inter[0] = index_outer[0]; index_inter[1] = index_outer[1];
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D2_(_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_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_],
             udiss_total[_MODEL_NVARS_],udiss_acoustic[_MODEL_NVARS_];

      /* 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[_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]);

      /* Compute total dissipation */
      _NavierStokes2DRoeAverage_        (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      alphaL = c + absolute(vel[dir]);
      betaL = absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      alphaR = c + absolute(vel[dir]);
      betaR = absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      alphaavg = c + absolute(vel[dir]);
      betaavg = absolute(vel[dir]);
      alpha  = kappa*max3(alphaL,alphaR,alphaavg);
      beta   = kappa*max3(betaL,betaR,betaavg);
      _ArraySetValue_(D,_MODEL_NVARS_*_MODEL_NVARS_,0.0);
      D[0]  = alpha;
      D[5]  = (dir == _XDIR_ ? alpha : beta);
      D[10] = (dir == _YDIR_ ? alpha : beta);
      D[15] = beta;
      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_MODEL_NVARS_,udiss_total,modA,udiff);

      /* Compute dissipation for the linearized acoustic modes */
      _NavierStokes2DRoeAverage_        (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param->gamma);
      _NavierStokes2DLeftEigenvectors_  (uavg,L,param->gamma,dir);
      _NavierStokes2DRightEigenvectors_ (uavg,R,param->gamma,dir);
      _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pL),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      alphaL = c + absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_((uref+_MODEL_NVARS_*pR),rho,vel[0],vel[1],E,P,param->gamma);
      c = sqrt(param->gamma*P/rho);
      alphaR = c + absolute(vel[dir]);
      _NavierStokes2DGetFlowVar_(uavg,rho,vel[0],vel[1],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);
      D[0]  = alpha;
      D[5]  = (dir == _YDIR_ ? 0.0 : alpha);
      D[10] = (dir == _XDIR_ ? 0.0 : alpha);
      D[15] = 0.0;
      MatMult4(_MODEL_NVARS_,DL,D,L);
      MatMult4(_MODEL_NVARS_,modA,R,DL);
      MatVecMult4(_MODEL_NVARS_,udiss_acoustic,modA,udiff);

      /* Compute dissipation for the entropy modes */
      _ArraySubtract1D_(udiss,udiss_total,udiss_acoustic,_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];
    }
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
  }

  return(0);
}

int NavierStokes2DGravityField	(	void *	s,
		void *	m
	)

This function computes the pressure and density variation functions for the hydrostatic balance of the type specified by the user (NavierStokes2D: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 (NavierStokes2D::rho0) and pressure (NavierStokes2D::p0). This function computes \(\varrho\) (NavierStokes2D::grav_field_f) and \(\varphi\) (NavierStokes2D::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 NavierStokes2DGravityField.c.

{
  HyPar           *solver = (HyPar*)          s;
  MPIVariables    *mpi    = (MPIVariables*)   m;
  NavierStokes2D  *param  = (NavierStokes2D*) 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 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);
      f[p] = exp( (gx*xcoord+gy*ycoord)/RT);
      g[p] = exp(-(gx*xcoord+gy*ycoord)/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);
      f[p] = raiseto((1.0-(gx*xcoord+gy*ycoord)/(Cp*T0)), (-1.0 /(gamma-1.0)));
      g[p] = raiseto((1.0-(gx*xcoord+gy*ycoord)/(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 ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
      if (gx != 0) {
        if (!mpi->rank) {
          fprintf(stderr,"Error in NavierStokes2DGravityField(): HB = 3 is implemented only for ");
          fprintf(stderr,"gravity force along the y-coordinate.\n");
        }
        return(1);
      }
      double Pexner = 1 + ((gy*gy)/(Cp*T0*Nbv*Nbv)) * (exp(-(Nbv*Nbv/gy)*ycoord)-1.0);
      f[p] = raiseto(Pexner, (-1.0 /(gamma-1.0))) * exp(Nbv*Nbv*ycoord/gy);
      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 NavierStokes2DModifiedSolution	(	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) \\ 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\right) \end{equation}

\(\varrho\) and \(\varphi\) are computed in NavierStokes2DGravityField(). 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 NavierStokes2DModifiedSolution.c.

{
  HyPar           *solver = (HyPar*)          s;
  NavierStokes2D  *param  = (NavierStokes2D*) 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, E, P;
    _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,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] = (P*inv_gamma_m1)*(1.0/param->grav_field_g[p]) + (0.5*rho*(uvel*uvel+vvel*vvel))*param->grav_field_f[p];
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

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

Pre-step function for the 2D 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 NavierStokes2D::solution for the linearized flux partitioning.
• The linearized "fast" Jacobian (representing the acoustic modes) NavierStokes2D::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,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 33 of file NavierStokes2DPreStep.c.

{
  HyPar             *solver = (HyPar*)   s;
  NavierStokes2D    *param  = (NavierStokes2D*) 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);

    dir = _XDIR_;
    A = (param->fast_jac + 2*JacSize*p + dir*JacSize);
    /* get the eigenvalues, and left & right eigenvectors */
    _NavierStokes2DEigenvalues_      ((u+_MODEL_NVARS_*p),D,param->gamma,dir);
    _NavierStokes2DLeftEigenvectors_ ((u+_MODEL_NVARS_*p),L,param->gamma,dir);
    _NavierStokes2DRightEigenvectors_((u+_MODEL_NVARS_*p),R,param->gamma,dir);
    /* remove the entropy modes (corresponding to eigenvalues u) */
    D[2*_MODEL_NVARS_+2] = D[3*_MODEL_NVARS_+3] = 0.0;
    /* assemble the Jacobian */
    MatMult4(_MODEL_NVARS_,DL,D,L );
    MatMult4(_MODEL_NVARS_,A ,R,DL);

    dir = _YDIR_;
    A = (param->fast_jac + 2*JacSize*p + dir*JacSize);
    /* get the eigenvalues, and left & right eigenvectors */
    _NavierStokes2DEigenvalues_      ((u+_MODEL_NVARS_*p),D,param->gamma,dir)
    _NavierStokes2DLeftEigenvectors_ ((u+_MODEL_NVARS_*p),L,param->gamma,dir);
    _NavierStokes2DRightEigenvectors_((u+_MODEL_NVARS_*p),R,param->gamma,dir);
    /* remove the entropy modes (corresponding to eigenvalues v) */
    D[1*_MODEL_NVARS_+1] = D[3*_MODEL_NVARS_+3] = 0.0;
    /* assemble the Jacobian */
    MatMult4(_MODEL_NVARS_,DL,D,L );
    MatMult4(_MODEL_NVARS_,A ,R,DL);

    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

int NavierStokes2DInitialize	(	void *	s,
		void *	m
	)

Initialize the 2D Navier-Stokes (NavierStokes2D) 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	NavierStokes2D::gamma	1.4
Pr	double	NavierStokes2D::Pr	0.72
Re	double	NavierStokes2D::Re	-1
Minf	double	NavierStokes2D::Minf	1.0
gravity	double,double	NavierStokes2D::grav_x,NavierStokes2D::grav_y	0.0,0.0
rho_ref	double	NavierStokes2D::rho0	1.0
p_ref	double	NavierStokes2D::p0	1.0
HB	int	NavierStokes2D::HB	1
R	double	NavierStokes2D::R	1.0
upwinding	char[]	NavierStokes2D::upw_choice	"roe" (_ROE_)

• If "HB" (NavierStokes2D::HB) is specified as 3, it should be followed by the the Brunt-Vaisala frequency (NavierStokes2D::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 103 of file NavierStokes2DInitialize.c.

{
  HyPar           *solver  = (HyPar*)          s;
  MPIVariables    *mpi     = (MPIVariables*)   m;
  NavierStokes2D  *physics = (NavierStokes2D*) solver->physics;
  int             ferr     = 0;

  static int count = 0;

  if (solver->nvars != _MODEL_NVARS_) {
    fprintf(stderr,"Error in NavierStokes2DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
    return(1);
  }
  if (solver->ndims != _MODEL_NDIMS_) {
    fprintf(stderr,"Error in NavierStokes2DInitialize(): 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->rho0   = 1.0;
  physics->p0     = 1.0;
  physics->HB     = 1;
  physics->R      = 1.0;
  physics->N_bv   = 0.0;
  strcpy(physics->upw_choice,"roe");

  /* 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) return(1);
      if (!strcmp(word, "begin")){
          while (strcmp(word, "end")){
              ferr = fscanf(in,"%s",word);                  if (ferr != 1) return(1);
          if (!strcmp(word, "gamma")) {
            ferr = fscanf(in,"%lf",&physics->gamma);    if (ferr != 1) return(1);
          } else if (!strcmp(word,"upwinding")) {
            ferr = fscanf(in,"%s",physics->upw_choice); if (ferr != 1) return(1);
          } else if (!strcmp(word,"Pr")) {
            ferr = fscanf(in,"%lf",&physics->Pr);       if (ferr != 1) return(1);
          } else if (!strcmp(word,"Re")) {
            ferr = fscanf(in,"%lf",&physics->Re);       if (ferr != 1) return(1);
          } else if (!strcmp(word,"Minf")) {
            ferr = fscanf(in,"%lf",&physics->Minf);     if (ferr != 1) return(1);
          } else if (!strcmp(word,"gravity")) {
            ferr = fscanf(in,"%lf",&physics->grav_x);   if (ferr != 1) return(1);
            ferr = fscanf(in,"%lf",&physics->grav_y);   if (ferr != 1) return(1);
          } else if (!strcmp(word,"rho_ref")) {
            ferr = fscanf(in,"%lf",&physics->rho0);     if (ferr != 1) return(1);
          } else if (!strcmp(word,"p_ref")) {
            ferr = fscanf(in,"%lf",&physics->p0);       if (ferr != 1) return(1);
          } else if (!strcmp(word,"HB")) {
            ferr = fscanf(in,"%d",&physics->HB);        if (ferr != 1) return(1);
            if (physics->HB==3) {
              ferr = fscanf(in,"%lf",&physics->N_bv);   if (ferr != 1) return(1);
            }
          } else if (!strcmp(word,"R")) {
            ferr = fscanf(in,"%lf",&physics->R);        if (ferr != 1) return(1);
          } 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_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->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_integer   (&physics->HB       ,1                ,0,&mpi->world); CHECKERR(ierr);

  /* Scaling the Reynolds number with the 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))
      && (strcmp(physics->upw_choice,_LLF_    ))
      && (strcmp(physics->upw_choice,_RUSANOV_))
      && (strcmp(physics->upw_choice,_ROE_    ))              ) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in NavierStokes2DInitialize: %s, %s or %s upwinding is needed for flows ",_LLF_,_ROE_,_RUSANOV_);
      fprintf(stderr,"with gravitational forces.\n");
    }
    return(1);
  }
  /* check that solver has the correct choice of diffusion formulation */
  if (strcmp(solver->spatial_type_par,_NC_2STAGE_) && (physics->Re > 0)) {
    if (!mpi->rank)
      fprintf(stderr,"Error in NavierStokes2DInitialize(): Parabolic term spatial discretization must be \"%s\"\n",_NC_2STAGE_);
    return(1);
  }

  /* initializing physical model-specific functions */
#if defined(HAVE_CUDA) && defined(CUDA_VAR_ORDERDING_AOS)
  if (solver->use_gpu) {
    solver->FFunction  = gpuNavierStokes2DFlux;
    solver->SFunction  = gpuNavierStokes2DSource;
    solver->UFunction  = gpuNavierStokes2DModifiedSolution;
  } else {
#endif
    solver->PreStep               = NavierStokes2DPreStep;
    solver->ComputeCFL            = NavierStokes2DComputeCFL;
    solver->FFunction             = NavierStokes2DFlux;
    solver->SFunction             = NavierStokes2DSource;
    solver->UFunction             = NavierStokes2DModifiedSolution;
    solver->AveragingFunction     = NavierStokes2DRoeAverage;
    solver->GetLeftEigenvectors   = NavierStokes2DLeftEigenvectors;
    solver->GetRightEigenvectors  = NavierStokes2DRightEigenvectors;
#if defined(HAVE_CUDA) && defined(CUDA_VAR_ORDERDING_AOS)
  }
#endif

#if defined(HAVE_CUDA) && defined(CUDA_VAR_ORDERDING_AOS)
  if (solver->use_gpu) {
    if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
      if (!mpi->rank) {
        fprintf(stderr,"Error in NavierStokes2DInitialize(): Not yet implemented on GPU!");
      }
      return 1;
    } else {
      solver->JFunction = NavierStokes2DJacobian;
      if (!strcmp(physics->upw_choice,_RUSANOV_)) solver->Upwind = gpuNavierStokes2DUpwindRusanov;
      else {
        if (!mpi->rank) {
          fprintf(stderr,"Error in NavierStokes2DInitialize(): %s is not implemented on GPU. ",
                  physics->upw_choice);
          fprintf(stderr,"Only choice is %s.\n",_RUSANOV_);
        }
        return 1;
      }
    }
  } else {
#endif
    if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
      solver->FdFFunction = NavierStokes2DNonStiffFlux;
      solver->dFFunction  = NavierStokes2DStiffFlux;
      solver->JFunction = NavierStokes2DStiffJacobian;
      if (!strcmp(physics->upw_choice,_ROE_)) {
        solver->Upwind    = NavierStokes2DUpwindRoe;
        solver->UpwinddF  = NavierStokes2DUpwinddFRoe;
        solver->UpwindFdF = NavierStokes2DUpwindFdFRoe;
      } else if (!strcmp(physics->upw_choice,_RF_)) {
        solver->Upwind    = NavierStokes2DUpwindRF;
        solver->UpwinddF  = NavierStokes2DUpwinddFRF;
        solver->UpwindFdF = NavierStokes2DUpwindFdFRF;
      } else if (!strcmp(physics->upw_choice,_LLF_)) {
        solver->Upwind    = NavierStokes2DUpwindLLF;
        solver->UpwinddF  = NavierStokes2DUpwinddFLLF;
        solver->UpwindFdF = NavierStokes2DUpwindFdFLLF;
      } else if (!strcmp(physics->upw_choice,_RUSANOV_)) {
        solver->Upwind    = NavierStokes2DUpwindRusanovModified;
        solver->UpwinddF  = NavierStokes2DUpwinddFRusanovModified;
        solver->UpwindFdF = NavierStokes2DUpwindFdFRusanovModified;
      } else {
        if (!mpi->rank) {
          fprintf(stderr,"Error in NavierStokes2DInitialize(): %s is not a valid upwinding scheme ",
                  physics->upw_choice);
          fprintf(stderr,"for use with split hyperbolic flux form. Use %s, %s, %s, or %s.\n",
                  _ROE_,_RF_,_LLF_,_RUSANOV_);
        }
        return(1);
      }
    } else {
      solver->JFunction = NavierStokes2DJacobian;
      if      (!strcmp(physics->upw_choice,_ROE_    )) solver->Upwind = NavierStokes2DUpwindRoe;
      else if (!strcmp(physics->upw_choice,_RF_     )) solver->Upwind = NavierStokes2DUpwindRF;
      else if (!strcmp(physics->upw_choice,_LLF_    )) solver->Upwind = NavierStokes2DUpwindLLF;
      else if (!strcmp(physics->upw_choice,_SWFS_   )) solver->Upwind = NavierStokes2DUpwindSWFS;
      else if (!strcmp(physics->upw_choice,_RUSANOV_)) solver->Upwind = NavierStokes2DUpwindRusanov;
      else {
        if (!mpi->rank) {
          fprintf(stderr,"Error in NavierStokes2DInitialize(): %s is not a valid upwinding scheme. ",
                  physics->upw_choice);
          fprintf(stderr,"Choices are %s, %s, %s, %s, and %s.\n",_ROE_,_RF_,_LLF_,_SWFS_,_RUSANOV_);
        }
        return(1);
      }
    }
#if defined(HAVE_CUDA) && defined(CUDA_VAR_ORDERDING_AOS)
  }
#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;

  /* 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) && defined(CUDA_VAR_ORDERDING_AOS)
  if (solver->use_gpu) {
    solver->ParabolicFunction = NavierStokes2DParabolicFunction;
  } else {
#endif
    solver->ParabolicFunction = NavierStokes2DParabolicFunction;
#if defined(HAVE_CUDA) && defined(CUDA_VAR_ORDERDING_AOS)
  }
#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 (2*size*_MODEL_NVARS_*_MODEL_NVARS_,sizeof(double));
  physics->solution     = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
  /* initialize the gravity fields */
  IERR NavierStokes2DGravityField(solver,mpi); CHECKERR(ierr);

#if defined(HAVE_CUDA) && defined(CUDA_VAR_ORDERDING_AOS)
  if (solver->use_gpu) gpuNavierStokes2DInitialize(s,m);
#endif

  count++;
  return(0);
}