TYCHO  1.3.0
 All Data Structures Files Functions Variables Enumerations Enumerator
/home/kapf/tycho_docu/dt_calc_first.c
Go to the documentation of this file.
1 /*
2  * dt_calc_first.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 
21 inline double max_dt_first(double a, double b) {
22  double tmp;
23 
24  if (a < b) tmp = b;
25  if (a == b)tmp = b;
26  if (a > b) tmp = a;
27 
28  return tmp;
29 }
30 
34 inline double min_dt_first(double a, double b) {
35  double tmp;
36 
37  if (a < b) tmp = a;
38  if (a == b)tmp = a;
39  if (a > b) tmp = b;
40 
41  return tmp;
42 }
43 
50 int dt_calc_first(int x, int y, int z) {
51  int i, j, k;
52  double rdt1, tmpdy, tmpdz, tmp;
53  double xvel, yvel, zvel, svel;
54  double max_tmp1, max_tmp2, max_tmp3, min_tmp1, dtx;
55  double kin_viscosity;
56  double s_visc;
57  double temperature;
58 
59  rdt1 = 0.0;
60  svel = 0.0;
61  kin_viscosity = 0.0;
62  s_visc = 0.0;
63 
64 
65  /*
66  The signal velocity is calculated and compared to the
67  gas velocity in each cell. The routine searches for
68  the maximum in the whole computational domain.
69  */
70  if (dimension == 1) {
71  j = 0;
72  k = 0;
73 
74  for (i = 0; i < x; i++) {
75  if (dom[i][j][k] == 0) {
76  svel = sqrt(Gamma1 * pre[i][j][k] / rho[i][j][k]) / zdx[i];
77  xvel = fabs(vx[i][j][k]) / zdx[i];
78  tmp = max_dt_first(svel, xvel);
79  rdt1 = max_dt_first(tmp, rdt1);
80 
81  if (viscosity_on_off == 1) {
82  temperature = pre[i][j][k] / (gasconstant * rho[i][j][k]);
83  kin_viscosity = kinematic_viscosity(temperature, rho[i][j][k]);
84  s_visc = (kin_viscosity / (tmp * tmp));
85  rdt1 = max_dt_first(s_visc, rdt1);
86  }
87  }
88  }
89  }
90 
91 
92 #ifdef _OPENMP
93 #pragma omp parallel default(none) \
94  private(j, i, k, tmpdy, tmp, svel, xvel, \
95  yvel, max_tmp2, max_tmp3, temperature, \
96  kin_viscosity,s_visc) \
97  shared(x, y, z, dom, rho, pre, Gamma1, vx, \
98  vy, zdx, zdy, dimension, rdt1, viscosity_on_off, \
99  courant, gasconstant)
100  {
101 #endif
102  if (dimension == 2) {
103 
104  k = 0;
105 
106  for (i = 0; i < x; i++) {
107 
108 #ifdef _OPENMP
109 #pragma omp for schedule(static)
110 #endif
111 
112  for (j = 0; j < y; j++) {
113  if (dom[i][j][k] == 0) {
114  tmpdy = zdy[j];
115  tmp = min_dt_first(zdx[i], tmpdy);
116  svel = sqrt(Gamma1 * pre[i][j][k] / rho[i][j][k]) / zdx[i];
117  xvel = fabs(vx[i][j][k]) / zdx[i];
118  yvel = fabs(vy[i][j][k]) / tmpdy;
119  max_tmp2 = max_dt_first(svel, xvel);
120  max_tmp3 = max_dt_first(max_tmp2, yvel);
121  rdt1 = max_dt_first(max_tmp3, rdt1);
122 
123  if (viscosity_on_off == 1) {
124  temperature = pre[i][j][k] / (gasconstant * rho[i][j][k]);
125  kin_viscosity = kinematic_viscosity(temperature, rho[i][j][k]);
126  s_visc = (kin_viscosity / (tmp * tmp));
127  rdt1 = max_dt_first(s_visc, rdt1);
128  }
129  }
130  }
131  }
132 
133  }
134 
135 #ifdef _OPENMP
136  }
137 #endif
138 
139 #ifdef _OPENMP
140 #pragma omp parallel default(none) \
141  private(j, i, k, tmpdy, tmp, tmpdz, svel, xvel, \
142  yvel, zvel, max_tmp1, max_tmp2, max_tmp3, \
143  min_tmp1,temperature, viscosity_on_off, \
144  kin_viscosity, s_visc) \
145  shared(x, y, z, dom, rho, pre, Gamma1, vx, \
146  vy, vz, zdx, zdy, zdz,dimension, rdt1, gasconstant)
147  {
148 #endif
149 
150  if (dimension == 3) {
151  for (i = 0; i < x; i++) {
152 
153 #ifdef _OPENMP
154 #pragma omp for schedule(static)
155 #endif
156 
157  for (j = 0; j < y; j++) {
158  for (k = 0; k < z; k++) {
159  if (dom[i][j][k] == 0) {
160  tmpdy = zdy[j];
161  tmpdz = zdz[k];
162  min_tmp1 = min_dt_first(tmpdy, tmpdz);
163  tmp = min_dt_first(zdx[i], min_tmp1);
164  svel = sqrt(Gamma1 * pre[i][j][k] / rho[i][j][k]) / zdx[i];
165  xvel = fabs(vx[i][j][k]) / zdx[i];
166  yvel = fabs(vy[i][j][k]) / tmpdy;
167  zvel = fabs(vz[i][j][k]) / tmpdz;
168  max_tmp1 = max_dt_first(xvel, yvel);
169  max_tmp2 = max_dt_first(zvel, svel);
170  max_tmp3 = max_dt_first(max_tmp1, max_tmp2);
171  rdt1 = max_dt_first(max_tmp3, rdt1);
172 
173  if (viscosity_on_off == 1) {
174  temperature = pre[i][j][k] / (gasconstant * rho[i][j][k]);
175  kin_viscosity = kinematic_viscosity(temperature, rho[i][j][k]);
176  s_visc = (kin_viscosity / (tmp * tmp));
177  rdt1 = max_dt_first(s_visc, rdt1);
178  }
179 
180  }
181  }
182  }
183  }
184  }
185 
186 #ifdef _OPENMP
187  }
188 #endif
189  // here the CFL condition is applied
190  dtx = courant / rdt1;
191  dt = dtx;
192 
193  intial_soundspeed = rdt1*spacing;
194 
195  if (dt == 0.0) {
196  printf("A time_simstep of zero ---> STOP\n");
197  exit(42);
198  }
199 
200  return 0;
201 }