HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
BCSubsonicAmbivalent.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <basic.h>
8 #include <arrayfunctions.h>
9 #include <boundaryconditions.h>
10 
11 #include <physicalmodels/euler2d.h>
12 #include <physicalmodels/navierstokes3d.h>
13 
28 int BCSubsonicAmbivalentU(
29  void *b,
30  void *m,
31  int ndims,
32  int nvars,
33  int *size,
34  int ghosts,
35  double *phi,
36  double waqt
37  )
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
BCSubsonicAmbivalentU
int BCSubsonicAmbivalentU(void *, void *, int, int, int *, int, double *, double)
Definition: BCSubsonicAmbivalent.c:28
navierstokes3d.h
3D Navier Stokes equations (compressible flows)
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
boundaryconditions.h
Containts the structures and definitions for boundary condition implementation.
_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
basic.h
Some basic definitions and macros.
arrayfunctions.h
Contains macros and function definitions for common array operations.
_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