NavierStokes2DFlux.c File Reference

Functions to compute the hyperbolic flux for 2D Navier-Stokes equations. More...

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

Go to the source code of this file.

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

Detailed Description

Functions to compute the hyperbolic flux for 2D Navier-Stokes equations.

Author

Debojyoti Ghosh

Definition in file NavierStokes2DFlux.c.

Function Documentation

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);
}