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

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

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

Go to the source code of this file.

Functions

double NavierStokes3DComputeCFL (void *, void *, double, double)
 
int NavierStokes3DFlux (double *, double *, int, void *, double)
 
int NavierStokes3DStiffFlux (double *, double *, int, void *, double)
 
int NavierStokes3DNonStiffFlux (double *, double *, int, void *, double)
 
int NavierStokes3DRoeAverage (double *, double *, double *, void *)
 
int NavierStokes3DParabolicFunction (double *, double *, void *, void *, double)
 
int NavierStokes3DSource (double *, double *, void *, void *, double)
 
int NavierStokes3DJacobian (double *, double *, void *, int, int, int)
 
int NavierStokes3DStiffJacobian (double *, double *, void *, int, int, int)
 
int NavierStokes3DLeftEigenvectors (double *, double *, void *, int)
 
int NavierStokes3DRightEigenvectors (double *, double *, void *, int)
 
int NavierStokes3DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwindRF (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwindLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwindRusanov (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwinddFRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwindFdFRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwindRusanovModified (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwinddFRusanovModified (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DUpwindFdFRusanovModified (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int NavierStokes3DGravityField (void *, void *)
 
int NavierStokes3DModifiedSolution (double *, double *, int, void *, void *, double)
 
int NavierStokes3DIBAdiabatic (void *, void *, double *, double)
 
int NavierStokes3DIBIsothermal (void *, void *, double *, double)
 
int NavierStokes3DPreStep (double *, void *, void *, double)
 
int NavierStokes3DIBForces (void *, void *, double)
 
int gpuNavierStokes3DInitialize (void *, void *)
 
int gpuNavierStokes3DFlux (double *, double *, int, void *, double)
 
int gpuNavierStokes3DParabolicFunction (double *, double *, void *, void *, double)
 
int gpuNavierStokes3DSource (double *, double *, void *, void *, double)
 
int gpuNavierStokes3DModifiedSolution (double *, double *, int, void *, void *, double)
 
int gpuNavierStokes3DUpwindRusanov (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int gpuNavierStokes3DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int gpuNavierStokes3DPreStep (double *, void *, void *, double)
 
int NavierStokes3DInitialize (void *s, void *m)
 

Detailed Description

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

Author
Debojyoti Ghosh

Definition in file NavierStokes3DInitialize.c.

Function Documentation

◆ NavierStokes3DComputeCFL()

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

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

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
dtTime step size for which to compute the CFL
tTime

Definition at line 16 of file NavierStokes3DComputeCFL.c.

22 {
23  HyPar *solver = (HyPar*) s;
24  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
25 
26  int *dim = solver->dim_local;
27  int ghosts = solver->ghosts;
28  int ndims = solver->ndims;
29  int index[ndims];
30  double *u = solver->u;
31 
32  double max_cfl = 0;
33  int done = 0; _ArraySetValue_(index,ndims,0);
34  while (!done) {
35  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
36  double rho, vx, vy, vz, e, P, c, dxinv, dyinv, dzinv, local_cfl[3];
38 
39  c = sqrt(param->gamma*P/rho); /* speed of sound */
40  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
41  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv); /* 1/dy */
42  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv); /* 1/dz */
43 
44  local_cfl[_XDIR_] = (absolute(vx)+c)*dt*dxinv; /* local cfl for this grid point (x) */
45  local_cfl[_YDIR_] = (absolute(vy)+c)*dt*dyinv; /* local cfl for this grid point (y) */
46  local_cfl[_ZDIR_] = (absolute(vz)+c)*dt*dzinv; /* local cfl for this grid point (z) */
47  if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
48  if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];
49  if (local_cfl[_ZDIR_] > max_cfl) max_cfl = local_cfl[_ZDIR_];
50 
51  _ArrayIncrementIndex_(ndims,dim,index,done);
52  }
53 
54  return(max_cfl);
55 }
#define absolute(a)
Definition: math_ops.h:32
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * u
Definition: hypar.h:116
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
double * dxinv
Definition: hypar.h:110
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DFlux()

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

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

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

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

Parameters
fArray to hold the computed flux vector (same layout as u)
uArray with the solution vector
dirSpatial dimension (x, y, or z) for which to compute the flux
sSolver object of type HyPar
tCurrent simulation time

Definition at line 23 of file NavierStokes3DFlux.c.

30 {
31  HyPar *solver = (HyPar*) s;
32  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
33  int i;
34 
35  int *dim = solver->dim_local;
36  int ghosts = solver->ghosts;
37  int ndims = solver->ndims;
38  int index[ndims], bounds[ndims], offset[ndims];
39 
40  /* set bounds for array index to include ghost points */
41  _ArrayCopy1D_(dim,bounds,ndims);
42  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
43 
44  /* set offset such that index is compatible with ghost point arrangement */
45  _ArraySetValue_(offset,ndims,-ghosts);
46 
47  int done = 0; _ArraySetValue_(index,ndims,0);
48  while (!done) {
49  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
50  double rho, vx, vy, vz, e, P;
53  _ArrayIncrementIndex_(ndims,bounds,index,done);
54  }
55 
56  return(0);
57 }
#define _NavierStokes3DSetFlux_(f, stride, rho, vx, vy, vz, e, P, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _ArrayCopy1D_(x, y, size)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DStiffFlux()

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

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

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

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

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

Parameters
fArray to hold the computed flux vector (same layout as u)
uArray with the solution vector
dirSpatial dimension (x,y, or z) for which to compute the flux
sSolver object of type HyPar
tCurrent simulation time

Definition at line 70 of file NavierStokes3DFlux.c.

77 {
78  HyPar *solver = (HyPar*) s;
79  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
80  int *dim = solver->dim_local;
81  int ghosts = solver->ghosts;
82  static const int JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
83  static int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
84 
85  /* set bounds for array index to include ghost points */
86  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,_MODEL_NDIMS_);
87  /* set offset such that index is compatible with ghost point arrangement */
88  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
89 
90  int done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
91  while (!done) {
92  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
93  double *Af = param->fast_jac+(_MODEL_NDIMS_*p+dir)*JacSize;
95  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
96  }
97 
98  return(0);
99 }
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _ArrayAddCopy1D_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)

◆ NavierStokes3DNonStiffFlux()

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

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

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

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

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

Parameters
fArray to hold the computed flux vector (same layout as u)
uArray with the solution vector
dirSpatial dimension (x,y, or z) for which to compute the flux
sSolver object of type HyPar
tCurrent simulation time

Definition at line 111 of file NavierStokes3DFlux.c.

118 {
119  HyPar *solver = (HyPar*) s;
120  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
121  int *dim = solver->dim_local;
122  int ghosts = solver->ghosts;
123  static const int JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
124  static int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
125  static double ftot[_MODEL_NVARS_], fstiff[_MODEL_NVARS_];
126 
127  /* set bounds for array index to include ghost points */
128  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,_MODEL_NDIMS_);
129  /* set offset such that index is compatible with ghost point arrangement */
130  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
131 
132  int done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
133  while (!done) {
134  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
135  /* compute total flux */
136  double rho, vx, vy, vz, e, P;
138  _NavierStokes3DSetFlux_(ftot,_NavierStokes3D_stride_,rho,vx,vy,vz,e,P,dir);
139  /* compute stiff stuff */
140  double *Af = param->fast_jac+(_MODEL_NDIMS_*p+dir)*JacSize;
141  MatVecMult5(_MODEL_NVARS_,fstiff,Af,(u+_MODEL_NVARS_*p));
142  /* subtract stiff flux from total flux */
143  _ArraySubtract1D_((f+_MODEL_NVARS_*p),ftot,fstiff,_MODEL_NVARS_);
144  /* Done */
145  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
146  }
147 
148  return(0);
149 }
#define _NavierStokes3DSetFlux_(f, stride, rho, vx, vy, vz, e, P, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayAddCopy1D_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DRoeAverage()

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

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

Parameters
uavgThe computed Roe-averaged state
uLLeft state (conserved variables)
uRRight state (conserved variables)
pObject of type NavierStokes3D with physics-related variables

Definition at line 18 of file NavierStokes3DFunctions.c.

24 {
25  NavierStokes3D *param = (NavierStokes3D*) p;
27  return(0);
28 }
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int _NavierStokes3D_stride_

◆ NavierStokes3DParabolicFunction()

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

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

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

where

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

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

Reference:

  • Tannehill, Anderson and Pletcher, Computational Fluid Mechanics and Heat Transfer, Chapter 5, Section 5.1.7 (However, the non-dimensional velocity and the Reynolds number is based on speed of sound, instead of the freestream velocity).
Parameters
parArray to hold the computed viscous terms
uSolution vector array
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time

Definition at line 50 of file NavierStokes3DParabolicFunction.c.

57 {
58  HyPar *solver = (HyPar*) s;
59  MPIVariables *mpi = (MPIVariables*) m;
60  NavierStokes3D *physics = (NavierStokes3D*) solver->physics;
61  int i,j,k,v;
63 
64  int ghosts = solver->ghosts;
65  int imax = solver->dim_local[0];
66  int jmax = solver->dim_local[1];
67  int kmax = solver->dim_local[2];
68  int *dim = solver->dim_local;
69  int size = solver->npoints_local_wghosts;
70 
71  _ArraySetValue_(par,size*_MODEL_NVARS_,0.0);
72  if (physics->Re <= 0) return(0); /* inviscid flow */
73  solver->count_par++;
74 
75  static double two_third = 2.0/3.0;
76  double inv_gamma_m1 = 1.0 / (physics->gamma-1.0);
77  double inv_Re = 1.0 / physics->Re;
78  double inv_Pr = 1.0 / physics->Pr;
79 
80 #if defined(CPU_STAT)
81  clock_t startEvent, stopEvent;
82 #endif
83 
84  double *Q; /* primitive variables */
85  Q = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
86 
87 #if defined(CPU_STAT)
88  startEvent = clock();
89 #endif
90 
91  for (i=-ghosts; i<(imax+ghosts); i++) {
92  for (j=-ghosts; j<(jmax+ghosts); j++) {
93  for (k=-ghosts; k<(kmax+ghosts); k++) {
94  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
95  double energy,pressure;
96  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
99  Q[p+0],
100  Q[p+1],
101  Q[p+2],
102  Q[p+3],
103  energy,
104  pressure,
105  physics->gamma);
106  Q[p+4] = physics->gamma*pressure/Q[p+0]; /* temperature */
107  }
108  }
109  }
110 
111 #if defined(CPU_STAT)
112  stopEvent = clock();
113  printf("%-50s CPU time (secs) = %.6f\n",
114  "NaverStokes3DParabolicFunction_Q",(double)(stopEvent-startEvent)/CLOCKS_PER_SEC);
115 #endif
116 
117  double *QDerivX = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
118  double *QDerivY = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
119  double *QDerivZ = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
120 
121  IERR solver->FirstDerivativePar(QDerivX,Q,_XDIR_,1,solver,mpi); CHECKERR(ierr);
122  IERR solver->FirstDerivativePar(QDerivY,Q,_YDIR_,1,solver,mpi); CHECKERR(ierr);
123  IERR solver->FirstDerivativePar(QDerivZ,Q,_ZDIR_,1,solver,mpi); CHECKERR(ierr);
124 
126  solver->ghosts,mpi,QDerivX); CHECKERR(ierr);
128  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
130  solver->ghosts,mpi,QDerivY); CHECKERR(ierr);
131 
132  for (i=-ghosts; i<(imax+ghosts); i++) {
133  for (j=-ghosts; j<(jmax+ghosts); j++) {
134  for (k=-ghosts; k<(kmax+ghosts); k++) {
135  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
136  double dxinv, dyinv, dzinv;
137  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
138  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
139  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
140  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
141  _ArrayScale1D_((QDerivX+p),dxinv,_MODEL_NVARS_);
142  _ArrayScale1D_((QDerivY+p),dyinv,_MODEL_NVARS_);
143  _ArrayScale1D_((QDerivZ+p),dzinv,_MODEL_NVARS_);
144  }
145  }
146  }
147 
148  double *FViscous = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
149  double *FDeriv = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
150 
151  /* Along X */
152  for (i=-ghosts; i<(imax+ghosts); i++) {
153  for (j=-ghosts; j<(jmax+ghosts); j++) {
154  for (k=-ghosts; k<(kmax+ghosts); k++) {
155  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
156  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
157 
158  double uvel, vvel, wvel, T, Tx,
159  ux, uy, uz, vx, vy, wx, wz;
160  uvel = (Q+p)[1];
161  vvel = (Q+p)[2];
162  wvel = (Q+p)[3];
163  T = (Q+p)[4];
164  Tx = (QDerivX+p)[4];
165  ux = (QDerivX+p)[1];
166  vx = (QDerivX+p)[2];
167  wx = (QDerivX+p)[3];
168  uy = (QDerivY+p)[1];
169  vy = (QDerivY+p)[2];
170  uz = (QDerivZ+p)[1];
171  wz = (QDerivZ+p)[3];
172 
173  /* calculate viscosity coeff based on Sutherland's law */
174  double mu = raiseto(T, 0.76);
175  /*
176  if (p/_MODEL_NVARS_ == 49841) {
177  printf("[CPU] 49841 mu = %.20e T = %.20e log(0.76) = %.20e T*log(0.76) = %.20e\n",
178  mu, T, log(0.76), T*log(0.76));
179  }
180  */
181 
182  double tau_xx, tau_xy, tau_xz, qx;
183  tau_xx = two_third * (mu*inv_Re) * (2*ux - vy - wz);
184  tau_xy = (mu*inv_Re) * (uy + vx);
185  tau_xz = (mu*inv_Re) * (uz + wx);
186  qx = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tx;
187 
188  (FViscous+p)[0] = 0.0;
189  (FViscous+p)[1] = tau_xx;
190  (FViscous+p)[2] = tau_xy;
191  (FViscous+p)[3] = tau_xz;
192  (FViscous+p)[4] = uvel*tau_xx + vvel*tau_xy + wvel*tau_xz + qx;
193  }
194  }
195  }
196 
197  IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
198  for (i=0; i<imax; i++) {
199  for (j=0; j<jmax; j++) {
200  for (k=0; k<kmax; k++) {
201  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
202  double dxinv;
203  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
204  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
205  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dxinv * (FDeriv+p)[v] );
206  }
207  }
208  }
209 
210  /* Along Y */
211  for (i=-ghosts; i<(imax+ghosts); i++) {
212  for (j=-ghosts; j<(jmax+ghosts); j++) {
213  for (k=-ghosts; k<(kmax+ghosts); k++) {
214  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
215  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
216 
217  double uvel, vvel, wvel, T, Ty,
218  ux, uy, vx, vy, vz, wy, wz;
219  uvel = (Q+p)[1];
220  vvel = (Q+p)[2];
221  wvel = (Q+p)[3];
222  T = (Q+p)[4];
223  Ty = (QDerivY+p)[4];
224  ux = (QDerivX+p)[1];
225  vx = (QDerivX+p)[2];
226  uy = (QDerivY+p)[1];
227  vy = (QDerivY+p)[2];
228  wy = (QDerivY+p)[3];
229  vz = (QDerivZ+p)[2];
230  wz = (QDerivZ+p)[3];
231 
232  /* calculate viscosity coeff based on Sutherland's law */
233  double mu = raiseto(T, 0.76);
234 
235  double tau_yx, tau_yy, tau_yz, qy;
236  tau_yx = (mu*inv_Re) * (uy + vx);
237  tau_yy = two_third * (mu*inv_Re) * (-ux + 2*vy - wz);
238  tau_yz = (mu*inv_Re) * (vz + wy);
239  qy = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Ty;
240 
241  (FViscous+p)[0] = 0.0;
242  (FViscous+p)[1] = tau_yx;
243  (FViscous+p)[2] = tau_yy;
244  (FViscous+p)[3] = tau_yz;
245  (FViscous+p)[4] = uvel*tau_yx + vvel*tau_yy + wvel*tau_yz + qy;
246  }
247  }
248  }
249 
250  IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
251  for (i=0; i<imax; i++) {
252  for (j=0; j<jmax; j++) {
253  for (k=0; k<kmax; k++) {
254  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
255  double dyinv;
256  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
257  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
258  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
259  }
260  }
261  }
262 
263  /* Along Z */
264  for (i=-ghosts; i<(imax+ghosts); i++) {
265  for (j=-ghosts; j<(jmax+ghosts); j++) {
266  for (k=-ghosts; k<(kmax+ghosts); k++) {
267  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
268  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
269 
270  double uvel, vvel, wvel, T, Tz,
271  ux, uz, vy, vz, wx, wy, wz;
272  uvel = (Q+p)[1];
273  vvel = (Q+p)[2];
274  wvel = (Q+p)[3];
275  T = (Q+p)[4];
276  Tz = (QDerivZ+p)[4];
277  ux = (QDerivX+p)[1];
278  wx = (QDerivX+p)[3];
279  vy = (QDerivY+p)[2];
280  wy = (QDerivY+p)[3];
281  uz = (QDerivZ+p)[1];
282  vz = (QDerivZ+p)[2];
283  wz = (QDerivZ+p)[3];
284 
285  /* calculate viscosity coeff based on Sutherland's law */
286  double mu = raiseto(T,0.76);
287 
288  double tau_zx, tau_zy, tau_zz, qz;
289  tau_zx = (mu*inv_Re) * (uz + wx);
290  tau_zy = (mu*inv_Re) * (vz + wy);
291  tau_zz = two_third * (mu*inv_Re) * (-ux - vy + 2*wz);
292  qz = ( mu*inv_Re * inv_gamma_m1 * inv_Pr ) * Tz;
293 
294  (FViscous+p)[0] = 0.0;
295  (FViscous+p)[1] = tau_zx;
296  (FViscous+p)[2] = tau_zy;
297  (FViscous+p)[3] = tau_zz;
298  (FViscous+p)[4] = uvel*tau_zx + vvel*tau_zy + wvel*tau_zz + qz;
299  }
300  }
301  }
302 
303  IERR solver->FirstDerivativePar(FDeriv,FViscous,_ZDIR_,-1,solver,mpi); CHECKERR(ierr);
304  for (i=0; i<imax; i++) {
305  for (j=0; j<jmax; j++) {
306  for (k=0; k<kmax; k++) {
307  int p,index[3]; index[0]=i; index[1]=j; index[2]=k;
308  double dzinv;
309  _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
310  _GetCoordinate_(_ZDIR_,index[_ZDIR_],dim,ghosts,solver->dxinv,dzinv);
311  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dzinv * (FDeriv+p)[v] );
312  }
313  }
314  }
315 
316  free(Q);
317  free(QDerivX);
318  free(QDerivY);
319  free(QDerivZ);
320  free(FViscous);
321  free(FDeriv);
322 
323  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,_MODEL_NVARS_);
324  return(0);
325 }
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
int flag_ib
Definition: hypar.h:441
double * iblank
Definition: hypar.h:436
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int count_par
Definition: hypar.h:418
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
int(* FirstDerivativePar)(double *, double *, int, int, void *, void *)
Definition: hypar.h:243
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define _ArrayScale1D_(x, a, size)
#define raiseto(x, a)
Definition: math_ops.h:37
#define _ArrayBlockMultiply_(x, a, n, bs)
#define _MODEL_NVARS_
Definition: euler1d.h:58
int npoints_local_wghosts
Definition: hypar.h:42
#define _DECLARE_IERR_
Definition: basic.h:17
double * dxinv
Definition: hypar.h:110
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DSource()

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

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

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

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

References:

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, AIAA Journal, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580.
Parameters
sourceArray to hold the computed source
uSolution vector array
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time

Definition at line 38 of file NavierStokes3DSource.c.

45 {
46  HyPar *solver = (HyPar* ) s;
47  MPIVariables *mpi = (MPIVariables*) m;
48  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
49 
50  if ((param->grav_x == 0.0) && (param->grav_y == 0.0) && (param->grav_z == 0.0))
51  return(0); /* no gravitational forces */
52 
53  int v, done, p, p1, p2, dir;
54  double *SourceI = solver->fluxI; /* interace source term */
55  double *SourceC = solver->fluxC; /* cell-centered source term */
56  double *SourceL = solver->fL;
57  double *SourceR = solver->fR;
58 
59  int ghosts = solver->ghosts;
60  int *dim = solver->dim_local;
61  double *x = solver->x;
62  double *dxinv = solver->dxinv;
63  double RT = param->p0 / param->rho0;
64  static int index[_MODEL_NDIMS_],index1[_MODEL_NDIMS_],index2[_MODEL_NDIMS_],dim_interface[_MODEL_NDIMS_];
65  static double grav[_MODEL_NDIMS_];
66 
67  grav[_XDIR_] = param->grav_x;
68  grav[_YDIR_] = param->grav_y;
69  grav[_ZDIR_] = param->grav_z;
70  for (dir = 0; dir < _MODEL_NDIMS_; dir++) {
71  if (grav[dir] != 0.0) {
72  /* set interface dimensions */
73  _ArrayCopy1D_(dim,dim_interface,_MODEL_NDIMS_); dim_interface[dir]++;
74  /* calculate the split source function exp(-phi/RT) */
75  IERR NavierStokes3DSourceFunction(SourceC,u,x,solver,mpi,t,dir); CHECKERR(ierr);
76  /* calculate the left and right interface source terms */
77  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,dir,solver,mpi,0); CHECKERR(ierr);
78  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,dir,solver,mpi,0); CHECKERR(ierr);
79  /* calculate the final interface source term */
80  IERR NavierStokes3DSourceUpwind(SourceI,SourceL,SourceR,u,dir,solver,t);
81  /* calculate the final cell-centered source term */
82  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
83  while (!done) {
84  _ArrayCopy1D_(index,index1,_MODEL_NDIMS_);
85  _ArrayCopy1D_(index,index2,_MODEL_NDIMS_); index2[dir]++;
86  _ArrayIndex1D_(_MODEL_NDIMS_,dim ,index ,ghosts,p );
87  _ArrayIndex1D_(_MODEL_NDIMS_,dim_interface,index1,0 ,p1);
88  _ArrayIndex1D_(_MODEL_NDIMS_,dim_interface,index2,0 ,p2);
89 
90  double dx_inverse; _GetCoordinate_(dir,index[dir],dim,ghosts,dxinv,dx_inverse);
91  double rho, vel[_MODEL_NDIMS_], e, P;
92  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],e,P,param->gamma);
93  double term[_MODEL_NVARS_] = {0.0, rho*RT*(dir==_XDIR_), rho*RT*(dir==_YDIR_), rho*RT*(dir==_ZDIR_), rho*RT*vel[dir]};
94  for (v=0; v<_MODEL_NVARS_; v++) {
95  source[_MODEL_NVARS_*p+v] += ( (term[v]*param->grav_field_f[p])
96  * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse );
97  }
98  vel[0] = P; /* useless statement to avoid compiler warnings */
99  _ArrayIncrementIndex_(_MODEL_NDIMS_,dim,index,done);
100  }
101  }
102  }
103 
104  return(0);
105 }
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * x
Definition: hypar.h:107
double * fluxI
Definition: hypar.h:136
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int NavierStokes3DSourceUpwind(double *, double *, double *, double *, int, void *, double)
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int(* InterpolateInterfacesHyp)(double *, double *, double *, double *, int, int, void *, void *, int)
Definition: hypar.h:224
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
double * fR
Definition: hypar.h:139
double * fL
Definition: hypar.h:139
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
int NavierStokes3DSourceFunction(double *, double *, double *, void *, void *, double, int)
Structure of MPI-related variables.
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _ArrayCopy1D_(x, y, size)
double * fluxC
Definition: hypar.h:128
double * grav_field_f
double * dxinv
Definition: hypar.h:110
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DJacobian()

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

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

Parameters
JacJacobian matrix: 1D array of size nvar^2 = 25
usolution at a grid point (array of size nvar = 5)
pobject containing the physics-related parameters
dirdimension (0 -> x, 1 -> y, 2 -> z)
nvarsnumber of vector components
upw0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 15 of file NavierStokes3DJacobian.c.

25 {
26  NavierStokes3D *param = (NavierStokes3D*) p;
29 
30  /* get the eigenvalues and left,right eigenvectors */
34 
35  int aupw = absolute(upw), k;
36  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]) );
37  k = 6; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
38  k = 12; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
39  k = 18; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
40  k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
41 
42  MatMult5(_MODEL_NVARS_,DL,D,L);
43  MatMult5(_MODEL_NVARS_,Jac,R,DL);
44 
45  return(0);
46 }
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define min(a, b)
Definition: math_ops.h:14
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
static const int _NavierStokes3D_stride_
#define MatMult5(N, A, X, Y)
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes3DStiffJacobian()

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

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

Parameters
JacJacobian matrix: 1D array of size nvar^2 = 25
usolution at a grid point (array of size nvar = 5)
pobject containing the physics-related parameters
dirdimension (0 -> x, 1 -> y, 2 -> z)
nvarsnumber of vector components
upw0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 53 of file NavierStokes3DJacobian.c.

63 {
64  NavierStokes3D *param = (NavierStokes3D*) p;
67 
68  /* get the eigenvalues and left,right eigenvectors */
72 
73  int aupw = absolute(upw), k;
74  if (dir == _XDIR_) {
75  k = 0; D[k] = 0.0;
76  k = 6; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
77  k = 12; D[k] = 0.0;
78  k = 18; D[k] = 0.0;
79  k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
80  } else if (dir == _YDIR_) {
81  k = 0; D[k] = 0.0;
82  k = 6; D[k] = 0.0;
83  k = 12; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
84  k = 18; D[k] = 0.0;
85  k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
86  } else if (dir == _ZDIR_) {
87  k = 0; D[k] = 0.0;
88  k = 6; D[k] = 0.0;
89  k = 12; D[k] = 0.0;
90  k = 18; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
91  k = 24; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
92  }
93 
94  MatMult5(_MODEL_NVARS_,DL,D,L);
95  MatMult5(_MODEL_NVARS_,Jac,R,DL);
96 
97  return(0);
98 }
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define min(a, b)
Definition: math_ops.h:14
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _YDIR_
Definition: euler2d.h:41
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
static const int _NavierStokes3D_stride_
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define max(a, b)
Definition: math_ops.h:18

◆ NavierStokes3DLeftEigenvectors()

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

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

Parameters
uConserved solution at a grid point
LArray of size nvar^2 = 5^2 to save the matrix of left eigenvectors in (row-major format).
pObject of type NavierStokes3D with physics-related variables
dirSpatial dimension (x, y, or z)

Definition at line 20 of file NavierStokes3DEigen.c.

27 {
28  NavierStokes3D *param = (NavierStokes3D*) p;
30  return(0);
31 }
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
static const int _NavierStokes3D_stride_

◆ NavierStokes3DRightEigenvectors()

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

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

Parameters
uConserved solution at a grid point
RArray of size nvar^2 = 5^2 to save the matrix of right eigenvectors in (row-major format).
pObject of type NavierStokes3D with physics-related variables
dirSpatial dimension (x, y, or z)

Definition at line 39 of file NavierStokes3DEigen.c.

46 {
47  NavierStokes3D *param = (NavierStokes3D*) p;
49  return(0);
50 }
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
static const int _NavierStokes3D_stride_

◆ NavierStokes3DUpwindRoe()

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

Roe's upwinding scheme.

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

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

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 40 of file NavierStokes3DUpwind.c.

51 {
52  HyPar *solver = (HyPar*) s;
53  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
54  int done;
55 
56  int *dim = solver->dim_local;
57 
58  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
59  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
60  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
64 
65  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
66 
67 #if defined(CPU_STAT)
68  clock_t startEvent, stopEvent;
69  startEvent = clock();
70 #endif
71 
72  while (!done) {
73  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
74  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
75  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
76  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
77  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
78  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
79  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
80  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
81 
82  /* Roe's upwinding scheme */
83 
84  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
85  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
86  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
87  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
88  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
89 
91  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
92  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
93  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
94 
95  /* Harten's Entropy Fix - Page 362 of Leveque */
96  int k;
97  double delta = 0.000001, delta2 = delta*delta;
98  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
99  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
100  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
101  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
102  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
103 
104  MatMult5(_MODEL_NVARS_,DL,D,L);
105  MatMult5(_MODEL_NVARS_,modA,R,DL);
106  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
107 
108  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
109  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
110  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
111  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
112  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
113  }
114  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
115  }
116 
117 #if defined(CPU_STAT)
118  stopEvent = clock();
119  printf("%-50s CPU time (secs) = %.6f dir = %d\n",
120  "NavierStokes3DUpwindRoe", (double)(stopEvent-startEvent)/CLOCKS_PER_SEC, dir);
121 #endif
122 
123  return(0);
124 }
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D3_(x, y, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _ArrayCopy1D_(x, y, size)

◆ NavierStokes3DUpwindRF()

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

Characteristic-based Roe-fixed upwinding scheme.

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

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

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

Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 140 of file NavierStokes3DUpwind.c.

151 {
152  HyPar *solver = (HyPar*) s;
153  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
154  int done,k;
155 
156  int *dim = solver->dim_local;
157 
158  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
159  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
160  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
162 
163  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
164  while (!done) {
165  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
166  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
167  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
168  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
169  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
170 
171  /* Roe-Fixed upwinding scheme */
172 
174 
175  _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
176  _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
177 
178  /* calculate characteristic fluxes and variables */
183 
184  double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
186  eigL[0] = D[0];
187  eigL[1] = D[6];
188  eigL[2] = D[12];
189  eigL[3] = D[18];
190  eigL[4] = D[24];
192  eigR[0] = D[0];
193  eigR[1] = D[6];
194  eigR[2] = D[12];
195  eigR[3] = D[18];
196  eigR[4] = D[24];
197  _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
198  eigC[0] = D[0];
199  eigC[1] = D[6];
200  eigC[2] = D[12];
201  eigC[3] = D[18];
202  eigC[4] = D[24];
203 
204  for (k = 0; k < _MODEL_NVARS_; k++) {
205  if ((eigL[k] > 0) && (eigC[k] > 0) && (eigR[k] > 0)) fc[k] = fcL[k];
206  else if ((eigL[k] < 0) && (eigC[k] < 0) && (eigR[k] < 0)) fc[k] = fcR[k];
207  else {
208  double alpha = max3(absolute(eigL[k]),absolute(eigC[k]),absolute(eigR[k]));
209  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
210  }
211  }
212 
213  /* calculate the interface flux from the characteristic flux */
215  }
216  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
217  }
218 
219  return(0);
220 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D3_(x, y, size)
static const int _NavierStokes3D_stride_
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)

◆ NavierStokes3DUpwindLLF()

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

Characteristic-based local Lax-Friedrich upwinding scheme.

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

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

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

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 244 of file NavierStokes3DUpwind.c.

255 {
256  HyPar *solver = (HyPar*) s;
257  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
258  int done;
259 
260  int *dim = solver->dim_local;
261 
262  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
263  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
264  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
266 
267  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
268  while (!done) {
269  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
270  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
271  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
272  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
273  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
274 
275  /* Local Lax-Friedrich upwinding scheme */
276 
278 
279  _NavierStokes3DLeftEigenvectors_(uavg,dummy,L,param->gamma,dir);
280  _NavierStokes3DRightEigenvectors_(uavg,dummy,R,param->gamma,dir);
281 
282  /* calculate characteristic fluxes and variables */
287 
288  double eigL[_MODEL_NVARS_],eigC[_MODEL_NVARS_],eigR[_MODEL_NVARS_];
290  eigL[0] = D[0];
291  eigL[1] = D[6];
292  eigL[2] = D[12];
293  eigL[3] = D[18];
294  eigL[4] = D[24];
296  eigR[0] = D[0];
297  eigR[1] = D[6];
298  eigR[2] = D[12];
299  eigR[3] = D[18];
300  eigR[4] = D[24];
301  _NavierStokes3DEigenvalues_(uavg,dummy,D,param->gamma,dir);
302  eigC[0] = D[0];
303  eigC[1] = D[6];
304  eigC[2] = D[12];
305  eigC[3] = D[18];
306  eigC[4] = D[24];
307 
308  double alpha;
309  alpha = max3(absolute(eigL[0]),absolute(eigC[0]),absolute(eigR[0]));
310  fc[0] = 0.5 * (fcL[0] + fcR[0] + alpha * (ucL[0]-ucR[0]));
311  alpha = max3(absolute(eigL[1]),absolute(eigC[1]),absolute(eigR[1]));
312  fc[1] = 0.5 * (fcL[1] + fcR[1] + alpha * (ucL[1]-ucR[1]));
313  alpha = max3(absolute(eigL[2]),absolute(eigC[2]),absolute(eigR[2]));
314  fc[2] = 0.5 * (fcL[2] + fcR[2] + alpha * (ucL[2]-ucR[2]));
315  alpha = max3(absolute(eigL[3]),absolute(eigC[3]),absolute(eigR[3]));
316  fc[3] = 0.5 * (fcL[3] + fcR[3] + alpha * (ucL[3]-ucR[3]));
317  alpha = max3(absolute(eigL[4]),absolute(eigC[4]),absolute(eigR[4]));
318  fc[4] = 0.5 * (fcL[4] + fcR[4] + alpha * (ucL[4]-ucR[4]));
319 
320  /* calculate the interface flux from the characteristic flux */
322  }
323  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
324  }
325 
326  return(0);
327 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D3_(x, y, size)
static const int _NavierStokes3D_stride_
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)

◆ NavierStokes3DUpwindRusanov()

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

Rusanov's upwinding scheme.

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

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

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

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

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 349 of file NavierStokes3DUpwind.c.

360 {
361  HyPar *solver = (HyPar*) s;
362  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
363  int *dim = solver->dim_local, done;
364 
365  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
366  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
367  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
368 
369  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
370 
371  static int index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
372  indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_];
373 
374  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
375  while (!done) {
376  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
377  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
378 
379  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
380  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
381  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
382  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
383  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
384  int q = p*_MODEL_NVARS_;
385 
386  /* Modified Rusanov's upwinding scheme */
387 
388  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
389  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
390  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
391  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
392  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
393 
395 
396  double c, vel[_MODEL_NDIMS_], rho,E,P;
397  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
398  c = sqrt(param->gamma*P/rho);
399  double alphaL = c + absolute(vel[dir]);
400  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
401  c = sqrt(param->gamma*P/rho);
402  double alphaR = c + absolute(vel[dir]);
403  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
404  c = sqrt(param->gamma*P/rho);
405  double alphaavg = c + absolute(vel[dir]);
406 
407  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
408  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
409 
410  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-alpha*udiff[0];
411  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-alpha*udiff[1];
412  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-alpha*udiff[2];
413  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-alpha*udiff[3];
414  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-alpha*udiff[4];
415  }
416  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
417  }
418 
419  return(0);
420 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _ArrayCopy1D_(x, y, size)
#define max(a, b)
Definition: math_ops.h:18
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DUpwinddFRoe()

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

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

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 430 of file NavierStokes3DUpwind.c.

441 {
442  HyPar *solver = (HyPar*) s;
443  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
444  int done;
445 
446  int *dim = solver->dim_local;
447  double *uref = param->solution;
448 
449  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
450  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
451  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
455 
456  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
457  while (!done) {
458  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
459  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
460  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
461  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
462  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
463  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
464  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
465  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
466 
467  /* Roe's upwinding scheme */
468 
469  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
470  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
471  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
472  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
473  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
474 
476  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
477  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
478  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
479 
480  /* Harten's Entropy Fix - Page 362 of Leveque */
481  int k;
482  double delta = 0.000001, delta2 = delta*delta;
483  if (dir == _XDIR_) {
484  k=0; D[k] = 0.0;
485  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
486  k=12; D[k] = 0.0;
487  k=18; D[k] = 0.0;
488  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
489  } else if (dir == _YDIR_) {
490  k=0; D[k] = 0.0;
491  k=6; D[k] = 0.0;
492  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
493  k=18; D[k] = 0.0;
494  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
495  } else if (dir == _ZDIR_) {
496  k=0; D[k] = 0.0;
497  k=6; D[k] = 0.0;
498  k=12; D[k] = 0.0;
499  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
500  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
501  }
502 
503  MatMult5(_MODEL_NVARS_,DL,D,L);
504  MatMult5(_MODEL_NVARS_,modA,R,DL);
505  MatVecMult5(_MODEL_NVARS_,udiss,modA,udiff);
506 
507  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
508  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
509  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
510  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
511  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
512  }
513  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
514  }
515 
516  return(0);
517 }
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D3_(x, y, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _ArrayCopy1D_(x, y, size)

◆ NavierStokes3DUpwindFdFRoe()

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

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

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x, y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 527 of file NavierStokes3DUpwind.c.

538 {
539  HyPar *solver = (HyPar*) s;
540  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
541  int done;
542 
543  int *dim = solver->dim_local;
544  double *uref = param->solution;
545 
546  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
547  _ArrayCopy1D3_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
548  _ArrayCopy1D3_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
552 
553  done = 0; int index_outer[3] = {0,0,0}, index_inter[3];
554  while (!done) {
555  _ArrayCopy1D3_(index_outer,index_inter,_MODEL_NDIMS_);
556  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
557  int p; _ArrayIndex1D3_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
558  int indexL[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
559  int indexR[_MODEL_NDIMS_]; _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
560  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
561  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
562  int k;
563  double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_],
564  udiss_total[_MODEL_NVARS_],udiss_stiff[_MODEL_NVARS_];
565  double delta = 0.000001, delta2 = delta*delta;
566 
567  /* Roe's upwinding scheme */
568 
569  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
570  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
571  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
572  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
573  udiff[4] = 0.5 * (uR[_MODEL_NVARS_*p+4] - uL[_MODEL_NVARS_*p+4]);
574 
575  /* Compute total dissipation */
577  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
578  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
579  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
580  k=0; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
581  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
582  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
583  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
584  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
585  MatMult5(_MODEL_NVARS_,DL,D,L);
586  MatMult5(_MODEL_NVARS_,modA,R,DL);
587  MatVecMult5(_MODEL_NVARS_,udiss_total,modA,udiff);
588 
589  /* Compute dissipation corresponding to acoustic modes */
591  _NavierStokes3DEigenvalues_ (uavg,dummy,D,param->gamma,dir);
592  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
593  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
594  if (dir == _XDIR_) {
595  k=0; D[k] = 0.0;
596  k=6; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
597  k=12; D[k] = 0.0;
598  k=18; D[k] = 0.0;
599  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
600  } else if (dir == _YDIR_) {
601  k=0; D[k] = 0.0;
602  k=6; D[k] = 0.0;
603  k=12; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
604  k=18; D[k] = 0.0;
605  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
606  } else if (dir == _ZDIR_) {
607  k=0; D[k] = 0.0;
608  k=6; D[k] = 0.0;
609  k=12; D[k] = 0.0;
610  k=18; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
611  k=24; D[k] = (absolute(D[k]) < delta ? (D[k]*D[k]+delta2)/(2*delta) : absolute(D[k]) );
612  }
613  MatMult5(_MODEL_NVARS_,DL,D,L);
614  MatMult5(_MODEL_NVARS_,modA,R,DL);
615  MatVecMult5(_MODEL_NVARS_,udiss_stiff,modA,udiff);
616 
617  /* Compute the dissipation term for the entropy modes */
618  _ArraySubtract1D_(udiss,udiss_total,udiss_stiff,_MODEL_NVARS_);
619 
620  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
621  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
622  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
623  fI[_MODEL_NVARS_*p+3] = 0.5 * (fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3]) - udiss[3];
624  fI[_MODEL_NVARS_*p+4] = 0.5 * (fL[_MODEL_NVARS_*p+4]+fR[_MODEL_NVARS_*p+4]) - udiss[4];
625  }
626  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
627  }
628 
629  return(0);
630 }
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
#define _ArraySubtract1D_(x, a, b, size)
#define _ArrayIndex1D3_(N, imax, i, ghost, index)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayCopy1D3_(x, y, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _ArrayCopy1D_(x, y, size)

◆ NavierStokes3DUpwindRusanovModified()

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

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

  • Ghosh, D., Constantinescu, E. M., "Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning", Submitted (http://arxiv.org/abs/1510.05751).
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 638 of file NavierStokes3DUpwind.c.

649 {
650  HyPar *solver = (HyPar*) s;
651  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
652  int *dim = solver->dim_local, done;
653 
654  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
655  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
656  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
657 
661 
662  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
663  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
664 
665  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
666 
667  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
668  while (!done) {
669  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
670  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
671 
672  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
673  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
674  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
675  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
676  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
677  int q = p*_MODEL_NVARS_;
678 
679  /* Modified Rusanov's upwinding scheme */
680 
681  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
682  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
683  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
684  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
685  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
686 
688  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
689  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
690 
691  double c, vel[_MODEL_NDIMS_], rho,E,P;
692 
693  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
694  c = sqrt(param->gamma*P/rho);
695  double alphaL = c + absolute(vel[dir]);
696  double betaL = absolute(vel[dir]);
697 
698  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
699  c = sqrt(param->gamma*P/rho);
700  double alphaR = c + absolute(vel[dir]);
701  double betaR = absolute(vel[dir]);
702 
703  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
704  c = sqrt(param->gamma*P/rho);
705  double alphaavg = c + absolute(vel[dir]);
706  double betaavg = absolute(vel[dir]);
707 
708  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
709  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
710  double beta = kappa*max3(betaL,betaR,betaavg);
711 
713  if (dir == _XDIR_) {
714  D[0] = beta;
715  D[6] = alpha;
716  D[12] = beta;
717  D[18] = beta;
718  D[24] = alpha;
719  } else if (dir == _YDIR_) {
720  D[0] = beta;
721  D[6] = beta;
722  D[12] = alpha;
723  D[18] = beta;
724  D[24] = alpha;
725  } else if (dir == _ZDIR_) {
726  D[0] = beta;
727  D[6] = beta;
728  D[12] = beta;
729  D[18] = alpha;
730  D[24] = alpha;
731  }
732  MatMult5 (_MODEL_NVARS_,DL,D,L);
733  MatMult5 (_MODEL_NVARS_,modA,R,DL);
734  MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
735 
736  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
737  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
738  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
739  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
740  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
741  }
742  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
743  }
744 
745  return(0);
746 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _ArrayCopy1D_(x, y, size)
#define max(a, b)
Definition: math_ops.h:18
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DUpwinddFRusanovModified()

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

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

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 756 of file NavierStokes3DUpwind.c.

767 {
768  HyPar *solver = (HyPar*) s;
769  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
770  int *dim = solver->dim_local, done;
771  double *uref = param->solution;
772 
773  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
774  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
775  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
776 
780 
781  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
782  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
783 
784  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
785 
786  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
787  while (!done) {
788  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
789  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
790 
791  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
792  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
793  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
794  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
795  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
796  int q = p*_MODEL_NVARS_;
797 
798  /* Modified Rusanov's upwinding scheme */
799 
800  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
801  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
802  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
803  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
804  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
805 
807  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
808  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
809 
810  double c, vel[_MODEL_NDIMS_], rho,E,P;
811 
812  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
813  c = sqrt(param->gamma*P/rho);
814  double alphaL = c + absolute(vel[dir]);
815 
816  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
817  c = sqrt(param->gamma*P/rho);
818  double alphaR = c + absolute(vel[dir]);
819 
820  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
821  c = sqrt(param->gamma*P/rho);
822  double alphaavg = c + absolute(vel[dir]);
823 
824  double kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
825  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
826 
828  if (dir == _XDIR_) {
829  D[6] = alpha;
830  D[24] = alpha;
831  } else if (dir == _YDIR_) {
832  D[12] = alpha;
833  D[24] = alpha;
834  } else if (dir == _ZDIR_) {
835  D[18] = alpha;
836  D[24] = alpha;
837  }
838  MatMult5 (_MODEL_NVARS_,DL,D,L);
839  MatMult5 (_MODEL_NVARS_,modA,R,DL);
840  MatVecMult5 (_MODEL_NVARS_,udiss,modA,udiff);
841 
842  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
843  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
844  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
845  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
846  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
847  }
848  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
849  }
850 
851  return(0);
852 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _ArrayCopy1D_(x, y, size)
#define max(a, b)
Definition: math_ops.h:18
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DUpwindFdFRusanovModified()

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

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

  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent solution time

Definition at line 862 of file NavierStokes3DUpwind.c.

873 {
874  HyPar *solver = (HyPar*) s;
875  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
876  int *dim = solver->dim_local, done;
877  double *uref = param->solution;
878 
879  static int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
880  bounds_outer[0] = dim[0]; bounds_outer[1] = dim[1]; bounds_outer[2] = dim[2]; bounds_outer[dir] = 1;
881  bounds_inter[0] = dim[0]; bounds_inter[1] = dim[1]; bounds_inter[2] = dim[2]; bounds_inter[dir]++;
882 
886 
887  static int indexL[_MODEL_NDIMS_], indexR[_MODEL_NDIMS_],
888  index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_];
889 
890  static double udiff[_MODEL_NVARS_],uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
891 
892  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
893  while (!done) {
894  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
895  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
896 
897  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
898  _ArrayCopy1D_(index_inter,indexL,_MODEL_NDIMS_); indexL[dir]--;
899  _ArrayCopy1D_(index_inter,indexR,_MODEL_NDIMS_);
900  int pL; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexL,solver->ghosts,pL);
901  int pR; _ArrayIndex1D_(_MODEL_NDIMS_,dim,indexR,solver->ghosts,pR);
902  int q = p*_MODEL_NVARS_;
903  double udiff[_MODEL_NVARS_], udiss[_MODEL_NVARS_], uavg[_MODEL_NVARS_],
904  udiss_total[_MODEL_NVARS_],udiss_acoustic[_MODEL_NVARS_];
905 
906  /* Modified Rusanov's upwinding scheme */
907  /* Modified Rusanov's upwinding scheme */
908  double c, vel[_MODEL_NDIMS_], rho,E,P, alphaL, alphaR, alphaavg,
909  betaL, betaR, betaavg, alpha, beta,
910  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
911 
912  udiff[0] = 0.5 * (uR[q+0] - uL[q+0]);
913  udiff[1] = 0.5 * (uR[q+1] - uL[q+1]);
914  udiff[2] = 0.5 * (uR[q+2] - uL[q+2]);
915  udiff[3] = 0.5 * (uR[q+3] - uL[q+3]);
916  udiff[4] = 0.5 * (uR[q+4] - uL[q+4]);
917 
918  /* Compute total dissipation */
920  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
921  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
922 
923  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
924  c = sqrt(param->gamma*P/rho);
925  alphaL = c + absolute(vel[dir]);
926  betaL = absolute(vel[dir]);
927 
928  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
929  c = sqrt(param->gamma*P/rho);
930  alphaR = c + absolute(vel[dir]);
931  betaR = absolute(vel[dir]);
932 
933  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
934  c = sqrt(param->gamma*P/rho);
935  alphaavg = c + absolute(vel[dir]);
936  betaavg = absolute(vel[dir]);
937 
938  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
939  alpha = kappa*max3(alphaL,alphaR,alphaavg);
940  beta = kappa*max3(betaL,betaR,betaavg);
941 
943  if (dir == _XDIR_) {
944  D[0] = beta;
945  D[6] = alpha;
946  D[12] = beta;
947  D[18] = beta;
948  D[24] = alpha;
949  } else if (dir == _YDIR_) {
950  D[0] = beta;
951  D[6] = beta;
952  D[12] = alpha;
953  D[18] = beta;
954  D[24] = alpha;
955  } else if (dir == _ZDIR_) {
956  D[0] = beta;
957  D[6] = beta;
958  D[12] = beta;
959  D[18] = alpha;
960  D[24] = alpha;
961  }
962  MatMult5 (_MODEL_NVARS_,DL,D,L);
963  MatMult5 (_MODEL_NVARS_,modA,R,DL);
964  MatVecMult5 (_MODEL_NVARS_,udiss_total,modA,udiff);
965 
966  /* Compute dissipation for the linearized acoustic modes */
968  _NavierStokes3DLeftEigenvectors_ (uavg,dummy,L,param->gamma,dir);
969  _NavierStokes3DRightEigenvectors_ (uavg,dummy,R,param->gamma,dir);
970 
971  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pL),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
972  c = sqrt(param->gamma*P/rho);
973  alphaL = c + absolute(vel[dir]);
974 
975  _NavierStokes3DGetFlowVar_((uref+_MODEL_NVARS_*pR),_NavierStokes3D_stride_,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
976  c = sqrt(param->gamma*P/rho);
977  alphaR = c + absolute(vel[dir]);
978 
979  _NavierStokes3DGetFlowVar_(uavg,dummy,rho,vel[0],vel[1],vel[2],E,P,param->gamma);
980  c = sqrt(param->gamma*P/rho);
981  alphaavg = c + absolute(vel[dir]);
982 
983  kappa = max(param->grav_field_g[pL],param->grav_field_g[pR]);
984  alpha = kappa*max3(alphaL,alphaR,alphaavg);
985 
987  if (dir == _XDIR_) {
988  D[6] = alpha;
989  D[24] = alpha;
990  } else if (dir == _YDIR_) {
991  D[12] = alpha;
992  D[24] = alpha;
993  } else if (dir == _ZDIR_) {
994  D[18] = alpha;
995  D[24] = alpha;
996  }
997  MatMult5 (_MODEL_NVARS_,DL,D,L);
998  MatMult5 (_MODEL_NVARS_,modA,R,DL);
999  MatVecMult5 (_MODEL_NVARS_,udiss_acoustic,modA,udiff);
1000 
1001  /* Compute dissipation for the entropy modes */
1002  _ArraySubtract1D_(udiss,udiss_total,udiss_acoustic,_MODEL_NVARS_);
1003 
1004  fI[q+0] = 0.5*(fL[q+0]+fR[q+0])-udiss[0];
1005  fI[q+1] = 0.5*(fL[q+1]+fR[q+1])-udiss[1];
1006  fI[q+2] = 0.5*(fL[q+2]+fR[q+2])-udiss[2];
1007  fI[q+3] = 0.5*(fL[q+3]+fR[q+3])-udiss[3];
1008  fI[q+4] = 0.5*(fL[q+4]+fR[q+4])-udiss[4];
1009  }
1010  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
1011  }
1012 
1013  return(0);
1014 }
#define max3(a, b, c)
Definition: math_ops.h:27
#define absolute(a)
Definition: math_ops.h:32
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
#define _NavierStokes3DRoeAverage_(uavg, stride, uL, uR, gamma)
static const int dummy
#define _ArraySubtract1D_(x, a, b, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1D_(N, imax, i, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define MatVecMult5(N, y, A, x)
#define _ArrayCopy1D_(x, y, size)
#define max(a, b)
Definition: math_ops.h:18
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DGravityField()

int NavierStokes3DGravityField ( void *  s,
void *  m 
)

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

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

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

References:

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

Definition at line 33 of file NavierStokes3DGravityField.c.

37 {
38  HyPar *solver = (HyPar*) s;
39  MPIVariables *mpi = (MPIVariables*) m;
40  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
41 
42  double *f = param->grav_field_f;
43  double *g = param->grav_field_g;
44  int *dim = solver->dim_local;
45  int ghosts = solver->ghosts;
46  int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_],
47  offset[_MODEL_NDIMS_], d, done;
48 
49  /* set bounds for array index to include ghost points */
50  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_);
51  for (d=0; d<_MODEL_NDIMS_; d++) bounds[d] += 2*ghosts;
52  /* set offset such that index is compatible with ghost point arrangement */
53  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
54 
55  double p0 = param->p0;
56  double rho0 = param->rho0;
57  double RT = p0 / rho0;
58  double gamma= param->gamma;
59  double R = param->R;
60  double Cp = gamma*R / (gamma-1.0);
61  double T0 = p0 / (rho0 * R);
62  double gx = param->grav_x;
63  double gy = param->grav_y;
64  double gz = param->grav_z;
65  double Nbv = param->N_bv;
66 
67  /* set the value of the gravity field */
68  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
69  if (param->HB == 1) {
70  while (!done) {
71  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
72  double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
73  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
74  double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
75  f[p] = exp( (gx*xcoord+gy*ycoord+gz*zcoord)/RT);
76  g[p] = exp(-(gx*xcoord+gy*ycoord+gz*zcoord)/RT);
77  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
78  }
79  } else if (param->HB == 2) {
80  while (!done) {
81  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
82  double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
83  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
84  double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
85  f[p] = raiseto((1.0-(gx*xcoord+gy*ycoord+gz*zcoord)/(Cp*T0)), (-1.0 /(gamma-1.0)));
86  g[p] = raiseto((1.0-(gx*xcoord+gy*ycoord+gz*zcoord)/(Cp*T0)), (gamma/(gamma-1.0)));
87  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
88  }
89  } else if (param->HB == 3) {
90  while (!done) {
91  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
92  double zcoord; _GetCoordinate_(_ZDIR_,index[_ZDIR_]-ghosts,dim,ghosts,solver->x,zcoord);
93  if ((gx != 0) || (gy != 0)) {
94  if (!mpi->rank) {
95  fprintf(stderr,"Error in NavierStokes3DGravityField(): HB = 3 is implemented only for ");
96  fprintf(stderr,"gravity force along the z-coordinate.\n");
97  }
98  return(1);
99  }
100  double Pexner = 1 + ((gz*gz)/(Cp*T0*Nbv*Nbv)) * (exp(-(Nbv*Nbv/gz)*zcoord)-1.0);
101  f[p] = raiseto(Pexner, (-1.0 /(gamma-1.0))) * exp(Nbv*Nbv*zcoord/gz);
102  g[p] = raiseto(Pexner, (gamma/(gamma-1.0)));
103  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
104  }
105  } else {
106  while (!done) {
107  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
108  f[p] = g[p] = 1.0;
109  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
110  }
111  }
112 
113  /* A sensible simulation will not specify peridic boundary conditions
114  * along a direction in which gravity acts.
115  * Gravity will be zero along the dimension periodic BCs are specified,
116  * so the value of the gravity potential will be the same along those
117  * grid lines.
118  * Thus, for both these cases, extrapolate the gravity field at the
119  * boundaries */
120  int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_];
121  for (d = 0; d < _MODEL_NDIMS_; d++) {
122  /* left boundary */
123  if (!mpi->ip[d]) {
124  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
125  _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
126  done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
127  while (!done) {
128  _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = ghosts-1-indexb[d];
129  int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
130  int p2; _ArrayIndex1D_ (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
131  f[p1] = f[p2];
132  g[p1] = g[p2];
133  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
134  }
135  }
136  /* right boundary */
137  if (mpi->ip[d] == mpi->iproc[d]-1) {
138  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
139  _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
140  done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
141  while (!done) {
142  _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = dim[d]-1-indexb[d];
143  int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
144  int p2; _ArrayIndex1D_ (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
145  f[p1] = f[p2];
146  g[p1] = g[p2];
147  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
148  }
149  }
150  }
151 
152  return(0);
153 }
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
double * x
Definition: hypar.h:107
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define raiseto(x, a)
Definition: math_ops.h:37
#define _ArrayCopy1D_(x, y, size)
double * grav_field_f

◆ NavierStokes3DModifiedSolution()

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

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

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

where

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

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

References:

  • Ghosh, D., Constantinescu, E.M., Well-Balanced Formulation of Gravitational Source Terms for Conservative Finite-Difference Atmospheric Flow Solvers, AIAA Paper 2015-2889, 7th AIAA Atmospheric and Space Environments Conference, June 22-26, 2015, Dallas, TX, http://dx.doi.org/10.2514/6.2015-2889
  • Ghosh, D., Constantinescu, E.M., A Well-Balanced, Conservative Finite-Difference Algorithm for Atmospheric Flows, AIAA Journal, 54 (4), 2016, pp. 1370-1385, http://dx.doi.org/10.2514/1.J054580.
Parameters
uCArray to hold the computed modified solution
uSolution vector array
dspatial dimension (not used)
sSolver object of type HyPar
mMPI object of time MPIVariables
waqtCurrent simulation time

Definition at line 31 of file NavierStokes3DModifiedSolution.c.

39 {
40  HyPar *solver = (HyPar*) s;
41  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
42 
43  int ghosts = solver->ghosts;
44  int *dim = solver->dim_local;
45  int ndims = solver->ndims;
46  int index[ndims], bounds[ndims], offset[ndims];
47 
48  /* set bounds for array index to include ghost points */
49  _ArrayCopy1D_(dim,bounds,ndims);
50  int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
51 
52  /* set offset such that index is compatible with ghost point arrangement */
53  _ArraySetValue_(offset,ndims,-ghosts);
54 
55  int done = 0; _ArraySetValue_(index,ndims,0);
56  double inv_gamma_m1 = 1.0 / (param->gamma-1.0);
57  while (!done) {
58  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
59  double rho, uvel, vvel, wvel, E, P;
60  _NavierStokes3DGetFlowVar_((u+_MODEL_NVARS_*p),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,E,P,param->gamma);
61  uC[_MODEL_NVARS_*p+0] = u[_MODEL_NVARS_*p+0] * param->grav_field_f[p];
62  uC[_MODEL_NVARS_*p+1] = u[_MODEL_NVARS_*p+1] * param->grav_field_f[p];
63  uC[_MODEL_NVARS_*p+2] = u[_MODEL_NVARS_*p+2] * param->grav_field_f[p];
64  uC[_MODEL_NVARS_*p+3] = u[_MODEL_NVARS_*p+3] * param->grav_field_f[p];
65  uC[_MODEL_NVARS_*p+4] = (P*inv_gamma_m1)*(1.0/param->grav_field_g[p]) + (0.5*rho*(uvel*uvel+vvel*vvel+wvel*wvel))*param->grav_field_f[p];
66  _ArrayIncrementIndex_(ndims,bounds,index,done);
67  }
68 
69  return(0);
70 }
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _ArrayCopy1D_(x, y, size)
double * grav_field_f
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DIBAdiabatic()

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

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

Parameters
sSolver object of type HyPar
mSolver object of type HyPar
uArray with the solution vector
tCurrent simulation time

Definition at line 19 of file NavierStokes3DImmersedBoundary.c.

24 {
25  HyPar *solver = (HyPar*) s;
26  MPIVariables *mpi = (MPIVariables*) m;
27  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
28  IBNode *boundary = IB->boundary;
29  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
30  static double v[_MODEL_NVARS_];
31  int n, j, k, nb = IB->n_boundary_nodes;
32 
33  if (!solver->flag_ib) return(0);
34 
35  /* Ideally, this shouldn't be here - But this function is called everywhere
36  (through ApplyIBConditions()) *before* MPIExchangeBoundariesnD is called! */
38 
39  double inv_gamma_m1 = 1.0 / (param->gamma - 1.0);
40 
41  double ramp_fac = 1.0;
42  if (param->t_ib_ramp > 0) {
43  double x = t/param->t_ib_ramp;
44  if (!strcmp(param->ib_ramp_type,_IB_RAMP_LINEAR_)) {
45  ramp_fac = x;
46  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_SMOOTHEDSLAB_)) {
47  double a = 0.0;
48  double b = 1.0;
49  double c = 0.5;
50  double r = param->t_ib_ramp/param->t_ib_width;
51  ramp_fac = (a*exp(c*r)+b*exp(r*x))/(exp(c*r)+exp(r*x));
52  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_DISABLE_)) {
53  ramp_fac = 0.0;
54  } else {
55  fprintf(stderr,"Error in NavierStokes3DImmersedBoundary():\n");
56  fprintf(stderr," Ramp type %s not recognized.\n", param->ib_ramp_type);
57  return 1;
58  }
59  }
60 
61  for (n=0; n<nb; n++) {
62 
63  int node_index = boundary[n].p;
64  double *alpha = &(boundary[n].interp_coeffs[0]);
65  int *nodes = &(boundary[n].interp_nodes[0]);
66  double factor = boundary[n].surface_distance / boundary[n].interp_node_distance;
67 
69  for (j=0; j<_IB_NNODES_; j++) {
70  for (k=0; k<_MODEL_NVARS_; k++) {
71  v[k] += ( alpha[j] * u[_MODEL_NVARS_*nodes[j]+k] );
72  }
73  }
74 
75  double rho, uvel, vvel, wvel, energy, pressure;
76  _NavierStokes3DGetFlowVar_(v,_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,param->gamma);
77 
78  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
81  rho_gpt,
82  uvel_gpt,
83  vvel_gpt,
84  wvel_gpt,
85  energy_gpt,
86  pressure_gpt,
87  param->gamma );
88 
89  double rho_ib_target, uvel_ib_target, vvel_ib_target, wvel_ib_target, pressure_ib_target;
90  rho_ib_target = rho;
91  pressure_ib_target = pressure;
92  uvel_ib_target = -uvel * factor;
93  vvel_ib_target = -vvel * factor;
94  wvel_ib_target = -wvel * factor;
95 
96  double rho_ib, uvel_ib, vvel_ib, wvel_ib, energy_ib, pressure_ib;
97  rho_ib = ramp_fac * rho_ib_target + (1.0-ramp_fac) * rho_gpt;
98  pressure_ib = ramp_fac * pressure_ib_target + (1.0-ramp_fac) * pressure_gpt;
99  uvel_ib = ramp_fac * uvel_ib_target + (1.0-ramp_fac) * uvel_gpt;
100  vvel_ib = ramp_fac * vvel_ib_target + (1.0-ramp_fac) * vvel_gpt;
101  wvel_ib = ramp_fac * wvel_ib_target + (1.0-ramp_fac) * wvel_gpt;
102  energy_ib = inv_gamma_m1*pressure_ib
103  + 0.5*rho_ib*(uvel_ib*uvel_ib+vvel_ib*vvel_ib+wvel_ib*wvel_ib);
104 
105  u[_MODEL_NVARS_*node_index+0] = rho_ib;
106  u[_MODEL_NVARS_*node_index+1] = rho_ib * uvel_ib;
107  u[_MODEL_NVARS_*node_index+2] = rho_ib * vvel_ib;
108  u[_MODEL_NVARS_*node_index+3] = rho_ib * wvel_ib;
109  u[_MODEL_NVARS_*node_index+4] = energy_ib;
110  }
111 
112  return 0;
113 }
#define _IB_NNODES_
Structure defining an immersed boundary node.
Structure containing variables for immersed boundary implementation.
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
void * ib
Definition: hypar.h:443
int flag_ib
Definition: hypar.h:441
#define _IB_RAMP_SMOOTHEDSLAB_
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _IB_RAMP_DISABLE_
#define _IB_RAMP_LINEAR_
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DIBIsothermal()

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

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

Parameters
sSolver object of type HyPar
mSolver object of type HyPar
uArray with the solution vector
tCurrent simulation time

Definition at line 119 of file NavierStokes3DImmersedBoundary.c.

124 {
125  HyPar *solver = (HyPar*) s;
126  MPIVariables *mpi = (MPIVariables*) m;
127  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
128  IBNode *boundary = IB->boundary;
129  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
130  static double v[_MODEL_NVARS_];
131  int n, j, k, nb = IB->n_boundary_nodes;
132 
133  if (!solver->flag_ib) return(0);
134 
135  /* Ideally, this shouldn't be here - But this function is called everywhere
136  (through ApplyIBConditions()) *before* MPIExchangeBoundariesnD is called! */
138 
139  double inv_gamma_m1 = 1.0 / (param->gamma - 1.0);
140 
141  double ramp_fac = 1.0;
142  if (param->t_ib_ramp > 0) {
143  double x = t/param->t_ib_ramp;
144  if (!strcmp(param->ib_ramp_type,_IB_RAMP_LINEAR_)) {
145  ramp_fac = x;
146  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_SMOOTHEDSLAB_)) {
147  double a = 0.0;
148  double b = 1.0;
149  double c = 0.5;
150  double r = param->t_ib_ramp/param->t_ib_width;
151  ramp_fac = (a*exp(c*r)+b*exp(r*x))/(exp(c*r)+exp(r*x));
152  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_DISABLE_)) {
153  ramp_fac = 0.0;
154  } else {
155  fprintf(stderr,"Error in NavierStokes3DImmersedBoundary():\n");
156  fprintf(stderr," Ramp type %s not recognized.\n", param->ib_ramp_type);
157  return 1;
158  }
159  }
160 
161  for (n=0; n<nb; n++) {
162 
163  int node_index = boundary[n].p;
164  double *alpha = &(boundary[n].interp_coeffs[0]);
165  int *nodes = &(boundary[n].interp_nodes[0]);
166  double factor = boundary[n].surface_distance / boundary[n].interp_node_distance;
167 
169  for (j=0; j<_IB_NNODES_; j++) {
170  for (k=0; k<_MODEL_NVARS_; k++) {
171  v[k] += ( alpha[j] * u[_MODEL_NVARS_*nodes[j]+k] );
172  }
173  }
174 
175  double rho, uvel, vvel, wvel, energy, pressure, temperature;
176  _NavierStokes3DGetFlowVar_(v,_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,param->gamma);
177  temperature = pressure / rho;
178 
179  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt, temperature_gpt;
182  rho_gpt,
183  uvel_gpt,
184  vvel_gpt,
185  wvel_gpt,
186  energy_gpt,
187  pressure_gpt,
188  param->gamma );
189  temperature_gpt = pressure_gpt / rho_gpt;
190 
191  double rho_ib_target,
192  uvel_ib_target, vvel_ib_target, wvel_ib_target,
193  pressure_ib_target,
194  temperature_ib_target;
195  temperature_ib_target = (1.0+factor)*param->T_ib_wall - factor * temperature;
196  if ( (temperature_ib_target < param->T_ib_wall/param->ib_T_tol)
197  || (temperature_ib_target > param->T_ib_wall*param->ib_T_tol) ) {
198  temperature_ib_target = param->T_ib_wall;
199  }
200  pressure_ib_target = pressure;
201  rho_ib_target = pressure_ib_target / temperature_ib_target;
202  uvel_ib_target = - factor * uvel;
203  vvel_ib_target = - factor * vvel;
204  wvel_ib_target = - factor * wvel;
205 
206  double rho_ib, uvel_ib, vvel_ib, wvel_ib, energy_ib, pressure_ib;
207  rho_ib = ramp_fac * rho_ib_target + (1.0-ramp_fac) * rho_gpt;
208  pressure_ib = ramp_fac * pressure_ib_target + (1.0-ramp_fac) * pressure_gpt;
209  uvel_ib = ramp_fac * uvel_ib_target + (1.0-ramp_fac) * uvel_gpt;
210  vvel_ib = ramp_fac * vvel_ib_target + (1.0-ramp_fac) * vvel_gpt;
211  wvel_ib = ramp_fac * wvel_ib_target + (1.0-ramp_fac) * wvel_gpt;
212  energy_ib = inv_gamma_m1*pressure_ib
213  + 0.5*rho_ib*(uvel_ib*uvel_ib+vvel_ib*vvel_ib+wvel_ib*wvel_ib);
214 
215  u[_MODEL_NVARS_*node_index+0] = rho_ib;
216  u[_MODEL_NVARS_*node_index+1] = rho_ib * uvel_ib;
217  u[_MODEL_NVARS_*node_index+2] = rho_ib * vvel_ib;
218  u[_MODEL_NVARS_*node_index+3] = rho_ib * wvel_ib;
219  u[_MODEL_NVARS_*node_index+4] = energy_ib;
220  }
221 
222  return(0);
223 }
#define _IB_NNODES_
Structure defining an immersed boundary node.
Structure containing variables for immersed boundary implementation.
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
void * ib
Definition: hypar.h:443
int flag_ib
Definition: hypar.h:441
#define _IB_RAMP_SMOOTHEDSLAB_
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MODEL_NDIMS_
Definition: euler1d.h:56
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _IB_RAMP_DISABLE_
#define _IB_RAMP_LINEAR_
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define _MODEL_NVARS_
Definition: euler1d.h:58
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)

◆ NavierStokes3DPreStep()

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

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

  • The solution at the beginning of the step is copied into NavierStokes3D::solution for the linearized flux partitioning.
  • The linearized "fast" Jacobian (representing the acoustic modes) NavierStokes3D::fast_jac is computed as:

    \begin{equation} A_f\left({\bf u}\right) = R\left({\bf u}\right)\Lambda_f\left({\bf u}\right)L\left({\bf u}\right) \end{equation}

    where \(R\) and \(L\) are the matrices of right and left eigenvectors, and,

    \begin{equation} \Lambda_f = diag\left[0,0,0,u+c,u-c \right] \end{equation}

    with \(c\) being the speed of sound.

    Reference:
  • Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.
Parameters
uSolution vector
sSolver object of type HyPar
mMPI object of type MPIVariables
waqtCurrent simulation time

Definition at line 32 of file NavierStokes3DPreStep.c.

38 {
39  HyPar *solver = (HyPar*) s;
40  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
41  int *dim = solver->dim_local;
42  int ghosts = solver->ghosts, dir, p;
43  double *A;
44 
45  static const int ndims = _MODEL_NDIMS_;
46  static const int JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
47  static int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
50 
51  /* set bounds for array index to include ghost points */
52  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
53  /* set offset such that index is compatible with ghost point arrangement */
54  _ArraySetValue_(offset,ndims,-ghosts);
55  /* copy the solution to act as a reference for linearization */
57 
58  int done = 0; _ArraySetValue_(index,ndims,0);
59  while (!done) {
60  _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
61  int q = _MODEL_NVARS_*p;
62 
63  dir = _XDIR_;
64  A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
65  /* get the eigenvalues, and left & right eigenvectors */
69  /* remove the entropy modes (corresponding to eigenvalues u) */
70  D[0] = D[12] = D[18] = 0.0;
71  /* assemble the Jacobian */
72  MatMult5(_MODEL_NVARS_,DL,D,L );
73  MatMult5(_MODEL_NVARS_,A ,R,DL);
74 
75  dir = _YDIR_;
76  A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
77  /* get the eigenvalues, and left & right eigenvectors */
81  /* remove the entropy modes (corresponding to eigenvalues v) */
82  D[0] = D[6] = D[18] = 0.0;
83  /* assemble the Jacobian */
84  MatMult5(_MODEL_NVARS_,DL,D,L );
85  MatMult5(_MODEL_NVARS_,A ,R,DL);
86 
87  dir = _ZDIR_;
88  A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
89  /* get the eigenvalues, and left & right eigenvectors */
93  /* remove the entropy modes (corresponding to eigenvalues v) */
94  D[0] = D[6] = D[12] = 0.0;
95  /* assemble the Jacobian */
96  MatMult5(_MODEL_NVARS_,DL,D,L );
97  MatMult5(_MODEL_NVARS_,A ,R,DL);
98 
99  _ArrayIncrementIndex_(ndims,bounds,index,done);
100  }
101  return(0);
102 }
#define _NavierStokes3DRightEigenvectors_(u, stride, R, ga, dir)
#define _NavierStokes3DLeftEigenvectors_(u, stride, L, ga, dir)
#define _ZDIR_
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
#define _ArrayAddCopy1D_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
#define _NavierStokes3DEigenvalues_(u, stride, D, gamma, dir)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
static const int _NavierStokes3D_stride_
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
#define MatMult5(N, A, X, Y)
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArrayCopy1D_(x, y, size)

◆ NavierStokes3DIBForces()

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

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

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
a_tCurrent simulation time

Definition at line 264 of file NavierStokes3DIBForces.c.

267 {
268  HyPar *solver = (HyPar*) s;
269  MPIVariables *mpi = (MPIVariables*) m;
270  NavierStokes3D *physics = (NavierStokes3D*) solver->physics;
271  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
272  int ierr;
273 
274  if (!solver->flag_ib) return(0);
275 
276  int npts = solver->npoints_local_wghosts;
277 
278  double* pressure = (double*) calloc (npts, sizeof(double));
279  ierr = NavierStokes3DComputePressure(pressure, solver->u, solver);
280  if (ierr) {
281  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
282  fprintf(stderr," NavierStokes3DComputePressure() returned with error.\n");
283  return 1;
284  }
285  double* temperature = (double*) calloc(npts, sizeof(double));
286  ierr = NavierStokes3DComputeTemperature(temperature, solver->u, solver);
287  if (ierr) {
288  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
289  fprintf(stderr," NavierStokes3DComputeTemperature() returned with error.\n");
290  return 1;
291  }
292 
293  /* Compute surface pressure */
294  double* p_surface = NULL;
295  ierr = IBComputeFacetVar(solver, mpi, pressure, 1, &p_surface);
296  if (ierr) {
297  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
298  fprintf(stderr," IBComputeFacetVar() returned with error.\n");
299  return 1;
300  }
301  /* Compute surface temperature */
302  double* T_surface = NULL;
303  ierr = IBComputeFacetVar(solver, mpi, temperature, 1, &T_surface);
304  if (ierr) {
305  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
306  fprintf(stderr," IBComputeFacetVar() returned with error.\n");
307  return 1;
308  }
309  /* Compute normal surface pressure gradient */
310  double *ngrad_p_surface = NULL;
311  ierr = IBComputeNormalGradient(solver, mpi, pressure, 1, &ngrad_p_surface);
312  if (ierr) {
313  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
314  fprintf(stderr," IBComputeNormalGradient() returned with error.\n");
315  return 1;
316  }
317  /* Compute normal temperature gradient */
318  double *ngrad_T_surface = NULL;
319  ierr = IBComputeNormalGradient(solver, mpi, temperature, 1, &ngrad_T_surface);
320  if (ierr) {
321  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
322  fprintf(stderr," IBComputeNormalGradient() returned with error.\n");
323  return 1;
324  }
325  /* Compute shear forces */
326  double *shear = NULL;
327  ierr = ComputeShear(solver, mpi, solver->u, &shear);
328  if (ierr) {
329  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
330  fprintf(stderr," ComputeShear() returned with error.\n");
331  return 1;
332  }
333 
334  char surface_filename[_MAX_STRING_SIZE_] = "surface";
335  if (solver->nsims == 1) {
336  if (!strcmp(solver->op_overwrite,"no")) {
337  strcat(surface_filename,solver->filename_index);
338  }
339  } else {
340  char index[_MAX_STRING_SIZE_];
341  GetStringFromInteger(solver->my_idx, index, (int)log10(solver->nsims)+1);
342  strcat(surface_filename, "_");
343  strcat(surface_filename, index);
344  strcat(surface_filename, "_");
345  }
346  strcat(surface_filename,".dat");
347  if (!mpi->rank) {
348  printf("Writing immersed body surface data file %s.\n",surface_filename);
349  }
350  ierr = WriteSurfaceData( mpi,
351  IB,
352  p_surface,
353  T_surface,
354  ngrad_p_surface,
355  ngrad_T_surface,
356  shear,
357  surface_filename );
358  if (ierr) {
359  fprintf(stderr,"Error in NavierStokes3DIBForces()\n");
360  fprintf(stderr," WriteSurfaceData() returned with error\n");
361  return 1;
362  }
363 
364  free(pressure);
365  free(temperature);
366  if (p_surface) free(p_surface);
367  if (T_surface) free(T_surface);
368  if (ngrad_p_surface) free(ngrad_p_surface);
369  if (ngrad_T_surface) free(ngrad_T_surface);
370  if (shear) free(shear);
371 
372  return 0;
373 }
Structure containing variables for immersed boundary implementation.
char * filename_index
Definition: hypar.h:197
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
int nsims
Definition: hypar.h:64
double * u
Definition: hypar.h:116
void * ib
Definition: hypar.h:443
int NavierStokes3DComputeTemperature(double *, const double *const, void *)
int flag_ib
Definition: hypar.h:441
int NavierStokes3DComputePressure(double *, const double *const, void *)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _MAX_STRING_SIZE_
Definition: basic.h:14
static int WriteSurfaceData(void *m, void *ib, const double *const p_surface, const double *const T_surface, const double *const ngrad_p_surface, const double *const ngrad_T_surface, const double *const shear, char *filename)
char op_overwrite[_MAX_STRING_SIZE_]
Definition: hypar.h:191
static int ComputeShear(void *s, void *m, const double *const u, double **const sf)
void GetStringFromInteger(int, char *, int)
int my_idx
Definition: hypar.h:61
int IBComputeNormalGradient(void *, void *, const double *const, int, double **const)
void * physics
Definition: hypar.h:266
Structure of MPI-related variables.
int npoints_local_wghosts
Definition: hypar.h:42
int IBComputeFacetVar(void *, void *, const double *const, int, double **const)

◆ gpuNavierStokes3DInitialize()

int gpuNavierStokes3DInitialize ( void *  s,
void *  m 
)

Initialize GPU-related arrays.

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables

Definition at line 10 of file NavierStokes3DInitialize_GPU.c.

14 {
15  HyPar *solver = (HyPar*) s;
16  NavierStokes3D *physics = (NavierStokes3D*) solver->physics;
17 
18  int *dim = solver->dim_local;
19  int ghosts = solver->ghosts;
20  int d, size = 1; for (d = 0; d <_MODEL_NDIMS_; d++) size *= (dim[d] + 2*ghosts);
21 
22  gpuMalloc((void**)&physics->gpu_Q, size*_MODEL_NVARS_*sizeof(double));
23  gpuMalloc((void**)&physics->gpu_QDerivX, size*_MODEL_NVARS_*sizeof(double));
24  gpuMalloc((void**)&physics->gpu_QDerivY, size*_MODEL_NVARS_*sizeof(double));
25  gpuMalloc((void**)&physics->gpu_QDerivZ, size*_MODEL_NVARS_*sizeof(double));
26  gpuMalloc((void**)&physics->gpu_FViscous, size*_MODEL_NVARS_*sizeof(double));
27  gpuMalloc((void**)&physics->gpu_FDeriv, size*_MODEL_NVARS_*sizeof(double));
28  gpuMalloc((void**)&physics->gpu_grav_field_f, size*sizeof(double));
29  gpuMalloc((void**)&physics->gpu_grav_field_g, size*sizeof(double));
30  gpuMalloc((void**)&physics->gpu_fast_jac, _MODEL_NDIMS_*size*_MODEL_NVARS_*_MODEL_NVARS_*sizeof(double));
31  gpuMalloc((void**)&physics->gpu_solution, size*_MODEL_NVARS_*sizeof(double));
32  gpuMemset(physics->gpu_Q, 0, size*_MODEL_NVARS_*sizeof(double));
33  gpuMemset(physics->gpu_QDerivX, 0, size*_MODEL_NVARS_*sizeof(double));
34  gpuMemset(physics->gpu_QDerivY, 0, size*_MODEL_NVARS_*sizeof(double));
35  gpuMemset(physics->gpu_QDerivZ, 0, size*_MODEL_NVARS_*sizeof(double));
36  gpuMemset(physics->gpu_FViscous, 0, size*_MODEL_NVARS_*sizeof(double));
37  gpuMemset(physics->gpu_FDeriv, 0, size*_MODEL_NVARS_*sizeof(double));
38  gpuMemcpy(physics->gpu_grav_field_f, physics->grav_field_f, size*sizeof(double), gpuMemcpyHostToDevice);
39  gpuMemcpy(physics->gpu_grav_field_g, physics->grav_field_g, size*sizeof(double), gpuMemcpyHostToDevice);
40  gpuMemset(physics->gpu_fast_jac, 0, _MODEL_NDIMS_*size*_MODEL_NVARS_*_MODEL_NVARS_*sizeof(double));
41  gpuMemset(physics->gpu_solution, 0, size*_MODEL_NVARS_*sizeof(double));
42 
43  return (0);
44 }
double * gpu_QDerivX
void gpuMemset(void *, int, size_t)
double * gpu_QDerivY
double * gpu_grav_field_f
double * gpu_solution
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
double * grav_field_g
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
void gpuMalloc(void **, size_t)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
double * gpu_fast_jac
int * dim_local
Definition: hypar.h:37
void gpuMemcpy(void *, const void *, size_t, enum gpuMemcpyKind)
double * gpu_grav_field_g
void * physics
Definition: hypar.h:266
double * gpu_FViscous
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
double * gpu_QDerivZ
double * gpu_FDeriv
double * grav_field_f

◆ gpuNavierStokes3DFlux()

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

◆ gpuNavierStokes3DParabolicFunction()

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

◆ gpuNavierStokes3DSource()

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

◆ gpuNavierStokes3DModifiedSolution()

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

◆ gpuNavierStokes3DUpwindRusanov()

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

◆ gpuNavierStokes3DUpwindRoe()

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

◆ gpuNavierStokes3DPreStep()

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

◆ NavierStokes3DInitialize()

int NavierStokes3DInitialize ( void *  s,
void *  m 
)

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

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

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

where the list of keywords are:

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

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

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables

Definition at line 116 of file NavierStokes3DInitialize.c.

119 {
120  HyPar *solver = (HyPar*) s;
121  MPIVariables *mpi = (MPIVariables*) m;
122  NavierStokes3D *physics = (NavierStokes3D*) solver->physics;
123  int ferr = 0;
124 
125  static int count = 0;
126 
127  if (solver->nvars != _MODEL_NVARS_) {
128  fprintf(stderr,"Error in NavierStokes3DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
129  return(1);
130  }
131  if (solver->ndims != _MODEL_NDIMS_) {
132  fprintf(stderr,"Error in NavierStokes3DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
133  return(1);
134  }
135 
136  /* default values */
137  physics->gamma = 1.4;
138  physics->Pr = 0.72;
139  physics->Re = -1;
140  physics->Minf = 1.0;
141  physics->C1 = 1.458e-6;
142  physics->C2 = 110.4;
143  physics->grav_x = 0.0;
144  physics->grav_y = 0.0;
145  physics->grav_z = 0.0;
146  physics->rho0 = 1.0;
147  physics->p0 = 1.0;
148  physics->HB = 1;
149  physics->R = 1.0;
150  physics->N_bv = 0.0;
151  physics->T_ib_wall = -DBL_MAX;
152  physics->t_ib_ramp = -1.0;
153  physics->t_ib_width= 0.05;
154  physics->ib_T_tol = 5;
155  strcpy(physics->upw_choice,"roe");
156  strcpy(physics->ib_write_surface_data,"yes");
157  strcpy(physics->ib_wall_type,"adiabatic");
158  strcpy(physics->ib_ramp_type,"linear");
159 
160  /* reading physical model specific inputs - all processes */
161  if (!mpi->rank) {
162  FILE *in;
163  if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
164  in = fopen("physics.inp","r");
165  if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
166  else {
167  char word[_MAX_STRING_SIZE_];
168  ferr = fscanf(in,"%s",word);
169  if (ferr != 1) {
170  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
171  return 1;
172  }
173  if (!strcmp(word, "begin")){
174  while (strcmp(word, "end")){
175  ferr = fscanf(in,"%s",word);
176  if (ferr != 1) {
177  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
178  return 1;
179  }
180  if (!strcmp(word, "gamma")) {
181  ferr = fscanf(in,"%lf",&physics->gamma);
182  if (ferr != 1) {
183  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
184  return 1;
185  }
186  } else if (!strcmp(word,"upwinding")) {
187  ferr = fscanf(in,"%s",physics->upw_choice);
188  if (ferr != 1) {
189  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
190  return 1;
191  }
192  } else if (!strcmp(word,"Pr")) {
193  ferr = fscanf(in,"%lf",&physics->Pr);
194  if (ferr != 1) {
195  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
196  return 1;
197  }
198  } else if (!strcmp(word,"Re")) {
199  ferr = fscanf(in,"%lf",&physics->Re);
200  if (ferr != 1) {
201  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
202  return 1;
203  }
204  } else if (!strcmp(word,"Minf")) {
205  ferr = fscanf(in,"%lf",&physics->Minf);
206  if (ferr != 1) {
207  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
208  return 1;
209  }
210  } else if (!strcmp(word,"gravity")) {
211  ferr = fscanf(in,"%lf",&physics->grav_x);
212  if (ferr != 1) {
213  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
214  return 1;
215  }
216  ferr = fscanf(in,"%lf",&physics->grav_y);
217  if (ferr != 1) {
218  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
219  return 1;
220  }
221  ferr = fscanf(in,"%lf",&physics->grav_z);
222  if (ferr != 1) {
223  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
224  return 1;
225  }
226  } else if (!strcmp(word,"rho_ref")) {
227  ferr = fscanf(in,"%lf",&physics->rho0);
228  if (ferr != 1) {
229  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
230  return 1;
231  }
232  } else if (!strcmp(word,"p_ref")) {
233  ferr = fscanf(in,"%lf",&physics->p0);
234  if (ferr != 1) {
235  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
236  return 1;
237  }
238  } else if (!strcmp(word,"HB")) {
239  ferr = fscanf(in,"%d",&physics->HB);
240  if (ferr != 1) {
241  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
242  return 1;
243  }
244  if (physics->HB==3) {
245  ferr = fscanf(in,"%lf",&physics->N_bv);
246  if (ferr != 1) {
247  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
248  return 1;
249  }
250  }
251  } else if (!strcmp(word,"R")) {
252  ferr = fscanf(in,"%lf",&physics->R);
253  if (ferr != 1) {
254  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
255  return 1;
256  }
257  } else if (!strcmp(word,"ib_surface_data")) {
258  ferr = fscanf(in,"%s",physics->ib_write_surface_data);
259  if (ferr != 1) {
260  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
261  return 1;
262  }
263  } else if (!strcmp(word,"ib_wall_type")) {
264  ferr = fscanf(in,"%s",physics->ib_wall_type);
265  if (ferr != 1) {
266  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
267  return 1;
268  }
269  if (!strcmp(physics->ib_wall_type,_IB_ISOTHERMAL_)) {
270  ferr = fscanf(in,"%lf",&physics->T_ib_wall);
271  if (ferr != 1) {
272  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
273  return 1;
274  }
275  }
276  if (!solver->flag_ib) {
277  printf("Warning: in NavierStokes3DInitialize().\n");
278  printf("Warning: no immersed body present; specification of ib_wall_type unnecessary.\n");
279  }
280  } else if (!strcmp(word,"ib_ramp_time")) {
281  ferr = fscanf(in,"%lf",&physics->t_ib_ramp);
282  if (ferr != 1) {
283  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
284  return 1;
285  }
286  if (!solver->flag_ib) {
287  printf("Warning: in NavierStokes3DInitialize().\n");
288  printf("Warning: no immersed body present; specification of ib_ramp_time unnecessary.\n");
289  }
290  } else if (!strcmp(word,"ib_ramp_width")) {
291  ferr = fscanf(in,"%lf",&physics->t_ib_width);
292  if (ferr != 1) {
293  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
294  return 1;
295  }
296  if (!solver->flag_ib) {
297  printf("Warning: in NavierStokes3DInitialize().\n");
298  printf("Warning: no immersed body present; specification of ib_ramp_width unnecessary.\n");
299  }
300  } else if (!strcmp(word,"ib_ramp_type")) {
301  ferr = fscanf(in,"%s",physics->ib_ramp_type);
302  if (ferr != 1) {
303  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
304  return 1;
305  }
306  if (!solver->flag_ib) {
307  printf("Warning: in NavierStokes3DInitialize().\n");
308  printf("Warning: no immersed body present; specification of ib_ramp_type unnecessary.\n");
309  }
310  } else if (!strcmp(word,"ib_T_tolerance")) {
311  ferr = fscanf(in,"%lf",&physics->ib_T_tol);
312  if (ferr != 1) {
313  fprintf(stderr, "Read error while reading physics.inp in NavierStokes3DInitialize().\n");
314  return 1;
315  }
316  if (!solver->flag_ib) {
317  printf("Warning: in NavierStokes3DInitialize().\n");
318  printf("Warning: no immersed body present; specification of ib_T_tolerance unnecessary.\n");
319  }
320  } else if (strcmp(word,"end")) {
321  char useless[_MAX_STRING_SIZE_];
322  ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
323  printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
324  printf("recognized or extraneous. Ignoring.\n");
325  }
326  }
327  } else {
328  fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
329  return(1);
330  }
331  }
332  fclose(in);
333  }
334 
339  IERR MPIBroadcast_double (&physics->gamma ,1 ,0,&mpi->world); CHECKERR(ierr);
340  IERR MPIBroadcast_double (&physics->Pr ,1 ,0,&mpi->world); CHECKERR(ierr);
341  IERR MPIBroadcast_double (&physics->Re ,1 ,0,&mpi->world); CHECKERR(ierr);
342  IERR MPIBroadcast_double (&physics->Minf ,1 ,0,&mpi->world); CHECKERR(ierr);
343  IERR MPIBroadcast_double (&physics->grav_x ,1 ,0,&mpi->world); CHECKERR(ierr);
344  IERR MPIBroadcast_double (&physics->grav_y ,1 ,0,&mpi->world); CHECKERR(ierr);
345  IERR MPIBroadcast_double (&physics->grav_z ,1 ,0,&mpi->world); CHECKERR(ierr);
346  IERR MPIBroadcast_double (&physics->rho0 ,1 ,0,&mpi->world); CHECKERR(ierr);
347  IERR MPIBroadcast_double (&physics->p0 ,1 ,0,&mpi->world); CHECKERR(ierr);
348  IERR MPIBroadcast_double (&physics->R ,1 ,0,&mpi->world); CHECKERR(ierr);
349  IERR MPIBroadcast_double (&physics->N_bv ,1 ,0,&mpi->world); CHECKERR(ierr);
350  IERR MPIBroadcast_double (&physics->T_ib_wall ,1 ,0,&mpi->world); CHECKERR(ierr);
351  IERR MPIBroadcast_double (&physics->t_ib_ramp ,1 ,0,&mpi->world); CHECKERR(ierr);
352  IERR MPIBroadcast_double (&physics->t_ib_width ,1 ,0,&mpi->world); CHECKERR(ierr);
353  IERR MPIBroadcast_double (&physics->ib_T_tol ,1 ,0,&mpi->world); CHECKERR(ierr);
354  IERR MPIBroadcast_integer (&physics->HB ,1 ,0,&mpi->world); CHECKERR(ierr);
355 
356  /* if file output is disabled in HyPar, respect that */
357  if (!strcmp(solver->op_file_format,"none")) {
358  if (!strcmp(physics->ib_write_surface_data,"yes")) {
359  if (!mpi->rank) {
360  printf("Warning from NavierStokes3DInitialize(): solver->op_file_format is set to \"none\", thus ");
361  printf("setting physics->ib_write_surface_data to \"no\" (no solution files will be written).\n");
362  }
363  }
364  strcpy(physics->ib_write_surface_data,"no");
365  }
366 
367  /* Scaling Re by M_inf */
368  physics->Re /= physics->Minf;
369 
370  /* check that a well-balanced upwinding scheme is being used for cases with gravity */
371  if ( ((physics->grav_x != 0.0) || (physics->grav_y != 0.0) || (physics->grav_z != 0.0) )
372  && (strcmp(physics->upw_choice,_RUSANOV_)) ) {
373  if (!mpi->rank) {
374  fprintf(stderr,"Error in NavierStokes3DInitialize: %s upwinding is needed for flows ",_RUSANOV_);
375  fprintf(stderr,"with gravitational forces.\n");
376  }
377  return(1);
378  }
379  /* check that solver has the correct choice of diffusion formulation, if viscous flow */
380  if (strcmp(solver->spatial_type_par,_NC_2STAGE_) && (physics->Re > 0)) {
381  if (!mpi->rank) {
382  fprintf(stderr,"Error in NavierStokes3DInitialize(): Parabolic term spatial discretization must be \"%s\"\n",_NC_2STAGE_);
383  }
384  return(1);
385  }
386 
387  /* initializing physical model-specific functions */
388 #if defined(HAVE_CUDA)
389  if (solver->use_gpu) {
394  } else {
395 #endif
396  solver->PreStep = NavierStokes3DPreStep;
398  solver->FFunction = NavierStokes3DFlux;
404 #if defined(HAVE_CUDA)
405  }
406 #endif
407 
408  if (solver->flag_ib) {
409  if (!strcmp(physics->ib_wall_type,_IB_ADIABATIC_)) {
411  } else if (!strcmp(physics->ib_wall_type,_IB_ISOTHERMAL_)) {
413  } else {
414  fprintf(stderr, "Error in NavierStokes3DInitialize()\n");
415  fprintf(stderr, " invalid value for IB wall type (%s).\n",
416  physics->ib_wall_type );
417  }
418  if (!strcmp(physics->ib_write_surface_data,"yes")) {
420  }
421  }
422 
423 #if defined(HAVE_CUDA)
424  if (solver->use_gpu) {
425 
426  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
427  if (!mpi->rank) {
428  fprintf(stderr,"Error in NavierStokes3DInitialize(): Not available on GPU!");
429  }
430  return 1;
431  } else {
433  if (!strcmp(physics->upw_choice,_ROE_ )) {
435  } else if (!strcmp(physics->upw_choice,_RUSANOV_)) {
437  } else {
438  if (!mpi->rank) {
439  fprintf(stderr,"Error in NavierStokes3DInitialize(): %s is not a valid upwinding scheme on GPU. ",
440  physics->upw_choice);
441  fprintf(stderr,"Choices are %s and %s.\n",_ROE_,_RUSANOV_);
442  }
443  return(1);
444  }
445  }
446 
447  } else {
448 #endif
449 
450  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
454  if (!strcmp(physics->upw_choice,_ROE_)) {
458  } else if (!strcmp(physics->upw_choice,_RUSANOV_)) {
462  } else {
463  if (!mpi->rank) {
464  fprintf(stderr,"Error in NavierStokes3DInitialize(): %s is not a valid upwinding scheme ",
465  physics->upw_choice);
466  fprintf(stderr,"for use with split hyperbolic flux form. Use %s or %s.\n",
467  _ROE_,_RUSANOV_);
468  }
469  return(1);
470  }
471  } else {
473  if (!strcmp(physics->upw_choice,_ROE_ )) solver->Upwind = NavierStokes3DUpwindRoe;
474  else if (!strcmp(physics->upw_choice,_RF_ )) solver->Upwind = NavierStokes3DUpwindRF;
475  else if (!strcmp(physics->upw_choice,_LLF_ )) solver->Upwind = NavierStokes3DUpwindLLF;
476  else if (!strcmp(physics->upw_choice,_RUSANOV_)) solver->Upwind = NavierStokes3DUpwindRusanov;
477  else {
478  if (!mpi->rank) {
479  fprintf(stderr,"Error in NavierStokes3DInitialize(): %s is not a valid upwinding scheme. ",
480  physics->upw_choice);
481  fprintf(stderr,"Choices are %s, %s, %s, and %s.\n",_ROE_,_RF_,_LLF_,_RUSANOV_);
482  }
483  return(1);
484  }
485  }
486 
487 #if defined(HAVE_CUDA)
488  }
489 #endif
490 
491  /* set the value of gamma in all the boundary objects */
492  int n;
493  DomainBoundary *boundary = (DomainBoundary*) solver->boundary;
494  for (n = 0; n < solver->nBoundaryZones; n++) boundary[n].gamma = physics->gamma;
495 
496  /* finally, hijack the main solver's dissipation function pointer
497  * to this model's own function, since it's difficult to express
498  * the dissipation terms in the general form */
499 #if defined(HAVE_CUDA)
500  if (solver->use_gpu) {
502  } else {
503 #endif
505 #if defined(HAVE_CUDA)
506  }
507 #endif
508 
509  /* allocate array to hold the gravity field */
510  int *dim = solver->dim_local;
511  int ghosts = solver->ghosts;
512  int d, size = 1; for (d=0; d<_MODEL_NDIMS_; d++) size *= (dim[d] + 2*ghosts);
513  physics->grav_field_f = (double*) calloc (size, sizeof(double));
514  physics->grav_field_g = (double*) calloc (size, sizeof(double));
515  /* allocate arrays to hold the fast Jacobian for split form of the hyperbolic flux */
516  physics->fast_jac = (double*) calloc (_MODEL_NDIMS_*size*_MODEL_NVARS_*_MODEL_NVARS_,sizeof(double));
517  physics->solution = (double*) calloc (size*_MODEL_NVARS_,sizeof(double));
518  /* initialize the gravity fields */
519  IERR NavierStokes3DGravityField(solver,mpi); CHECKERR(ierr);
520 
521 #if defined(HAVE_CUDA)
522  if (solver->use_gpu) gpuNavierStokes3DInitialize(s,m);
523 #endif
524 
525  count++;
526  return(0);
527 }
int NavierStokes3DGravityField(void *, void *)
int gpuNavierStokes3DUpwindRusanov(double *, double *, double *, double *, double *, double *, int, void *, double)
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
int NavierStokes3DRightEigenvectors(double *, double *, void *, int)
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
int gpuNavierStokes3DParabolicFunction(double *, double *, void *, void *, double)
#define CHECKERR(ierr)
Definition: basic.h:18
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
void * boundary
Definition: hypar.h:159
int gpuNavierStokes3DInitialize(void *, void *)
Structure containing variables and parameters specific to the 3D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 3D Navier-Stokes equations.
int NavierStokes3DUpwinddFRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23
int NavierStokes3DRoeAverage(double *, double *, double *, void *)
int gpuNavierStokes3DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
int NavierStokes3DUpwinddFRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
int(* AveragingFunction)(double *, double *, double *, void *)
Definition: hypar.h:354
double * grav_field_g
int NavierStokes3DStiffFlux(double *, double *, int, void *, double)
int flag_ib
Definition: hypar.h:441
char ib_ramp_type[_MAX_STRING_SIZE_]
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
int(* SFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:317
int(* IBFunction)(void *, void *, double *, double)
Definition: hypar.h:446
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
Structure containing the variables and function pointers defining a boundary.
int NavierStokes3DUpwindFdFRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
int ndims
Definition: hypar.h:26
int NavierStokes3DJacobian(double *, double *, void *, int, int, int)
int MPIBroadcast_character(char *, int, int, void *)
Definition: MPIBroadcast.c:37
int NavierStokes3DUpwindFdFRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
int gpuNavierStokes3DPreStep(double *, void *, void *, double)
int NavierStokes3DIBForces(void *, void *, double)
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int MPIBroadcast_double(double *, int, int, void *)
Definition: MPIBroadcast.c:9
double(* ComputeCFL)(void *, void *, double, double)
Definition: hypar.h:269
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int NavierStokes3DLeftEigenvectors(double *, double *, void *, int)
#define _MAX_STRING_SIZE_
Definition: basic.h:14
double NavierStokes3DComputeCFL(void *, void *, double, double)
int gpuNavierStokes3DFlux(double *, double *, int, void *, double)
#define _ROE_
Definition: euler1d.h:62
int NavierStokes3DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
int NavierStokes3DNonStiffFlux(double *, double *, int, void *, double)
int NavierStokes3DStiffJacobian(double *, double *, void *, int, int, int)
int NavierStokes3DPreStep(double *, void *, void *, double)
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int NavierStokes3DFlux(double *, double *, int, void *, double)
#define _RUSANOV_
Definition: euler1d.h:70
int NavierStokes3DParabolicFunction(double *, double *, void *, void *, double)
int NavierStokes3DSource(double *, double *, void *, void *, double)
MPI_Comm world
char upw_choice[_MAX_STRING_SIZE_]
int(* PreStep)(double *, void *, void *, double)
Definition: hypar.h:339
int NavierStokes3DUpwindRusanovModified(double *, double *, double *, double *, double *, double *, int, void *, double)
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
char spatial_type_par[_MAX_STRING_SIZE_]
Definition: hypar.h:96
int NavierStokes3DIBIsothermal(void *, void *, double *, double)
int * dim_local
Definition: hypar.h:37
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
char ib_write_surface_data[_MAX_STRING_SIZE_]
int NavierStokes3DUpwindRF(double *, double *, double *, double *, double *, double *, int, void *, double)
void * physics
Definition: hypar.h:266
int(* JFunction)(double *, double *, void *, int, int, int)
Definition: hypar.h:326
int NavierStokes3DModifiedSolution(double *, double *, int, void *, void *, double)
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
int use_gpu
Definition: hypar.h:449
#define _MODEL_NVARS_
Definition: euler1d.h:58
char ib_wall_type[_MAX_STRING_SIZE_]
int nBoundaryZones
Definition: hypar.h:157
int NavierStokes3DUpwindRusanov(double *, double *, double *, double *, double *, double *, int, void *, double)
#define _IB_ADIABATIC_
char op_file_format[_MAX_STRING_SIZE_]
Definition: hypar.h:186
int NavierStokes3DIBAdiabatic(void *, void *, double *, double)
int(* GetRightEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:359
int gpuNavierStokes3DModifiedSolution(double *, double *, int, void *, void *, double)
#define _RF_
Definition: euler1d.h:64
int(* GetLeftEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:357
int gpuNavierStokes3DSource(double *, double *, void *, void *, double)
int NavierStokes3DUpwindLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
#define _IB_ISOTHERMAL_
double * grav_field_f
int(* UFunction)(double *, double *, int, void *, void *, double)
Definition: hypar.h:321
#define _NC_2STAGE_
Definition: hypar.h:477
#define _LLF_
Definition: euler1d.h:66