TYCHO  1.3.0
 All Data Structures Files Functions Variables Enumerations Enumerator
/home/kapf/tycho_docu/diffusion.c
Go to the documentation of this file.
1 /*
2  * diffusion.c
3  *
4  * Author: Wolfgang Kapferer
5  */
6 
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <math.h>
10 
11 #ifdef _OPENMP
12 #include <omp.h>
13 #endif
14 
15 #include "variables_global.h"
16 #include "prototypes.h"
17 
25 int diffusion(int x, int y, int z) {
26 
27  int i, j, k, n;
28  double dxfac, dyfac, dzfac;
29  double coefficient_solid;
30  double alpha_solid_x, alpha_solid_y, alpha_solid_z;
31  double *temperature_1D, *temperature_1D_future;
32  double *mass;
33  double energy_left, energy_right;
34  double temperature_right_gas, temerpature_left_gas;
35  double A;
36  double heat_transfer_coefficient;
37  double volume;
38  int *dom_1D;
39  int dom_counter;
40  int i_local_min, i_local_max;
41  int condition;
42 
43  dom_counter = 0;
44  i_local_min = 0;
45  i_local_max = 0;
46  temperature_right_gas = 0.0;
47  temerpature_left_gas = 0.0;
48  A = 0.0;
49  heat_transfer_coefficient = 0.0;
50  condition = 0;
51 
52  // spacing in all three directions
53  dxfac = (xmax - xmin) / (double) x;
54  dyfac = (ymax - ymin) / (double) y;
55  dzfac = (zmax - zmin) / (double) z;
56  if (dimension == 1) volume = dxfac;
57  if (dimension == 2) volume = dxfac * dyfac;
58  if (dimension == 3) volume = dxfac * dyfac * dzfac;
59 
60 
61  dxfac = dxfac*dxfac;
62  dyfac = dyfac*dyfac;
63  dzfac = dzfac*dzfac;
64 
65  coefficient_solid = obstacle_heat_conductivity;
66 
67  if (dimension > 0) alpha_solid_x = coefficient_solid * dt / (dxfac);
68  if (dimension > 1) alpha_solid_y = coefficient_solid * dt / (dyfac);
69  if (dimension > 2) alpha_solid_z = coefficient_solid * dt / (dzfac);
70 
71 
72  if ((alpha_solid_x >= 0.3) && (dimension > 0)) {
73  dt = (0.3 * dxfac) / coefficient_solid;
74  }
75  if ((alpha_solid_y >= 0.3) && (dimension > 1)) {
76  dt = (0.3 * dyfac) / coefficient_solid;
77  }
78  if ((alpha_solid_z >= 0.3) && (dimension > 2)) {
79  dt = (0.3 * dzfac) / coefficient_solid;
80  }
81 
82 
83  //THE X-SWEEP=========================================================================================
84 #ifdef _OPENMP
85 #pragma omp parallel default(none) \
86  private(j, i, k, n, temperature_1D, dom_1D, \
87  temperature_1D_future, mass, dom_counter, \
88  i_local_min, i_local_max, energy_left, energy_right, \
89  condition, heat_transfer_coefficient) \
90  shared(x, y, z, dom, rho, pre, dt, dxfac, \
91  dyfac, dzfac, alpha_solid_x, gasconstant, A, \
92  dimension, volume, \
93  specific_heat_capacity_gas, \
94  specific_heat_capacity_obstacle)
95  {
96 #endif
97 
98  init_diffusion(&temperature_1D, &temperature_1D_future, &mass, &dom_1D);
99 
100  if (dimension == 1) A = 1;
101  if (dimension == 2) A = dyfac;
102  if (dimension == 3) A = dyfac * dzfac;
103 
104 
105  //the x-sweep
106  for (k = 0; k < z; k++) {
107 
108 
109 #ifdef _OPENMP
110 #pragma omp for schedule(static)
111 #endif
112 
113  for (j = 0; j < y; j++) {
114 
115  condition = 0;
116  i_local_min = 0;
117  i_local_max = 0;
118 
119  while (condition == 0) {
120 
121  dom_counter = 0;
122  energy_left = 0.0;
123  energy_right = 0.0;
124 
125  for (i = i_local_min; i < x; i++) {
126 
127  temperature_1D[i] = pre[i][j][k] / (gasconstant * rho[i][j][k]);
128  mass[i] = rho[i][j][k] * volume;
129  dom_1D[i] = dom[i][j][k];
130 
131  if (dom_1D[i] == 1) {
132  //Now we have the beginning- and the end-index of the solid rod
133  if (dom_counter == 0) i_local_min = i;
134  dom_counter++;
135  } else {
136  if (dom_counter > 0) {
137  i_local_max = i - 1;
138  break;
139  }
140  }
141  if (i == (x - 1)) {
142  i_local_max = i;
143  if (dom_counter == 0) i_local_min = i;
144  }
145  }
146 
147  if (i == x) {
148  condition = 1;
149  }
150 
151  //in the boundary of the obstacle
152  //case if i_local_min > 0
153  if (i_local_min != 0) {
154  heat_transfer_coefficient = alpha_heat_transfer(i_local_min - 1, j, k, 0);
155  energy_left = heat_transfer_coefficient * A * (temperature_1D[i_local_min - 1] - temperature_1D[i_local_min]) * dt;
156  }
157  if (i_local_max != (x - 1)) {
158  heat_transfer_coefficient = alpha_heat_transfer(i_local_max + 1, j, k, 0);
159  energy_right = heat_transfer_coefficient * A * (temperature_1D[i_local_max + 1] - temperature_1D[i_local_max]) * dt;
160  }
161 
162  //temperature in the obstacle higher than in the gas left of it
163  if ((energy_left <= 0) && (i_local_min != 0)) {
164  temperature_1D[i_local_min - 1] = temperature_1D[i_local_min - 1] + fabs(energy_left / (specific_heat_capacity_gas * mass[i_local_min - 1]));
165  temperature_1D[i_local_min] = temperature_1D[i_local_min] - fabs(energy_left / (specific_heat_capacity_obstacle * mass[i_local_min]));
166 
167  pre[i_local_min - 1][j][k] = (gasconstant * rho[i_local_min - 1][j][k]) * temperature_1D[i_local_min - 1];
168  pre[i_local_min][j][k] = (gasconstant * rho[i_local_min][j][k]) * temperature_1D[i_local_min];
169 
170  }
171 
172  //temperature in the obstacle lower than in the gas left of it
173  if ((energy_left > 0) && (i_local_min != 0)) {
174  temperature_1D[i_local_min - 1] = temperature_1D[i_local_min - 1] - fabs(energy_left / (specific_heat_capacity_gas * mass[i_local_min - 1]));
175  temperature_1D[i_local_min] = temperature_1D[i_local_min] + fabs(energy_left / (specific_heat_capacity_obstacle * mass[i_local_min]));
176 
177  pre[i_local_min - 1][j][k] = (gasconstant * rho[i_local_min - 1][j][k]) * temperature_1D[i_local_min - 1];
178  pre[i_local_min][j][k] = (gasconstant * rho[i_local_min][j][k]) * temperature_1D[i_local_min];
179 
180  }
181 
182  //temperature in the gas is higher than in the gas right of it
183  if ((energy_right > 0) && (i_local_max != (x - 1))) {
184  temperature_1D[i_local_max] = temperature_1D[i_local_max] + fabs(energy_right / (specific_heat_capacity_obstacle * mass[i_local_max]));
185  temperature_1D[i_local_max + 1] = temperature_1D[i_local_max + 1] - fabs(energy_right / (specific_heat_capacity_gas * mass[i_local_max + 1]));
186 
187  pre[i_local_max][j][k] = (gasconstant * rho[i_local_max][j][k]) * temperature_1D[i_local_max];
188  pre[i_local_max + 1][j][k] = (gasconstant * rho[i_local_max + 1][j][k]) * temperature_1D[i_local_max + 1];
189 
190  }
191 
192  //temperature in the gas lower than in the obstacle left of it
193  if ((energy_right <= 0) && (i_local_max != (x - 1))) {
194  temperature_1D[i_local_max] = temperature_1D[i_local_max] - fabs(energy_right / (specific_heat_capacity_obstacle * mass[i_local_max]));
195  temperature_1D[i_local_max + 1] = temperature_1D[i_local_max + 1] + fabs(energy_right / (specific_heat_capacity_gas * mass[i_local_max + 1]));
196 
197  pre[i_local_max][j][k] = (gasconstant * rho[i_local_max][j][k]) * temperature_1D[i_local_max];
198  pre[i_local_max + 1][j][k] = (gasconstant * rho[i_local_max + 1][j][k]) * temperature_1D[i_local_max + 1];
199 
200  }
201 
202  if (dom_counter > 2) {
203  for (n = i_local_min + 1; n < i_local_max; n++) {
204  temperature_1D_future[n] = alpha_solid_x * temperature_1D[n + 1] +
205  (1 - 2 * alpha_solid_x) * temperature_1D[n] +
206  alpha_solid_x * temperature_1D[n - 1];
207  }
208 
209  for (i = i_local_min + 1; i < i_local_max; i++) {
210  pre[i][j][k] = (gasconstant * rho[i][j][k]) * temperature_1D_future[i];
211  }
212  }
213 
214  for (i = i_local_min; i < i_local_max; i++) {
215  temperature_1D[i] = 0.0;
216  temperature_1D_future[i] = 0.0;
217  mass[i] = 0.0;
218  dom_1D[i] = 0;
219  }
220 
221  i_local_min = i_local_max + 1;
222  }
223  }
224  }
225  free_diffusion(&temperature_1D, &temperature_1D_future, &mass, &dom_1D);
226 
227 #ifdef _OPENMP
228  }
229 #endif
230 
231 
232  //THE Y-SWEEP=========================================================================================
233 #ifdef _OPENMP
234 #pragma omp parallel default(none) \
235  private(j, i, k, n, temperature_1D, dom_1D, \
236  temperature_1D_future, mass, dom_counter, \
237  i_local_min, i_local_max, energy_left, \
238  condition, energy_right, \
239  heat_transfer_coefficient) \
240  shared(x, y, z, dom, rho, pre, dt, dxfac, \
241  dyfac, dzfac, alpha_solid_y, \
242  gasconstant, A, dimension, volume, \
243  specific_heat_capacity_gas, \
244  specific_heat_capacity_obstacle)
245  {
246 #endif
247 
248  init_diffusion(&temperature_1D, &temperature_1D_future, &mass, &dom_1D);
249 
250  if (dimension == 1) A = 1;
251  if (dimension == 2) A = dxfac;
252  if (dimension == 3) A = dxfac * dzfac;
253 
254  //the y-sweep
255  for (k = 0; k < z; k++) {
256 
257 #ifdef _OPENMP
258 #pragma omp for schedule(static)
259 #endif
260 
261  for (i = 0; i < x; i++) {
262 
263  condition = 0;
264  i_local_min = 0;
265  i_local_max = 0;
266 
267  while (condition == 0) {
268 
269  dom_counter = 0;
270  energy_left = 0.0;
271  energy_right = 0.0;
272 
273 
274  for (j = i_local_min; j < y; j++) {
275 
276  temperature_1D[j] = pre[i][j][k] / (gasconstant * rho[i][j][k]);
277  mass[j] = rho[i][j][k] * volume;
278  dom_1D[j] = dom[i][j][k];
279 
280  if (dom_1D[j] == 1) {
281  //Now we have the beginning- and the end-index of the solid rod
282  if (dom_counter == 0) i_local_min = j;
283  dom_counter++;
284  } else {
285  if (dom_counter > 0) {
286  i_local_max = j - 1;
287  break;
288  }
289  }
290  if (j == (y - 1)) {
291  i_local_max = j;
292  if (dom_counter == 0) i_local_min = j;
293  }
294  }
295 
296  if (j == y) {
297  condition = 1;
298  }
299 
300  //Now the heat transfer is calcuated to set the the right temerpature
301  //in the boundary of the obstacle
302  //case if i_local_min > 0
303  if (i_local_min != 0) {
304  heat_transfer_coefficient = alpha_heat_transfer(i, i_local_min - 1, k, 1);
305  energy_left = heat_transfer_coefficient * A * (temperature_1D[i_local_min - 1] - temperature_1D[i_local_min]) * dt;
306  }
307  if (i_local_max != (y - 1)) {
308  heat_transfer_coefficient = alpha_heat_transfer(i, i_local_max + 1, k, 1);
309  energy_right = heat_transfer_coefficient * A * (temperature_1D[i_local_max + 1] - temperature_1D[i_local_max]) * dt;
310  }
311 
312  //temperature in the obstacle higher than in the gas left of it
313  if ((energy_left <= 0) && (i_local_min != 0)) {
314  temperature_1D[i_local_min - 1] = temperature_1D[i_local_min - 1] + fabs(energy_left / (specific_heat_capacity_gas * mass[i_local_min - 1]));
315  temperature_1D[i_local_min] = temperature_1D[i_local_min] - fabs(energy_left / (specific_heat_capacity_obstacle * mass[i_local_min]));
316 
317  pre[i][i_local_min - 1][k] = (gasconstant * rho[i][i_local_min - 1][k]) * temperature_1D[i_local_min - 1];
318  pre[i][i_local_min][k] = (gasconstant * rho[i][i_local_min][k]) * temperature_1D[i_local_min];
319  }
320 
321  //temperature in the obstacle lower than in the gas left of it
322  if ((energy_left > 0) && (i_local_min != 0)) {
323  temperature_1D[i_local_min - 1] = temperature_1D[i_local_min - 1] - fabs(energy_left / (specific_heat_capacity_gas * mass[i_local_min - 1]));
324  temperature_1D[i_local_min] = temperature_1D[i_local_min] + fabs(energy_left / (specific_heat_capacity_obstacle * mass[i_local_min]));
325 
326  pre[i][i_local_min - 1][k] = (gasconstant * rho[i][i_local_min - 1][k]) * temperature_1D[i_local_min - 1];
327  pre[i][i_local_min][k] = (gasconstant * rho[i][i_local_min][k]) * temperature_1D[i_local_min];
328  }
329 
330  //temperature in the obstacle higher than in the gas right of it
331  if ((energy_right > 0) && (i_local_max != (y - 1))) {
332  temperature_1D[i_local_max] = temperature_1D[i_local_max] + fabs(energy_right / (specific_heat_capacity_obstacle * mass[i_local_max]));
333  temperature_1D[i_local_max + 1] = temperature_1D[i_local_max + 1] - fabs(energy_right / (specific_heat_capacity_gas * mass[i_local_max + 1]));
334 
335  pre[i][i_local_max][k] = (gasconstant * rho[i][i_local_max][k]) * temperature_1D[i_local_max];
336  pre[i][i_local_max + 1][k] = (gasconstant * rho[i][i_local_max + 1][k]) * temperature_1D[i_local_max + 1];
337  }
338 
339  //temperature in the obstacle lower than in the gas right of it
340  if ((energy_right <= 0) && (i_local_max != (y - 1))) {
341  temperature_1D[i_local_max] = temperature_1D[i_local_max] - fabs(energy_right / (specific_heat_capacity_obstacle * mass[i_local_max]));
342  temperature_1D[i_local_max + 1] = temperature_1D[i_local_max + 1] + fabs(energy_right / (specific_heat_capacity_gas * mass[i_local_max + 1]));
343 
344  pre[i][i_local_max][k] = (gasconstant * rho[i][i_local_max][k]) * temperature_1D[i_local_max];
345  pre[i][i_local_max + 1][k] = (gasconstant * rho[i][i_local_max + 1][k]) * temperature_1D[i_local_max + 1];
346  }
347 
348  //printf("energy_left: %g energy_right: %g dt: %g\n", energy_left, energy_right, dt);
349 
350  if (dom_counter > 2) {
351 
352  for (n = i_local_min + 1; n < i_local_max; n++) {
353  //here the heat conductivity constant are determined
354  temperature_1D_future[n] = alpha_solid_y * temperature_1D[n + 1] +
355  (1 - 2 * alpha_solid_y) * temperature_1D[n] +
356  alpha_solid_y * temperature_1D[n - 1];
357  }
358 
359  for (j = i_local_min + 1; j < i_local_max; j++) {
360  pre[i][j][k] = (gasconstant * rho[i][j][k]) * temperature_1D_future[j];
361  }
362 
363 
364  for (j = i_local_min; j < i_local_max; j++) {
365  temperature_1D[j] = 0.0;
366  temperature_1D_future[j] = 0.0;
367  mass[j] = 0.0;
368  dom_1D[j] = 0;
369  }
370  }
371 
372  i_local_min = i_local_max + 1;
373 
374  }
375  }
376  }
377  free_diffusion(&temperature_1D, &temperature_1D_future, &mass, &dom_1D);
378 
379 #ifdef _OPENMP
380  }
381 #endif
382 
383  //THE Z-SWEEP=========================================================================================
384 #ifdef _OPENMP
385 #pragma omp parallel default(none) \
386  private(j, i, k, n, temperature_1D, dom_1D, \
387  temperature_1D_future, mass, dom_counter, \
388  i_local_min, i_local_max, energy_left, \
389  condition, energy_right, \
390  heat_transfer_coefficient) \
391  shared(x, y, z, dom, rho, pre, dt, dxfac, \
392  dyfac, dzfac, alpha_solid_z, gasconstant, A, \
393  dimension, volume, \
394  specific_heat_capacity_gas, \
395  specific_heat_capacity_obstacle)
396  {
397 #endif
398  init_diffusion(&temperature_1D, &temperature_1D_future, &mass, &dom_1D);
399 
400  if (dimension == 1) A = 1;
401  if (dimension == 2) A = dxfac;
402  if (dimension == 3) A = dxfac * dyfac;
403 
404  //the z-sweep
405  for (j = 0; j < y; j++) {
406 
407 #ifdef _OPENMP
408 #pragma omp for schedule(static)
409 #endif
410 
411  for (i = 0; i < x; i++) {
412 
413  condition = 0;
414  i_local_min = 0;
415  i_local_max = 0;
416 
417  while (condition == 0) {
418 
419  dom_counter = 0;
420  energy_left = 0.0;
421  energy_right = 0.0;
422 
423  for (k = i_local_min; k < z; k++) {
424 
425  temperature_1D[k] = pre[i][j][k] / (gasconstant * rho[i][j][k]);
426  mass[k] = rho[i][j][k] * volume;
427  dom_1D[k] = dom[i][j][k];
428 
429  if (dom_1D[k] == 1) {
430  //Now we have the beginning- and the end-index of the solid rod
431  if (dom_counter == 0) i_local_min = k;
432  dom_counter++;
433  } else {
434  if (dom_counter > 0) {
435  i_local_max = k - 1;
436  break;
437  }
438  }
439  if (k == (z - 1)) {
440  i_local_max = k;
441  if (dom_counter == 0) i_local_min = k;
442  }
443  }
444 
445  if (k == z) {
446  condition = 1;
447  }
448 
449  //Now the heat transfer is calculated to set the the right temperature
450  //in the boundary of the obstacle
451  //case if i_local_min > 0
452  if (i_local_min != 0) {
453  heat_transfer_coefficient = alpha_heat_transfer(i, j, i_local_min - 1, 2);
454  energy_left = heat_transfer_coefficient * A * (temperature_1D[i_local_min - 1] - temperature_1D[i_local_min]) * dt;
455  }
456  if (i_local_max != (z - 1)) {
457  heat_transfer_coefficient = alpha_heat_transfer(i, j, i_local_max + 1, 2);
458  energy_right = heat_transfer_coefficient * A * (temperature_1D[i_local_max + 1] - temperature_1D[i_local_max]) * dt;
459  }
460 
461  //temperature in the obstacle higher than in the gas left of it
462  if ((energy_left <= 0) && (i_local_min != 0)) {
463  temperature_1D[i_local_min - 1] = temperature_1D[i_local_min - 1] + fabs(energy_left / (specific_heat_capacity_gas * mass[i_local_min - 1]));
464  temperature_1D[i_local_min] = temperature_1D[i_local_min] - fabs(energy_left / (specific_heat_capacity_obstacle * mass[i_local_min]));
465 
466  pre[i][j][i_local_min - 1] = (gasconstant * rho[i][j][i_local_min - 1]) * temperature_1D[i_local_min - 1];
467  pre[i][j][i_local_min] = (gasconstant * rho[i][j][i_local_min]) * temperature_1D[i_local_min ];
468  }
469 
470  //temperature in the obstacle lower than in the gas left of it
471  if ((energy_left > 0) && (i_local_min != 0)) {
472  temperature_1D[i_local_min - 1] = temperature_1D[i_local_min - 1] - fabs(energy_left / (specific_heat_capacity_gas * mass[i_local_min - 1]));
473  temperature_1D[i_local_min] = temperature_1D[i_local_min] + fabs(energy_left / (specific_heat_capacity_obstacle * mass[i_local_min]));
474 
475  pre[i][j][i_local_min - 1] = (gasconstant * rho[i][j][i_local_min - 1]) * temperature_1D[i_local_min - 1];
476  pre[i][j][i_local_min] = (gasconstant * rho[i][j][i_local_min]) * temperature_1D[i_local_min ];
477  }
478 
479  //temperature in the obstacle higher than in the gas right of it
480  if ((energy_right > 0) && (i_local_max != (z - 1))) {
481  temperature_1D[i_local_max] = temperature_1D[i_local_max] + fabs(energy_right / (specific_heat_capacity_obstacle * mass[i_local_max]));
482  temperature_1D[i_local_max + 1] = temperature_1D[i_local_max + 1] - fabs(energy_right / (specific_heat_capacity_gas * mass[i_local_max + 1]));
483 
484  pre[i][j][i_local_max] = (gasconstant * rho[i][j][i_local_max]) * temperature_1D[i_local_max ];
485  pre[i][j][i_local_max + 1] = (gasconstant * rho[i][j][i_local_max + 1]) * temperature_1D[i_local_max + 1];
486  }
487 
488  //temperature in the obstacle lower than in the gas right of it
489  if ((energy_right <= 0) && (i_local_max != (z - 1))) {
490  temperature_1D[i_local_max] = temperature_1D[i_local_max] - fabs(energy_right / (specific_heat_capacity_obstacle * mass[i_local_max]));
491  temperature_1D[i_local_max + 1] = temperature_1D[i_local_max + 1] + fabs(energy_right / (specific_heat_capacity_gas * mass[i_local_max + 1]));
492 
493  pre[i][j][i_local_max] = (gasconstant * rho[i][j][i_local_max]) * temperature_1D[i_local_max ];
494  pre[i][j][i_local_max + 1] = (gasconstant * rho[i][j][i_local_max + 1]) * temperature_1D[i_local_max + 1];
495  }
496 
497  //to debug
498  //printf("energy_left: %g energy_right: %g dt: %g\n", energy_left, energy_right, dt);
499 
500  if (dom_counter > 2) {
501 
502  for (n = i_local_min + 1; n < i_local_max; n++) {
503  //here the heat conductivity constants are determined
504  temperature_1D_future[n] = alpha_solid_z * temperature_1D[n + 1] +
505  (1 - 2 * alpha_solid_z) * temperature_1D[n] +
506  alpha_solid_z * temperature_1D[n - 1];
507  }
508 
509  for (k = i_local_min + 1; k < i_local_max; k++) {
510  pre[i][j][k] = (gasconstant * rho[i][j][k]) * temperature_1D_future[k];
511  }
512 
513  for (k = i_local_min; k < i_local_max; k++) {
514  temperature_1D[k] = 0.0;
515  temperature_1D_future[k] = 0.0;
516  mass[k] = 0.0;
517  dom_1D[k] = 0;
518  }
519  }
520 
521  i_local_min = i_local_max + 1;
522 
523  }
524  }
525  }
526 
527  free_diffusion(&temperature_1D, &temperature_1D_future, &mass, &dom_1D);
528 
529 #ifdef _OPENMP
530  }
531 #endif
532 
533 
534  return 0;
535 }
Generated on Thu Oct 10 2013 17:15:52 for TYCHO by   doxygen 1.8.1.1