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

Compute the gravitational source term for the 3D Navier Stokes system. More...

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <physicalmodels/navierstokes3d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions

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

Detailed Description

Compute the gravitational source term for the 3D Navier Stokes system.

Author
Debojyoti Ghosh

Definition in file NavierStokes3DSource.c.

Function Documentation

◆ NavierStokes3DSourceFunction()

int NavierStokes3DSourceFunction ( double *  f,
double *  u,
double *  x,
void *  s,
void *  m,
double  t,
int  dir 
)

Compute the source function in the well-balanced treatment of the source terms. The source function is:

\begin{equation} dir = x \rightarrow \left[\begin{array}{c}0 \\ \varphi\left(x,y,z\right) \\ 0 \\ 0 \\ \varphi\left(x,y,z\right) \end{array}\right], \ dir = y \rightarrow \left[\begin{array}{c}0 \\ 0 \\ \varphi\left(x,y,z\right) \\ 0 \\ \varphi\left(x,y,z\right) \end{array}\right], \ dir = z \rightarrow \left[\begin{array}{c}0 \\ 0 \\ 0 \\ \varphi\left(x,y,z\right) \\ \varphi\left(x,y,z\right) \end{array}\right] \end{equation}

where \(\varphi\) (NavierStokes3D::grav_field_g) is computed in NavierStokes3D::GravityField().

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, http://dx.doi.org/10.2514/1.J054580.
Parameters
fArray to hold the computed source function
uSolution vector array
xArray of spatial coordinates (grid)
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time
dirSpatial dimension (x, y, or z)

Definition at line 124 of file NavierStokes3DSource.c.

133 {
134  HyPar *solver = (HyPar* ) s;
135  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
136 
137  int ghosts = solver->ghosts;
138  int *dim = solver->dim_local;
139  static int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
140 
141  /* set bounds for array index to include ghost points */
142  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_);
143  int i; for (i=0; i<_MODEL_NDIMS_; i++) bounds[i] += 2*ghosts;
144 
145  /* set offset such that index is compatible with ghost point arrangement */
146  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
147 
148  int done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
149  while (!done) {
150  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
151  (f+_MODEL_NVARS_*p)[0] = 0.0;
152  (f+_MODEL_NVARS_*p)[1] = param->grav_field_g[p] * (dir == _XDIR_);
153  (f+_MODEL_NVARS_*p)[2] = param->grav_field_g[p] * (dir == _YDIR_);
154  (f+_MODEL_NVARS_*p)[3] = param->grav_field_g[p] * (dir == _ZDIR_);
155  (f+_MODEL_NVARS_*p)[4] = param->grav_field_g[p];
156  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
157  }
158 
159  return(0);
160 }
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
NavierStokes3D
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.
Definition: navierstokes3d.h:482
NavierStokes3D::grav_field_g
double * grav_field_g
Definition: navierstokes3d.h:497
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ NavierStokes3DSourceUpwind()

int NavierStokes3DSourceUpwind ( double *  fI,
double *  fL,
double *  fR,
double *  u,
int  dir,
void *  s,
double  t 
)

Compute the "upwind" source function value at the interface: the upwinding is just the arithmetic average of the left and right biased interpolated values of the source function at the interface.

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, Submitted
Parameters
fIArray to hold the computed "upwind" interface source function
fLInterface source function value computed using left-biased interpolation
fRInterface source function value computed using right-biased interpolation
uSolution vector array
dirSpatial dimension (x,y, or z)
sSolver object of type HyPar
tCurrent simulation time

Definition at line 174 of file NavierStokes3DSource.c.

183 {
184  HyPar *solver = (HyPar*) s;
185  int done,k, *dim = solver->dim_local;
186  _DECLARE_IERR_;
187 
188 
189  static int index_outer[_MODEL_NDIMS_], index_inter[_MODEL_NDIMS_],
190  bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
191  _ArrayCopy1D_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
192  _ArrayCopy1D_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
193 
194  done = 0; _ArraySetValue_(index_outer,_MODEL_NDIMS_,0);
195  while (!done) {
196  _ArrayCopy1D_(index_outer,index_inter,_MODEL_NDIMS_);
197  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
198  int p; _ArrayIndex1D_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
199  for (k = 0; k < _MODEL_NVARS_; k++)
200  (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
201  }
202  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
203  }
204 
205  return(0);
206 }
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ 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 }
IERR
#define IERR
Definition: basic.h:16
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
NavierStokes3D
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.
Definition: navierstokes3d.h:482
HyPar::x
double * x
Definition: hypar.h:107
HyPar::fluxI
double * fluxI
Definition: hypar.h:136
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
NavierStokes3DSourceUpwind
int NavierStokes3DSourceUpwind(double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes3DSource.c:174
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
HyPar::InterpolateInterfacesHyp
int(* InterpolateInterfacesHyp)(double *, double *, double *, double *, int, int, void *, void *, int)
Definition: hypar.h:224
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_NavierStokes3D_stride_
static const int _NavierStokes3D_stride_
Definition: navierstokes3d.h:559
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::fR
double * fR
Definition: hypar.h:139
HyPar::fL
double * fL
Definition: hypar.h:139
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
HyPar::ghosts
int ghosts
Definition: hypar.h:52
NavierStokes3DSourceFunction
int NavierStokes3DSourceFunction(double *, double *, double *, void *, void *, double, int)
Definition: NavierStokes3DSource.c:124
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::fluxC
double * fluxC
Definition: hypar.h:128
NavierStokes3D::grav_field_f
double * grav_field_f
Definition: navierstokes3d.h:497
HyPar::dxinv
double * dxinv
Definition: hypar.h:110
_NavierStokes3DGetFlowVar_
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
Definition: navierstokes3d.h:98