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

Subsonic "ambivalent" boundary conditions (specific to Euler/Navier-Stokes systems) More...

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <boundaryconditions.h>
#include <physicalmodels/euler2d.h>
#include <physicalmodels/navierstokes3d.h>

Go to the source code of this file.

Functions

int BCSubsonicAmbivalentU (void *b, void *m, int ndims, int nvars, int *size, int ghosts, double *phi, double waqt)
 

Detailed Description

Subsonic "ambivalent" boundary conditions (specific to Euler/Navier-Stokes systems)

Author
Debojyoti Ghosh

Definition in file BCSubsonicAmbivalent.c.

Function Documentation

int BCSubsonicAmbivalentU ( void *  b,
void *  m,
int  ndims,
int  nvars,
int *  size,
int  ghosts,
double *  phi,
double  waqt 
)

Applies the subsonic "ambivalent" boundary condition: The velocity at the boundary face is extrapolated from the interior and its dot product with the boundary normal (pointing into the domain) is computed.

  • If it is positive, subsonic inflow boundary conditions are applied: the density and velocity at the physical boundary ghost points are specified, while the pressure is extrapolated from the interior of the domain.
  • If it is negative, subsonic outflow boundary conditions are applied: the pressure at the physical boundary ghost points is specified, while the density and velocity are extrapolated from the interior.

This boundary condition is specific to two and three dimension Euler and Navier-Stokes systems (Euler2D, NavierStokes2D, NavierStokes3D).

Parameters
bBoundary object of type DomainBoundary
mMPI object of type MPIVariables
ndimsNumber of spatial dimensions
nvarsNumber of variables/DoFs per grid point
sizeInteger array with the number of grid points in each spatial dimension
ghostsNumber of ghost points
phiThe solution array on which to apply the boundary condition
waqtCurrent solution time

Definition at line 28 of file BCSubsonicAmbivalent.c.

38 {
39  DomainBoundary *boundary = (DomainBoundary*) b;
40 
41  int dim = boundary->dim;
42  int face = boundary->face;
43 
44  if (ndims == 2) {
45 
46  /* create a fake physics object */
47  Euler2D physics;
48  double gamma;
49  gamma = physics.gamma = boundary->gamma;
50  double inv_gamma_m1 = 1.0/(gamma-1.0);
51 
52  /* boundary normal (pointing into the domain) */
53  double nx, ny;
54  if (dim == 0) {
55  nx = 1.0;
56  ny = 0.0;
57  } else if (dim == 1) {
58  nx = 0.0;
59  ny = 1.0;
60  }
61  nx *= (double) face;
62  ny *= (double) face;
63 
64  if (boundary->on_this_proc) {
65  int bounds[ndims], indexb[ndims], indexi[ndims], indexj[ndims];
66  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
67  _ArraySetValue_(indexb,ndims,0);
68  int done = 0;
69  while (!done) {
70  int p1, p2;
71  double rho, uvel, vvel, energy, pressure;
72 
73  /* compute boundary face velocity - 2nd order */
74  _ArrayCopy1D_(indexb,indexi,ndims);
75  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
76  _ArrayCopy1D_(indexi,indexj,ndims);
77  if (face == 1) {
78  indexi[dim] = 0;
79  indexj[dim] = indexi[dim] + 1;
80  } else if (face == -1) {
81  indexi[dim] = size[dim]-1;
82  indexj[dim] = indexi[dim] - 1;
83  }
84  _ArrayIndex1D_(ndims,size,indexi,ghosts,p1);
85  _ArrayIndex1D_(ndims,size,indexj,ghosts,p2);
86  double uvel1, uvel2, uvelb,
87  vvel1, vvel2, vvelb;
88  _Euler2DGetFlowVar_((phi+nvars*p1),rho,uvel1,vvel1,energy,pressure,(&physics));
89  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel2,vvel2,energy,pressure,(&physics));
90  uvelb = 1.5*uvel1 - 0.5*uvel2;
91  vvelb = 1.5*vvel1 - 0.5*vvel2;
92  double vel_normal = uvelb*nx + vvelb*ny;
93 
94  _ArrayCopy1D_(indexb,indexi,ndims);
95  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
96  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
97  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
98  else return(1);
99  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
100  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
101 
102  /* flow variables in the interior */
103  _Euler2DGetFlowVar_((phi+nvars*p2),rho,uvel,vvel,energy,pressure,(&physics));
104 
105  /* set the ghost point values */
106  double rho_gpt, uvel_gpt, vvel_gpt, energy_gpt, pressure_gpt;
107  if (vel_normal > 0) {
108  /* inflow */
109  rho_gpt = boundary->FlowDensity;
110  pressure_gpt = pressure;
111  uvel_gpt = boundary->FlowVelocity[0];
112  vvel_gpt = boundary->FlowVelocity[1];
113  } else {
114  /* outflow */
115  rho_gpt = rho;
116  pressure_gpt = boundary->FlowPressure;
117  uvel_gpt = uvel;
118  vvel_gpt = vvel;
119  }
120  energy_gpt = inv_gamma_m1*pressure_gpt
121  + 0.5 * rho_gpt * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt);
122 
123  phi[nvars*p1+0] = rho_gpt;
124  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
125  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
126  phi[nvars*p1+3] = energy_gpt;
127 
128  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
129  }
130  }
131 
132  } else if (ndims == 3) {
133 
134  /* create a fake physics object */
135  double gamma;
136  gamma = boundary->gamma;
137  double inv_gamma_m1 = 1.0/(gamma-1.0);
138 
139  /* boundary normal (pointing into the domain) */
140  double nx, ny, nz;
141  if (dim == 0) {
142  nx = 1.0;
143  ny = 0.0;
144  nz = 0.0;
145  } else if (dim == 1) {
146  nx = 0.0;
147  ny = 1.0;
148  nz = 0.0;
149  } else if (dim == 2) {
150  nx = 0.0;
151  ny = 0.0;
152  nz = 1.0;
153  }
154  nx *= (double) face;
155  ny *= (double) face;
156  nz *= (double) face;
157 
158  if (boundary->on_this_proc) {
159  int bounds[ndims], indexb[ndims], indexi[ndims], indexj[ndims];
160  _ArraySubtract1D_(bounds,boundary->ie,boundary->is,ndims);
161  _ArraySetValue_(indexb,ndims,0);
162  int done = 0;
163  while (!done) {
164  int p1, p2;
165  double rho, uvel, vvel, wvel, energy, pressure;
166 
167  /* compute boundary face velocity - 2nd order */
168  _ArrayCopy1D_(indexb,indexi,ndims);
169  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
170  _ArrayCopy1D_(indexi,indexj,ndims);
171  if (face == 1) {
172  indexi[dim] = 0;
173  indexj[dim] = indexi[dim] + 1;
174  } else if (face == -1) {
175  indexi[dim] = size[dim]-1;
176  indexj[dim] = indexi[dim] - 1;
177  }
178  _ArrayIndex1D_(ndims,size,indexi,ghosts,p1);
179  _ArrayIndex1D_(ndims,size,indexj,ghosts,p2);
180  double uvel1, uvel2, uvelb,
181  vvel1, vvel2, vvelb,
182  wvel1, wvel2, wvelb;
183  _NavierStokes3DGetFlowVar_((phi+nvars*p1),_NavierStokes3D_stride_,rho,uvel1,vvel1,wvel1,energy,pressure,gamma);
184  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel2,vvel2,wvel2,energy,pressure,gamma);
185  uvelb = 1.5*uvel1 - 0.5*uvel2;
186  vvelb = 1.5*vvel1 - 0.5*vvel2;
187  wvelb = 1.5*wvel1 - 0.5*wvel2;
188  double vel_normal = uvelb*nx + vvelb*ny + wvelb*nz;
189 
190  _ArrayCopy1D_(indexb,indexi,ndims);
191  _ArrayAdd1D_(indexi,indexi,boundary->is,ndims);
192  if (face == 1) indexi[dim] = ghosts-1-indexb[dim];
193  else if (face == -1) indexi[dim] = size[dim]-indexb[dim]-1;
194  else return(1);
195  _ArrayIndex1DWO_(ndims,size,indexb,boundary->is,ghosts,p1);
196  _ArrayIndex1D_(ndims,size,indexi,ghosts,p2);
197 
198  /* flow variables in the interior */
199  _NavierStokes3DGetFlowVar_((phi+nvars*p2),_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,gamma);
200 
201  /* set the ghost point values */
202  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
203  if (vel_normal > 0) {
204  /* inflow */
205  rho_gpt = boundary->FlowDensity;
206  pressure_gpt = pressure;
207  uvel_gpt = boundary->FlowVelocity[0];
208  vvel_gpt = boundary->FlowVelocity[1];
209  wvel_gpt = boundary->FlowVelocity[2];
210  } else {
211  /* outflow */
212  rho_gpt = rho;
213  pressure_gpt = boundary->FlowPressure;
214  uvel_gpt = uvel;
215  vvel_gpt = vvel;
216  wvel_gpt = wvel;
217  }
218  energy_gpt = inv_gamma_m1*pressure_gpt
219  + 0.5 * rho_gpt
220  * (uvel_gpt*uvel_gpt + vvel_gpt*vvel_gpt + wvel_gpt*wvel_gpt);
221 
222  phi[nvars*p1+0] = rho_gpt;
223  phi[nvars*p1+1] = rho_gpt * uvel_gpt;
224  phi[nvars*p1+2] = rho_gpt * vvel_gpt;
225  phi[nvars*p1+3] = rho_gpt * wvel_gpt;
226  phi[nvars*p1+4] = energy_gpt;
227 
228  _ArrayIncrementIndex_(ndims,bounds,indexb,done);
229  }
230  }
231 
232  }
233  return(0);
234 }
_NavierStokes3DGetFlowVar_
#define _NavierStokes3DGetFlowVar_(u, stride, rho, vx, vy, vz, e, P, gamma)
Definition: navierstokes3d.h:98
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_Euler2DGetFlowVar_
#define _Euler2DGetFlowVar_(u, rho, vx, vy, e, P, p)
Definition: euler2d.h:44
_ArraySubtract1D_
#define _ArraySubtract1D_(x, a, b, size)
Definition: arrayfunctions.h:201
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
DomainBoundary
Structure containing the variables and function pointers defining a boundary.
Definition: boundaryconditions.h:68
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
_ArrayAdd1D_
#define _ArrayAdd1D_(x, a, b, size)
Definition: arrayfunctions.h:212
Euler2D::gamma
double gamma
Definition: euler2d.h:245
_NavierStokes3D_stride_
static const int _NavierStokes3D_stride_
Definition: navierstokes3d.h:559