HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
NavierStokes3DImmersedBoundary.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <string.h>
8 #include <basic.h>
9 #include <arrayfunctions.h>
10 #include <immersedboundaries.h>
11 #include <physicalmodels/navierstokes3d.h>
12 #include <mpivars.h>
13 #include <hypar.h>
14 
19 int NavierStokes3DIBAdiabatic(void *s,
20  void *m,
21  double *u,
22  double t
23  )
24 {
25  HyPar *solver = (HyPar*) s;
26  MPIVariables *mpi = (MPIVariables*) m;
27  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
28  IBNode *boundary = IB->boundary;
29  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
30  static double v[_MODEL_NVARS_];
31  int n, j, k, nb = IB->n_boundary_nodes;
32 
33  if (!solver->flag_ib) return(0);
34 
35  /* Ideally, this shouldn't be here - But this function is called everywhere
36  (through ApplyIBConditions()) *before* MPIExchangeBoundariesnD is called! */
37  MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,solver->ghosts,mpi,u);
38 
39  double inv_gamma_m1 = 1.0 / (param->gamma - 1.0);
40 
41  double ramp_fac = 1.0;
42  if (param->t_ib_ramp > 0) {
43  double x = t/param->t_ib_ramp;
44  if (!strcmp(param->ib_ramp_type,_IB_RAMP_LINEAR_)) {
45  ramp_fac = x;
46  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_SMOOTHEDSLAB_)) {
47  double a = 0.0;
48  double b = 1.0;
49  double c = 0.5;
50  double r = param->t_ib_ramp/param->t_ib_width;
51  ramp_fac = (a*exp(c*r)+b*exp(r*x))/(exp(c*r)+exp(r*x));
52  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_DISABLE_)) {
53  ramp_fac = 0.0;
54  } else {
55  fprintf(stderr,"Error in NavierStokes3DImmersedBoundary():\n");
56  fprintf(stderr," Ramp type %s not recognized.\n", param->ib_ramp_type);
57  return 1;
58  }
59  }
60 
61  for (n=0; n<nb; n++) {
62 
63  int node_index = boundary[n].p;
64  double *alpha = &(boundary[n].interp_coeffs[0]);
65  int *nodes = &(boundary[n].interp_nodes[0]);
66  double factor = boundary[n].surface_distance / boundary[n].interp_node_distance;
67 
68  _ArraySetValue_(v,_MODEL_NVARS_,0.0);
69  for (j=0; j<_IB_NNODES_; j++) {
70  for (k=0; k<_MODEL_NVARS_; k++) {
71  v[k] += ( alpha[j] * u[_MODEL_NVARS_*nodes[j]+k] );
72  }
73  }
74 
75  double rho, uvel, vvel, wvel, energy, pressure;
76  _NavierStokes3DGetFlowVar_(v,_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,param->gamma);
77 
78  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt;
79  _NavierStokes3DGetFlowVar_( (u+_MODEL_NVARS_*node_index),
80  _NavierStokes3D_stride_,
81  rho_gpt,
82  uvel_gpt,
83  vvel_gpt,
84  wvel_gpt,
85  energy_gpt,
86  pressure_gpt,
87  param->gamma );
88 
89  double rho_ib_target, uvel_ib_target, vvel_ib_target, wvel_ib_target, pressure_ib_target;
90  rho_ib_target = rho;
91  pressure_ib_target = pressure;
92  uvel_ib_target = -uvel * factor;
93  vvel_ib_target = -vvel * factor;
94  wvel_ib_target = -wvel * factor;
95 
96  double rho_ib, uvel_ib, vvel_ib, wvel_ib, energy_ib, pressure_ib;
97  rho_ib = ramp_fac * rho_ib_target + (1.0-ramp_fac) * rho_gpt;
98  pressure_ib = ramp_fac * pressure_ib_target + (1.0-ramp_fac) * pressure_gpt;
99  uvel_ib = ramp_fac * uvel_ib_target + (1.0-ramp_fac) * uvel_gpt;
100  vvel_ib = ramp_fac * vvel_ib_target + (1.0-ramp_fac) * vvel_gpt;
101  wvel_ib = ramp_fac * wvel_ib_target + (1.0-ramp_fac) * wvel_gpt;
102  energy_ib = inv_gamma_m1*pressure_ib
103  + 0.5*rho_ib*(uvel_ib*uvel_ib+vvel_ib*vvel_ib+wvel_ib*wvel_ib);
104 
105  u[_MODEL_NVARS_*node_index+0] = rho_ib;
106  u[_MODEL_NVARS_*node_index+1] = rho_ib * uvel_ib;
107  u[_MODEL_NVARS_*node_index+2] = rho_ib * vvel_ib;
108  u[_MODEL_NVARS_*node_index+3] = rho_ib * wvel_ib;
109  u[_MODEL_NVARS_*node_index+4] = energy_ib;
110  }
111 
112  return 0;
113 }
114 
119 int NavierStokes3DIBIsothermal( void *s,
120  void *m,
121  double *u,
122  double t
123  )
124 {
125  HyPar *solver = (HyPar*) s;
126  MPIVariables *mpi = (MPIVariables*) m;
127  ImmersedBoundary *IB = (ImmersedBoundary*) solver->ib;
128  IBNode *boundary = IB->boundary;
129  NavierStokes3D *param = (NavierStokes3D*) solver->physics;
130  static double v[_MODEL_NVARS_];
131  int n, j, k, nb = IB->n_boundary_nodes;
132 
133  if (!solver->flag_ib) return(0);
134 
135  /* Ideally, this shouldn't be here - But this function is called everywhere
136  (through ApplyIBConditions()) *before* MPIExchangeBoundariesnD is called! */
137  MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,solver->dim_local,solver->ghosts,mpi,u);
138 
139  double inv_gamma_m1 = 1.0 / (param->gamma - 1.0);
140 
141  double ramp_fac = 1.0;
142  if (param->t_ib_ramp > 0) {
143  double x = t/param->t_ib_ramp;
144  if (!strcmp(param->ib_ramp_type,_IB_RAMP_LINEAR_)) {
145  ramp_fac = x;
146  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_SMOOTHEDSLAB_)) {
147  double a = 0.0;
148  double b = 1.0;
149  double c = 0.5;
150  double r = param->t_ib_ramp/param->t_ib_width;
151  ramp_fac = (a*exp(c*r)+b*exp(r*x))/(exp(c*r)+exp(r*x));
152  } else if (!strcmp(param->ib_ramp_type,_IB_RAMP_DISABLE_)) {
153  ramp_fac = 0.0;
154  } else {
155  fprintf(stderr,"Error in NavierStokes3DImmersedBoundary():\n");
156  fprintf(stderr," Ramp type %s not recognized.\n", param->ib_ramp_type);
157  return 1;
158  }
159  }
160 
161  for (n=0; n<nb; n++) {
162 
163  int node_index = boundary[n].p;
164  double *alpha = &(boundary[n].interp_coeffs[0]);
165  int *nodes = &(boundary[n].interp_nodes[0]);
166  double factor = boundary[n].surface_distance / boundary[n].interp_node_distance;
167 
168  _ArraySetValue_(v,_MODEL_NVARS_,0.0);
169  for (j=0; j<_IB_NNODES_; j++) {
170  for (k=0; k<_MODEL_NVARS_; k++) {
171  v[k] += ( alpha[j] * u[_MODEL_NVARS_*nodes[j]+k] );
172  }
173  }
174 
175  double rho, uvel, vvel, wvel, energy, pressure, temperature;
176  _NavierStokes3DGetFlowVar_(v,_NavierStokes3D_stride_,rho,uvel,vvel,wvel,energy,pressure,param->gamma);
177  temperature = pressure / rho;
178 
179  double rho_gpt, uvel_gpt, vvel_gpt, wvel_gpt, energy_gpt, pressure_gpt, temperature_gpt;
180  _NavierStokes3DGetFlowVar_( (u+_MODEL_NVARS_*node_index),
181  _NavierStokes3D_stride_,
182  rho_gpt,
183  uvel_gpt,
184  vvel_gpt,
185  wvel_gpt,
186  energy_gpt,
187  pressure_gpt,
188  param->gamma );
189  temperature_gpt = pressure_gpt / rho_gpt;
190 
191  double rho_ib_target,
192  uvel_ib_target, vvel_ib_target, wvel_ib_target,
193  pressure_ib_target,
194  temperature_ib_target;
195  temperature_ib_target = (1.0+factor)*param->T_ib_wall - factor * temperature;
196  if ( (temperature_ib_target < param->T_ib_wall/param->ib_T_tol)
197  || (temperature_ib_target > param->T_ib_wall*param->ib_T_tol) ) {
198  temperature_ib_target = param->T_ib_wall;
199  }
200  pressure_ib_target = pressure;
201  rho_ib_target = pressure_ib_target / temperature_ib_target;
202  uvel_ib_target = - factor * uvel;
203  vvel_ib_target = - factor * vvel;
204  wvel_ib_target = - factor * wvel;
205 
206  double rho_ib, uvel_ib, vvel_ib, wvel_ib, energy_ib, pressure_ib;
207  rho_ib = ramp_fac * rho_ib_target + (1.0-ramp_fac) * rho_gpt;
208  pressure_ib = ramp_fac * pressure_ib_target + (1.0-ramp_fac) * pressure_gpt;
209  uvel_ib = ramp_fac * uvel_ib_target + (1.0-ramp_fac) * uvel_gpt;
210  vvel_ib = ramp_fac * vvel_ib_target + (1.0-ramp_fac) * vvel_gpt;
211  wvel_ib = ramp_fac * wvel_ib_target + (1.0-ramp_fac) * wvel_gpt;
212  energy_ib = inv_gamma_m1*pressure_ib
213  + 0.5*rho_ib*(uvel_ib*uvel_ib+vvel_ib*vvel_ib+wvel_ib*wvel_ib);
214 
215  u[_MODEL_NVARS_*node_index+0] = rho_ib;
216  u[_MODEL_NVARS_*node_index+1] = rho_ib * uvel_ib;
217  u[_MODEL_NVARS_*node_index+2] = rho_ib * vvel_ib;
218  u[_MODEL_NVARS_*node_index+3] = rho_ib * wvel_ib;
219  u[_MODEL_NVARS_*node_index+4] = energy_ib;
220  }
221 
222  return(0);
223 }
_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
navierstokes3d.h
3D Navier Stokes equations (compressible flows)
IBNode
Structure defining an immersed boundary node.
Definition: immersedboundaries.h:125
ImmersedBoundary
Structure containing variables for immersed boundary implementation.
Definition: immersedboundaries.h:153
immersedboundaries.h
Structures and function definitions for immersed boundaries.
HyPar::ib
void * ib
Definition: hypar.h:443
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_IB_RAMP_SMOOTHEDSLAB_
#define _IB_RAMP_SMOOTHEDSLAB_
Definition: navierstokes3d.h:88
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::flag_ib
int flag_ib
Definition: hypar.h:441
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
mpivars.h
MPI related function definitions.
_MODEL_NDIMS_
#define _MODEL_NDIMS_
Definition: euler1d.h:56
_IB_RAMP_LINEAR_
#define _IB_RAMP_LINEAR_
Definition: navierstokes3d.h:86
HyPar::ghosts
int ghosts
Definition: hypar.h:52
hypar.h
Contains structure definition for hypar.
NavierStokes3DIBIsothermal
int NavierStokes3DIBIsothermal(void *s, void *m, double *u, double t)
Definition: NavierStokes3DImmersedBoundary.c:119
_IB_RAMP_DISABLE_
#define _IB_RAMP_DISABLE_
Definition: navierstokes3d.h:90
NavierStokes3DIBAdiabatic
int NavierStokes3DIBAdiabatic(void *s, void *m, double *u, double t)
Definition: NavierStokes3DImmersedBoundary.c:19
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
basic.h
Some basic definitions and macros.
arrayfunctions.h
Contains macros and function definitions for common array operations.
MPIExchangeBoundariesnD
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Definition: MPIExchangeBoundariesnD.c:42
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_IB_NNODES_
#define _IB_NNODES_
Definition: immersedboundaries.h:13
_NavierStokes3D_stride_
static const int _NavierStokes3D_stride_
Definition: navierstokes3d.h:559