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
NavierStokes2DGravityField.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/navierstokes2d.h>
#include <hypar.h>
#include <mpivars.h>

Go to the source code of this file.

Functions

int NavierStokes2DGravityField (void *s, void *m)
 

Detailed Description

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

Author
Debojyoti Ghosh

Definition in file NavierStokes2DGravityField.c.

Function Documentation

int NavierStokes2DGravityField ( void *  s,
void *  m 
)

This function computes the pressure and density variation functions for the hydrostatic balance of the type specified by the user (NavierStokes2D: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 (NavierStokes2D::rho0) and pressure (NavierStokes2D::p0). This function computes \(\varrho\) (NavierStokes2D::grav_field_f) and \(\varphi\) (NavierStokes2D::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 NavierStokes2DGravityField.c.

37 {
38  HyPar *solver = (HyPar*) s;
39  MPIVariables *mpi = (MPIVariables*) m;
40  NavierStokes2D *param = (NavierStokes2D*) 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 Nbv = param->N_bv;
65 
66  /* set the value of the gravity field */
67  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
68  if (param->HB == 1) {
69  while (!done) {
70  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
71  double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
72  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
73  f[p] = exp( (gx*xcoord+gy*ycoord)/RT);
74  g[p] = exp(-(gx*xcoord+gy*ycoord)/RT);
75  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
76  }
77  } else if (param->HB == 2) {
78  while (!done) {
79  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
80  double xcoord; _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->x,xcoord);
81  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
82  f[p] = raiseto((1.0-(gx*xcoord+gy*ycoord)/(Cp*T0)), (-1.0 /(gamma-1.0)));
83  g[p] = raiseto((1.0-(gx*xcoord+gy*ycoord)/(Cp*T0)), (gamma/(gamma-1.0)));
84  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
85  }
86  } else if (param->HB == 3) {
87  while (!done) {
88  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
89  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
90  if (gx != 0) {
91  if (!mpi->rank) {
92  fprintf(stderr,"Error in NavierStokes2DGravityField(): HB = 3 is implemented only for ");
93  fprintf(stderr,"gravity force along the y-coordinate.\n");
94  }
95  return(1);
96  }
97  double Pexner = 1 + ((gy*gy)/(Cp*T0*Nbv*Nbv)) * (exp(-(Nbv*Nbv/gy)*ycoord)-1.0);
98  f[p] = raiseto(Pexner, (-1.0 /(gamma-1.0))) * exp(Nbv*Nbv*ycoord/gy);
99  g[p] = raiseto(Pexner, (gamma/(gamma-1.0)));
100  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
101  }
102  } else {
103  while (!done) {
104  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
105  f[p] = g[p] = 1.0;
106  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
107  }
108  }
109 
110  /* A sensible simulation will not specify peridic boundary conditions
111  * along a direction in which gravity acts.
112  * Gravity will be zero along the dimension periodic BCs are specified,
113  * so the value of the gravity potential will be the same along those
114  * grid lines.
115  * Thus, for both these cases, extrapolate the gravity field at the
116  * boundaries */
117  int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_];
118  for (d = 0; d < _MODEL_NDIMS_; d++) {
119  /* left boundary */
120  if (!mpi->ip[d]) {
121  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
122  _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
123  done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
124  while (!done) {
125  _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = ghosts-1-indexb[d];
126  int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
127  int p2; _ArrayIndex1D_ (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
128  f[p1] = f[p2];
129  g[p1] = g[p2];
130  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
131  }
132  }
133  /* right boundary */
134  if (mpi->ip[d] == mpi->iproc[d]-1) {
135  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
136  _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
137  done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
138  while (!done) {
139  _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = dim[d]-1-indexb[d];
140  int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
141  int p2; _ArrayIndex1D_ (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
142  f[p1] = f[p2];
143  g[p1] = g[p2];
144  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
145  }
146  }
147  }
148 
149  return(0);
150 }
_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
NavierStokes2D::grav_field_f
double * grav_field_f
Definition: navierstokes2d.h:388
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
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