navierstokes2d.h File Reference

2D Navier Stokes equations (compressible flows) More...

#include <basic.h>

Go to the source code of this file.

Data Structures
struct	NavierStokes2D
	Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations. More...

Macros
#define	_NAVIER_STOKES_2D_   "navierstokes2d"
#define	_MODEL_NDIMS_   2
#define	_MODEL_NVARS_   4
#define	_ROE_   "roe"
#define	_RF_   "rf-char"
#define	_LLF_   "llf-char"
#define	_SWFS_   "steger-warming"
#define	_RUSANOV_   "rusanov"
#define	_XDIR_   0
#define	_YDIR_   1
#define	_NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
#define	_NavierStokes2DSetFlux_(f, rho, vx, vy, e, P, dir)
#define	_NavierStokes2DSetStiffFlux_(f, rho, vx, vy, e, P, dir, gamma)
#define	_NavierStokes2DSetNonStiffFlux_(f, rho, vx, vy, e, P, dir, gamma)
#define	_NavierStokes2DRoeAverage_(uavg, uL, uR, gamma)
#define	_NavierStokes2DEigenvalues_(u, D, gamma, dir)
#define	_NavierStokes2DLeftEigenvectors_(u, L, ga, dir)
#define	_NavierStokes2DRightEigenvectors_(u, R, ga, dir)

Functions
int	NavierStokes2DInitialize (void *, void *)
int	NavierStokes2DCleanup (void *)

Detailed Description

2D Navier Stokes equations (compressible flows)

Author

Debojyoti Ghosh

2D Navier-Stokes equations for viscous and inviscid compressible flows (with gravitational terms)

\begin{equation} \frac {\partial} {\partial t} \left[\begin{array}{c} \rho \\ \rho u \\ \rho v \\ e \end{array}\right] + \frac {\partial} {\partial x} \left[\begin{array}{c} \rho u \\ \rho u^2 + p \\ \rho u v \\ (e+p) u\end{array}\right] + \frac {\partial} {\partial y} \left[\begin{array}{c} \rho v \\ \rho u v \\ \rho v^2 + p \\ (e+p) v \end{array}\right] = \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] + \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] \end{equation}

where \({\bf g}\) is the gravitational force vector per unit mass, \({\bf \hat{i}},{\bf \hat{j}}\) are the unit vectors along the x and y, the viscous terms are given by

\begin{align} \tau_{ij} &= \frac{\mu}{Re_\infty} \left[ \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) - \frac{2}{3}\frac{\partial u_k}{\partial x_k} \delta_{ij} \right], \\ q_i &= - \frac{\mu}{\left(\gamma-1\right)Re_\infty Pr} \frac{\partial T}{\partial x_i} \end{align}

with \(\mu\) being the viscosity coefficient (computed using Sutherland's law), and the equation of state is

\begin{equation} e = \frac {p} {\gamma-1} + \frac{1}{2} \rho \left(u^2 + v^2\right) \end{equation}

References for the governing equations (as well as non-dimensional form):-

• 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).

Reference for the well-balanced treatment of gravitational source term:

• 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

Reference for the partitioning of the flux into its stiff (acoustic) and non-stiff (convective) components:

• 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.

Definition in file navierstokes2d.h.

---

Data Structure Documentation

NavierStokes2D

struct NavierStokes2D

Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.

Definition at line 374 of file navierstokes2d.h.

Data Fields
double	gamma	Ratio of heat capacities
char	upw_choice[_MAX_STRING_SIZE_]	choice of upwinding
double	grav_x	acceleration due to gravity in x
double	grav_y	acceleration due to gravity in y
double	rho0	reference density at zero altitude for flows with gravity
double	p0	reference pressure at zero altitude for flows with gravity
double	Re	Reynolds number
double	Pr	Prandtl number
double	Minf	Freestream Mach number
double	C1	Sutherlands law constants
double	C2	Sutherlands law constants
double	R	universal Gas constant
double *	grav_field_f	density variation function ( \(\varrho\)) for hydrostatic equilibrium for flows with gravity
double *	grav_field_g	pressure variation function ( \(\varrho\)) for hydrostatic equilibrium for flows with gravity
double *	fast_jac	"Fast" Jacobian of the flux function (comprising the acoustic modes)
double *	solution	array to store the solution at the beginning of each time step
int	HB	< Choice of hydrostatic balance for flows with gravity (1 - isothermal equilibrium, 2 - constant potential temperature 3 - stratified atmosphere with a Brunt-Vaisala frequency)
double	N_bv	the Brunt-Vaisala frequency for NavierStokes2D::HB = 3
double *	gpu_grav_field_f
double *	gpu_grav_field_g
double *	gpu_fast_jac
double *	gpu_solution

Macro Definition Documentation

_NAVIER_STOKES_2D_

#define _NAVIER_STOKES_2D_   "navierstokes2d"

2D Navier Stokes equations

Definition at line 47 of file navierstokes2d.h.

_MODEL_NDIMS_

#define _MODEL_NDIMS_   2

Number of spatial dimensions

Definition at line 53 of file navierstokes2d.h.

_MODEL_NVARS_

#define _MODEL_NVARS_   4

Number of variables per grid point

Definition at line 55 of file navierstokes2d.h.

_ROE_

#define _ROE_   "roe"

Roe's upwinding scheme

Definition at line 59 of file navierstokes2d.h.

_RF_

#define _RF_   "rf-char"

Characteristic-based Roe-fixed scheme

Definition at line 61 of file navierstokes2d.h.

_LLF_

#define _LLF_   "llf-char"

Characteristic-based local Lax-Friedrich scheme

Definition at line 63 of file navierstokes2d.h.

_SWFS_

#define _SWFS_   "steger-warming"

Steger-Warming flux splitting scheme

Definition at line 65 of file navierstokes2d.h.

_RUSANOV_

#define _RUSANOV_   "rusanov"

Rusanov's upwinding scheme

Definition at line 67 of file navierstokes2d.h.

_XDIR_

#define _XDIR_   0

dimension corresponding to the x spatial dimension

Definition at line 71 of file navierstokes2d.h.

_YDIR_

#define _YDIR_   1

dimension corresponding to the y spatial dimension

Definition at line 73 of file navierstokes2d.h.

_NavierStokes2DGetFlowVar_

#define _NavierStokes2DGetFlowVar_	(		u,
			rho,
			vx,
			vy,
			e,
			P,
			gamma
	)

Value:

{ \
    double  vsq; \
    rho = u[0]; \
    vx  = (rho==0) ? 0 : u[1] / rho; \
    vy  = (rho==0) ? 0 : u[2] / rho; \
    e   = u[3]; \
    vsq  = (vx*vx) + (vy*vy); \
    P   = (e - 0.5*rho*vsq) * (gamma-1.0); \
  }

Get the flow variables from the conserved solution vector.

\begin{equation} {\bf u} = \left[\begin{array}{c} \rho \\ \rho u \\ \rho v \\ e \end{array}\right] \end{equation}

Definition at line 81 of file navierstokes2d.h.

_NavierStokes2DSetFlux_

#define _NavierStokes2DSetFlux_	(		f,
			rho,
			vx,
			vy,
			e,
			P,
			dir
	)

Value:

{ \
    if (dir == _XDIR_) { \
      f[0] = rho * vx; \
      f[1] = rho * vx * vx + P; \
      f[2] = rho * vx * vy; \
      f[3] = (e + P) * vx; \
    } else if (dir == _YDIR_) { \
      f[0] = rho * vy; \
      f[1] = rho * vy * vx; \
      f[2] = rho * vy * vy + P; \
      f[3] = (e + P) * vy; \
    } \
  }

Compute the flux vector, given the flow variables

\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}

Definition at line 99 of file navierstokes2d.h.

_NavierStokes2DSetStiffFlux_

#define _NavierStokes2DSetStiffFlux_	(		f,
			rho,
			vx,
			vy,
			e,
			P,
			dir,
			gamma
	)

Value:

{ \
    double gamma_inv = 1.0/gamma; \
    if (dir == _XDIR_) { \
      f[0] = gamma_inv * rho * vx; \
      f[1] = gamma_inv * rho * vx * vx + P; \
      f[2] = gamma_inv * rho * vx * vy; \
      f[3] = (e + P) * vx - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vx; \
    } else if (dir == _YDIR_) { \
      f[0] = gamma_inv * rho * vy; \
      f[1] = gamma_inv * rho * vy * vx; \
      f[2] = gamma_inv * rho * vy * vy + P; \
      f[3] = (e + P) * vy - 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vy; \
    } \
  }

Compute the stiff flux vector (comprising the acoustic modes only), given the flow variables

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

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.

Definition at line 125 of file navierstokes2d.h.

_NavierStokes2DSetNonStiffFlux_

#define _NavierStokes2DSetNonStiffFlux_	(		f,
			rho,
			vx,
			vy,
			e,
			P,
			dir,
			gamma
	)

Value:

{ \
    double gamma_inv = 1.0/gamma; \
    if (dir == _XDIR_) { \
      f[0] = (gamma-1.0) * gamma_inv * rho * vx; \
      f[1] = (gamma-1.0) * gamma_inv * rho * vx * vx; \
      f[2] = (gamma-1.0) * gamma_inv * rho * vx * vy; \
      f[3] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vx; \
    } else if (dir == _YDIR_) { \
      f[0] = (gamma-1.0) * gamma_inv * rho * vy; \
      f[1] = (gamma-1.0) * gamma_inv * rho * vy * vx; \
      f[2] = (gamma-1.0) * gamma_inv * rho * vy * vy; \
      f[3] = 0.5 * gamma_inv * (gamma-1.0) * rho * (vx*vx+vy*vy) * vy; \
    } \
  }

Compute the non-stiff flux vector (comprising the entropy modes only), given the flow variables

\begin{eqnarray} dir = x, & {\bf f}\left({\bf u}\right) = \left[\begin{array}{c} \frac{\gamma-1}{\gamma}\rho u \\ \frac{\gamma-1}{\gamma}\rho u^2 \\ \frac{\gamma-11}{\gamma}\rho u v \\ \frac{1}{2} \frac{\gamma-1}{\gamma}\rho\left(u^2+v^2\right)u \end{array}\right], \\ dir = y, & {\bf f}\left({\bf u}\right) = \left[\begin{array}{c} \frac{\gamma-1}{\gamma}\rho v \\ \frac{\gamma-1}{\gamma}\rho u v \\ \frac{\gamma-1}{\gamma}\rho v^2 \\ \frac{1}{2} \frac{\gamma-1}{\gamma}\rho\left(u^2+v^2\right)v \end{array}\right] \end{eqnarray}

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.

Definition at line 152 of file navierstokes2d.h.

_NavierStokes2DRoeAverage_

#define _NavierStokes2DRoeAverage_	(		uavg,
			uL,
			uR,
			gamma
	)

Compute the Roe-average of two solutions.

Definition at line 171 of file navierstokes2d.h.

_NavierStokes2DEigenvalues_

#define _NavierStokes2DEigenvalues_	(		u,
			D,
			gamma,
			dir
	)

Value:

{ \
    double  rho,vx,vy,e,P,c,vn,vsq; \
    rho = u[0]; \
    vx  = u[1] / rho; \
    vy  = u[2] / rho; \
    e   = u[3]; \
    vsq  = (vx*vx) + (vy*vy); \
    P   = (e - 0.5*rho*vsq) * (gamma-1.0); \
    c    = sqrt(gamma*P/rho); \
    if      (dir == _XDIR_) vn = vx; \
    else if (dir == _YDIR_) vn = vy; \
    else               vn = 0.0; \
    if (dir == _XDIR_) {\
      D[0*_MODEL_NVARS_+0] = vn-c;   D[0*_MODEL_NVARS_+1] = 0;    D[0*_MODEL_NVARS_+2] = 0;      D[0*_MODEL_NVARS_+3] = 0; \
      D[1*_MODEL_NVARS_+0] = 0;      D[1*_MODEL_NVARS_+1] = vn+c; D[1*_MODEL_NVARS_+2] = 0;      D[1*_MODEL_NVARS_+3] = 0; \
      D[2*_MODEL_NVARS_+0] = 0;      D[2*_MODEL_NVARS_+1] = 0;    D[2*_MODEL_NVARS_+2] = vn;     D[2*_MODEL_NVARS_+3] = 0; \
      D[3*_MODEL_NVARS_+0] = 0;      D[3*_MODEL_NVARS_+1] = 0;    D[3*_MODEL_NVARS_+2] = 0;      D[3*_MODEL_NVARS_+3] = vn; \
    } else if (dir == _YDIR_) { \
      D[0*_MODEL_NVARS_+0] = vn-c;   D[0*_MODEL_NVARS_+1] = 0;    D[0*_MODEL_NVARS_+2] = 0;      D[0*_MODEL_NVARS_+3] = 0; \
      D[1*_MODEL_NVARS_+0] = 0;      D[1*_MODEL_NVARS_+1] = vn;   D[1*_MODEL_NVARS_+2] = 0;      D[1*_MODEL_NVARS_+3] = 0; \
      D[2*_MODEL_NVARS_+0] = 0;      D[2*_MODEL_NVARS_+1] = 0;    D[2*_MODEL_NVARS_+2] = vn+c;   D[2*_MODEL_NVARS_+3] = 0; \
      D[3*_MODEL_NVARS_+0] = 0;      D[3*_MODEL_NVARS_+1] = 0;    D[3*_MODEL_NVARS_+2] = 0;      D[3*_MODEL_NVARS_+3] = vn; \
    }\
  }

Compute the eigenvalues, given a solution vector in terms of the conserved variables. The eigenvalues are returned as a matrix D whose diagonal values are the eigenvalues. Admittedly, this is inefficient. The matrix D is stored in a row-major format.

Definition at line 213 of file navierstokes2d.h.

_NavierStokes2DLeftEigenvectors_

#define _NavierStokes2DLeftEigenvectors_	(		u,
			L,
			ga,
			dir
	)

Compute the left eigenvectors, given a solution vector in terms of the conserved variables. The eigenvectors are returned as a matrix L whose rows correspond to each eigenvector. The matrix L is stored in the row-major format.

Reference:

• Rohde, "Eigenvalues and eigenvectors of the Euler equations in general geometries", AIAA Paper 2001-2609, http://dx.doi.org/10.2514/6.2001-2609

Definition at line 247 of file navierstokes2d.h.

_NavierStokes2DRightEigenvectors_

#define _NavierStokes2DRightEigenvectors_	(		u,
			R,
			ga,
			dir
	)

Compute the right eigenvectors, given a solution vector in terms of the conserved variables. The eigenvectors are returned as a matrix R whose columns correspond to each eigenvector. The matrix R is stored in the row-major format.

Reference:

• Rohde, "Eigenvalues and eigenvectors of the Euler equations in general geometries", AIAA Paper 2001-2609, http://dx.doi.org/10.2514/6.2001-2609

Definition at line 309 of file navierstokes2d.h.

Function Documentation

NavierStokes2DInitialize()

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

NavierStokes2DCleanup()

int NavierStokes2DCleanup

(

void *

s

)

Function to clean up all allocations in the 2D Navier Stokes module.

Parameters

s	Object of type NavierStokes2D

Definition at line 11 of file NavierStokes2DCleanup.c.

{
  NavierStokes2D  *param  = (NavierStokes2D*) s;

  free(param->grav_field_f);
  free(param->grav_field_g);
  free(param->fast_jac);
  free(param->solution);
  return(0);
}