HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
NavierStokes2DSource.c File Reference

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

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

Go to the source code of this file.

Functions

static int NavierStokes2DSourceFunction (double *, double *, double *, void *, void *, double, int)
 
static int NavierStokes2DSourceUpwind (double *, double *, double *, double *, int, void *, double)
 
int NavierStokes2DSource (double *source, double *u, void *s, void *m, double t)
 

Detailed Description

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

Author
Debojyoti Ghosh

Definition in file NavierStokes2DSource.c.

Function Documentation

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

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\right) \\ 0 \\ \varphi\left(x,y\right) \end{array}\right], \ dir = y \rightarrow \left[\begin{array}{c}0 \\ 0 \\ \varphi\left(x,y\right) \\ \varphi\left(x,y\right) \end{array}\right] \end{equation}

where \(\varphi\) (NavierStokes2D::grav_field_g) is computed in NavierStokes2D::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, Submitted
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 or y)

Definition at line 151 of file NavierStokes2DSource.c.

160 {
161  HyPar *solver = (HyPar* ) s;
162  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
163 
164  int ghosts = solver->ghosts;
165  int *dim = solver->dim_local;
166  int ndims = solver->ndims;
167  int index[ndims], bounds[ndims], offset[ndims];
168 
169  /* set bounds for array index to include ghost points */
170  _ArrayCopy1D_(dim,bounds,ndims);
171  int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
172 
173  /* set offset such that index is compatible with ghost point arrangement */
174  _ArraySetValue_(offset,ndims,-ghosts);
175 
176  int done = 0; _ArraySetValue_(index,ndims,0);
177  while (!done) {
178  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
179  (f+_MODEL_NVARS_*p)[0] = 0.0;
180  (f+_MODEL_NVARS_*p)[1] = param->grav_field_g[p] * (dir == _XDIR_);
181  (f+_MODEL_NVARS_*p)[2] = param->grav_field_g[p] * (dir == _YDIR_);
182  (f+_MODEL_NVARS_*p)[3] = param->grav_field_g[p];
183  _ArrayIncrementIndex_(ndims,bounds,index,done);
184  }
185 
186  return(0);
187 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
NavierStokes2D::grav_field_g
double * grav_field_g
Definition: navierstokes2d.h:388
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
int NavierStokes2DSourceUpwind ( double *  fI,
double *  fL,
double *  fR,
double *  u,
int  dir,
void *  s,
double  t 
)
static

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 or y)
sSolver object of type HyPar
tCurrent simulation time

Definition at line 200 of file NavierStokes2DSource.c.

209 {
210  HyPar *solver = (HyPar*) s;
211  int done,k;
212  _DECLARE_IERR_;
213 
214  int ndims = solver->ndims;
215  int *dim = solver->dim_local;
216 
217  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
218  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
219  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
220 
221  done = 0; _ArraySetValue_(index_outer,ndims,0);
222  while (!done) {
223  _ArrayCopy1D_(index_outer,index_inter,ndims);
224  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
225  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
226  for (k = 0; k < _MODEL_NVARS_; k++)
227  (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
228  }
229  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
230  }
231 
232  return(0);
233 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
int NavierStokes2DSource ( double *  source,
double *  u,
void *  s,
void *  m,
double  t 
)

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

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

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

References:

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

44 {
45  HyPar *solver = (HyPar* ) s;
46  MPIVariables *mpi = (MPIVariables*) m;
47  NavierStokes2D *param = (NavierStokes2D*) solver->physics;
48 
49  if ((param->grav_x == 0.0) && (param->grav_y == 0.0))
50  return(0); /* no gravitational forces */
51 
52  int v, done, p, p1, p2;
53  double *SourceI = solver->fluxI; /* interace source term */
54  double *SourceC = solver->fluxC; /* cell-centered source term */
55  double *SourceL = solver->fL;
56  double *SourceR = solver->fR;
57 
58  int ndims = solver->ndims;
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  int index[ndims],index1[ndims],index2[ndims],dim_interface[ndims];
65 
66  /* Along X-direction */
67  if (param->grav_x != 0.0) {
68  /* set interface dimensions */
69  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;
70  /* calculate the split source function exp(-phi/RT) */
71  IERR NavierStokes2DSourceFunction(SourceC,u,x,solver,mpi,t,_XDIR_); CHECKERR(ierr);
72  /* calculate the left and right interface source terms */
73  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
74  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
75  /* calculate the final interface source term */
76  IERR NavierStokes2DSourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
77  /* calculate the final cell-centered source term */
78  done = 0; _ArraySetValue_(index,ndims,0);
79  while (!done) {
80  _ArrayCopy1D_(index,index1,ndims);
81  _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
82  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
83  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
84  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
85 
86  double dx_inverse; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
87  double rho, uvel, vvel, e, P;
88  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,e,P,param->gamma);
89  double term[_MODEL_NVARS_] = {0.0, rho*RT, 0.0, rho*RT*uvel};
90  for (v=0; v<_MODEL_NVARS_; v++) {
91  source[_MODEL_NVARS_*p+v] += ( (term[v]*param->grav_field_f[p])
92  * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse );
93  }
94  uvel = P; /* useless statement to avoid compiler warnings */
95  _ArrayIncrementIndex_(ndims,dim,index,done);
96  }
97  }
98 
99  /* Along Y-direction */
100  if (param->grav_y != 0.0) {
101  /* set interface dimensions */
102  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_YDIR_]++;
103  /* calculate the split source function exp(-phi/RT) */
104  IERR NavierStokes2DSourceFunction(SourceC,u,x,solver,mpi,t,_YDIR_); CHECKERR(ierr);
105  /* calculate the left and right interface source terms */
106  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
107  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
108  /* calculate the final interface source term */
109  IERR NavierStokes2DSourceUpwind(SourceI,SourceL,SourceR,u,_YDIR_,solver,t);
110 
111  /* calculate the final cell-centered source term */
112  done = 0; _ArraySetValue_(index,ndims,0);
113  while (!done) {
114  _ArrayCopy1D_(index,index1,ndims);
115  _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
116  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
117  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
118  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
119 
120  double dy_inverse; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
121  double rho, uvel, vvel, e, P;
122  _NavierStokes2DGetFlowVar_((u+_MODEL_NVARS_*p),rho,uvel,vvel,e,P,param->gamma);
123  double term[_MODEL_NVARS_] = {0.0, 0.0, rho*RT, rho*RT*vvel};
124  for (v=0; v<_MODEL_NVARS_; v++) {
125  source[_MODEL_NVARS_*p+v] += ( (term[v]*param->grav_field_f[p])
126  * (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dy_inverse );
127  }
128  uvel = P; /* useless statement to avoid compiler warnings */
129  _ArrayIncrementIndex_(ndims,dim,index,done);
130  }
131  }
132 
133  return(0);
134 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::fL
double * fL
Definition: hypar.h:139
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
HyPar::fluxI
double * fluxI
Definition: hypar.h:136
NavierStokes2D::grav_field_f
double * grav_field_f
Definition: navierstokes2d.h:388
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::fluxC
double * fluxC
Definition: hypar.h:128
HyPar::InterpolateInterfacesHyp
int(* InterpolateInterfacesHyp)(double *, double *, double *, double *, int, int, void *, void *, int)
Definition: hypar.h:224
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
_NavierStokes2DGetFlowVar_
#define _NavierStokes2DGetFlowVar_(u, rho, vx, vy, e, P, gamma)
Definition: navierstokes2d.h:81
HyPar::dxinv
double * dxinv
Definition: hypar.h:110
NavierStokes2DSourceFunction
static int NavierStokes2DSourceFunction(double *, double *, double *, void *, void *, double, int)
Definition: NavierStokes2DSource.c:151
NavierStokes2D
Structure containing variables and parameters specific to the 2D Navier Stokes equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D Navier-Stokes equations.
Definition: navierstokes2d.h:374
HyPar::ndims
int ndims
Definition: hypar.h:26
IERR
#define IERR
Definition: basic.h:16
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::x
double * x
Definition: hypar.h:107
NavierStokes2DSourceUpwind
static int NavierStokes2DSourceUpwind(double *, double *, double *, double *, int, void *, double)
Definition: NavierStokes2DSource.c:200
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
HyPar::fR
double * fR
Definition: hypar.h:139