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
NavierStokes3DGravityField.c File Reference

Function to compute the hydrostatically balanced pressure and density variation functions. More...

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

Go to the source code of this file.

Functions

int NavierStokes3DGravityField (void *s, void *m)
 

Detailed Description

Function to compute the hydrostatically balanced pressure and density variation functions.

Author
Debojyoti Ghosh

Definition in file NavierStokes3DGravityField.c.

Function Documentation

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 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
raiseto
#define raiseto(x, a)
Definition: math_ops.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
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
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ZDIR_
#define _ZDIR_
Definition: navierstokes3d.h:76
NavierStokes3D::grav_field_f
double * grav_field_f
Definition: navierstokes3d.h:497
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
NavierStokes3D::grav_field_g
double * grav_field_g
Definition: navierstokes3d.h:497
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
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::x
double * x
Definition: hypar.h:107
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75