HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
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

◆ NavierStokes2DGravityField()

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 }
HyPar::x
double * x
Definition: hypar.h:107
NavierStokes2D::grav_field_g
double * grav_field_g
Definition: navierstokes2d.h:388
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
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
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
NavierStokes2D::grav_field_f
double * grav_field_f
Definition: navierstokes2d.h:388
_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
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
raiseto
#define raiseto(x, a)
Definition: math_ops.h:37
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319