TYCHO  1.3.0
 All Data Structures Files Functions Variables Enumerations Enumerator
/home/kapf/tycho_docu/make_ic.c
Go to the documentation of this file.
1 /*
2  * make_ic.c
3  *
4  * Author: Wolfgang Kapferer
5  */
6 
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <math.h>
10 
11 #include "variables_global.h"
12 #include "prototypes.h"
13 
19 int make_ic(int x, int y, int z) {
20  int i, j, k;
21  double velo_x, velo_y, velo_z;
22 
23  double T0, T, A, rho0, pre0;
24 
25  //T(z)=T_0 + A*(z-z0)
26  T0 = 288.15; //[K]
27  A = -6.5E-3; //[K/m]
28 
29  rho0 = 1.229; // Density at sea level [kg/m^3]
30  pre0 = 1.013E5; //pressure at sea level [N/m^2]
31 
32  velo_x = small;
33  velo_y = small;
34  velo_z = small;
35 
36  if (starting_flow == 1) {
37  if (bound.down == 4) velo_y = inflow_velocity;
38  if (bound.up == 4) velo_y = inflow_velocity;
39  if (bound.left == 4) velo_x = inflow_velocity;
40  if (bound.right == 4) velo_x = inflow_velocity;
41  if (bound.front == 4) velo_z = inflow_velocity;
42  if (bound.back == 4) velo_z = inflow_velocity;
43  }
44 
45  if (strat_const_atmos == 1) {
46  for (i = 0; i < x; i++) {
47  for (j = 0; j < y; j++) {
48  for (k = 0; k < z; k++) {
49 
50  T = T0 + A * ((ymax / y) * (y - j));
51 
52  rho[i][j][k] = rho0 * pow((T / T0), (-grav_acc / (gasconstant * A) + 1));
53  pre[i][j][k] = pre0 * pow((T / T0), (-grav_acc / (gasconstant * A)));
54 
55  vx[i][j][k] = velo_x;
56  vy[i][j][k] = velo_y;
57  vz[i][j][k] = velo_z;
58 
59  }
60  }
61  }
62  }
63 
64  if (strat_const_atmos == 0) {
65  for (i = 0; i < x; i++) {
66  for (j = 0; j < y; j++) {
67  for (k = 0; k < z; k++) {
68 
69  rho[i][j][k] = 1.229;
70  pre[i][j][k] = 300 * gasconstant * rho[i][j][k];
71 
72  vx[i][j][k] = velo_x;
73  vy[i][j][k] = velo_y;
74  vz[i][j][k] = velo_z;
75 
76  }
77  }
78  }
79  }
80 
81  // If you need a conversion from temperature to pressure
82  // calculate_pressure(x, y, z, temperature);
83 
84  //just to add explosions
85  //insert_pressure(x, y, z);
86 
87  return 0;
88 }
89 
93 int calculate_pressure(int x, int y, int z, double temperature) {
94  int i, j, k;
95 
96 
97  for (i = 0; i < x; i++) {
98  for (j = 0; j < y; j++) {
99  for (k = 0; k < z; k++) {
100  pre[i][j][k] = temperature * gasconstant * rho[i][j][k];
101  }
102  }
103  }
104 
105  return 0;
106 }
107 
112 int initiate_domain(int x, int y, int z) {
113  int i, j, k;
114 
115  for (i = 0; i < x; i++) {
116  for (j = 0; j < y; j++) {
117  for (k = 0; k < z; k++) {
118  int i2 = i - x / 2.0;
119  int j2 = j - y / 2.0;
120  int xy = sqrt(i2 * i2 + j2 * j2);
121 
122  if (xy < x / 20) {
123  dom[i][j][k] = 1;
124  rho[i][j][k] = obstacle_density;
125  pre[i][j][k] = obstacle_temperature;
126 
127  vx[i][j][k] = 0.0;
128  vy[i][j][k] = 0.0;
129  vz[i][j][k] = 0.0;
130  } else {
131  dom[i][j][k] = 0;
132  }
133 
134  }
135  }
136  }
137 
138  printf("Initiate domain done\n");
139 
140  return 0;
141 }
142 
147 int initiate_domain_dom(int x, int y, int z) {
148  int i, j, k;
149 
150  for (i = 0; i < x; i++) {
151  for (j = 0; j < y; j++) {
152  for (k = 0; k < z; k++) {
153 
154  if (dom[i][j][k] == 1) {
155  rho[i][j][k] = obstacle_density;
156  pre[i][j][k] = obstacle_temperature;
157 
158  vx[i][j][k] = 0.0;
159  vy[i][j][k] = 0.0;
160  vz[i][j][k] = 0.0;
161  }
162  }
163  }
164  }
165 
166  printf("Initiate domain done\n");
167 
168  return 0;
169 }
170 
174 int initiate_domain_marker(int x, int y, int z) {
175  int i, j, k;
176 
177  for (i = 0; i < x; i++) {
178  for (j = 0; j < y; j++) {
179  for (k = 0; k < z; k++) {
180 
181  if (marker[i][j][k] == 1) {
182  marker[i][j][k] = marker_density;
183  }
184  }
185  }
186  }
187 
188  printf("Initiate domain done\n");
189 
190  return 0;
191 }
192 
196 int make_sod_ic(int x, int y, int z) {
197  int i, j, k;
198  float radius, radius1;
199 
200  if (dimension == 2) radius = sqrt(pow(x, 2) + pow(y, 2));
201  if (dimension == 3) radius = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
202 
203  for (i = 0; i < x; i++) {
204  for (j = 0; j < y; j++) {
205  for (k = 0; k < z; k++) {
206 
207  if (dimension == 1) {
208  if (i < x / 2) {
209  rho[i][j][k] = 1.0;
210  pre[i][j][k] = 1.0 * 5E4;
211  vx[i][j][k] = 0.0;
212  vy[i][j][k] = 0.0;
213  vz[i][j][k] = 0.0;
214  } else {
215  rho[i][j][k] = 0.15;
216  pre[i][j][k] = 0.10 * 5E4;
217  vx[i][j][k] = 0.0;
218  vy[i][j][k] = 0.0;
219  vz[i][j][k] = 0.0;
220  }
221  }
222 
223  if (dimension == 2) {
224 
225  radius1 = sqrt(pow(i, 2) + pow(j, 2));
226 
227  if (radius1 < radius / 2) {
228  rho[i][j][k] = 1.0;
229  pre[i][j][k] = 1.0 * 1E4;
230  vx[i][j][k] = 0.0;
231  vy[i][j][k] = 0.0;
232  vz[i][j][k] = 0.0;
233  } else {
234  rho[i][j][k] = 0.15;
235  pre[i][j][k] = 0.10 * 1E2;
236  vx[i][j][k] = 0.0;
237  vy[i][j][k] = 0.0;
238  vz[i][j][k] = 0.0;
239  }
240  }
241 
242  if (dimension == 3) {
243 
244  radius1 = sqrt(pow(i, 2) + pow(j, 2) + pow(k, 2));
245 
246  if (radius1 <= radius / 2) {
247  rho[i][j][k] = 1.0;
248  pre[i][j][k] = 1.0;
249  vx[i][j][k] = 0.0;
250  vy[i][j][k] = 0.0;
251  vz[i][j][k] = 0.0;
252  } else {
253  rho[i][j][k] = 0.15;
254  pre[i][j][k] = 0.10;
255  vx[i][j][k] = 0.0;
256  vy[i][j][k] = 0.0;
257  vz[i][j][k] = 0.0;
258  }
259  }
260 
261  }
262  }
263  }
264 
265  return 0;
266 }
267 
271 int make_kh_instabilities(int x, int y, int z) {
272 
273  int i, j, k;
274 
275  for (i = 0; i < x; i++) {
276  for (j = 0; j < y; j++) {
277  for (k = 0; k < z; k++) {
278 
279  if ((j < y / 4) || (j > 3 * y / 4)) {
280  rho[i][j][k] = 1;
281  pre[i][j][k] = 1E5;
282  vx[i][j][k] = -50;
283  vy[i][j][k] = 0.0;
284  vz[i][j][k] = 0.0;
285  } else {
286  rho[i][j][k] = 2;
287  pre[i][j][k] = 1E5;
288  vx[i][j][k] = 0.0;
289  vy[i][j][k] = 0.0;
290  vz[i][j][k] = 0.0;
291  }
292 
293  if (dimension == 3) {
294  if (j == y / 4) vy[i][j][k] = sin(xmax / x * i) * cos(zmax / z * k);
295  if (j == 3 * y / 4) vy[i][j][k] = -1 * sin(xmax / x * i) * cos(zmax / z * k);
296  }
297 
298  if (dimension == 2) {
299  if (j == y / 4) vy[i][j][k] = sin(xmax / x * i);
300  if (j == 3 * y / 4) vy[i][j][k] = -1 * sin(xmax / x * i);
301  }
302  }
303  }
304  }
305 
306 
307 
308  return 0;
309 }
310 
314 int insert_pressure(int x, int y, int z) {
315 
316  int i, j, k;
317 
318  for (i = 0; i < x; i++) {
319  for (j = 0; j < y; j++) {
320  for (k = 0; k < z; k++) {
321  if (dom[i][j][k] != 1) {
322 
323  if (j < y / 4) {
324  if (round((double) rand() / (double) RAND_MAX) == 1) {
325  pre[i][j][k] = 2000 * gasconstant * rho[i][j][k];
326  }
327  }
328  }
329  }
330  }
331  }
332 
333 
334  return 0;
335 }
Generated on Thu Oct 10 2013 17:15:52 for TYCHO by   doxygen 1.8.1.1