MFIX  2016-1
calc_vort_out.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: CALC_VORTICITY C
4 ! Purpose: Computes the vorticity C
5 ! C
6 ! Author: Jeff Dietiker Date: 21-Feb-08 C
7 ! Reviewer: Date: C
8 ! C
9 ! Revision Number # Date: ##-###-## C
10 ! Author: # C
11 ! Purpose: # C
12 ! C
13 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14  SUBROUTINE calc_vorticity
15 
16  USE param
17  USE param1
18  USE parallel
19  USE constant
20  USE run
21  USE toleranc
22  USE geometry
23  USE indices
24  USE compar
25  USE sendrecv
26  USE fldvar
27  USE quadric
28  USE cutcell
29  USE functions
30 
31  IMPLICIT NONE
32  INTEGER :: I,J,K,IJK,IP,IM,JP,JM,KP,KM,IJKE,IJKN,IJKT,IJKW,IJKS,IJKB
33  DOUBLE PRECISION :: DU_DX,DU_DY,DU_DZ,DV_DX,DV_DY,DV_DZ,DW_DX,DW_DY,DW_DZ
34  DOUBLE PRECISION :: OMEGA_X,OMEGA_Y,OMEGA_Z
35  DOUBLE PRECISION :: LAMBDA_1,LAMBDA_2,LAMBDA_3
36  DOUBLE PRECISION,DIMENSION(3,3) :: OMEGA,SS,AA
37  DOUBLE PRECISION,DIMENSION(4) :: POLY
38 
39  DO ijk = ijkstart3, ijkend3
40 
41  i = i_of(ijk)
42  j = j_of(ijk)
43  k = k_of(ijk)
44 
45  im = i - 1
46  ip = i + 1
47 
48  jm = j - 1
49  jp = j + 1
50 
51  km = k - 1
52  kp = k + 1
53 
54 
55  ijke = funijk(ip,j,k)
56  ijkn = funijk(i,jp,k)
57  ijkt = funijk(i,j,kp)
58 
59  ijkw = funijk(im,j,k)
60  ijks = funijk(i,jm,k)
61  ijkb = funijk(i,j,km)
62 
63 
64 
65  IF (fluid_at(ijk).AND.interior_cell_at(ijk)) THEN
66 
67  du_dx = zero
68  du_dy = zero
69  du_dz = zero
70 
71  dv_dx = zero
72  dv_dy = zero
73  dv_dz = zero
74 
75  dw_dx = zero
76  dw_dy = zero
77  dw_dz = zero
78 
79 !======================================================================
80 ! Get Velocity derivatives
81 !======================================================================
82 
83 !======================================================================
84 ! du/dx
85 !======================================================================
86 
87  IF ((fluid_at(ijke)).AND.(fluid_at(ijkw))) THEN
88 
89  du_dx = (u_g(ijke) - u_g(ijkw)) / (x_u(ijke) - x_u(ijkw))
90 
91  ELSE IF ((fluid_at(ijke)).AND.(.NOT.fluid_at(ijkw))) THEN
92 
93  du_dx = (u_g(ijke) - u_g(ijk)) / (x_u(ijke) - x_u(ijk))
94 
95  ELSE IF ((fluid_at(ijkw)).AND.(.NOT.fluid_at(ijke))) THEN
96 
97  du_dx = (u_g(ijk) - u_g(ijkw)) / (x_u(ijk) - x_u(ijkw))
98 
99  END IF
100 
101 !======================================================================
102 ! du/dy
103 !======================================================================
104 
105  IF ((fluid_at(ijkn)).AND.(fluid_at(ijks))) THEN
106 
107  du_dy = (u_g(ijkn) - u_g(ijks)) / (y_u(ijkn) - y_u(ijks))
108 
109  ELSE IF ((fluid_at(ijkn)).AND.(.NOT.fluid_at(ijks))) THEN
110 
111  du_dy = (u_g(ijkn) - u_g(ijk)) / (y_u(ijkn) - y_u(ijk))
112 
113  ELSE IF ((fluid_at(ijks)).AND.(.NOT.fluid_at(ijkn))) THEN
114 
115  du_dy = (u_g(ijk) - u_g(ijks)) / (y_u(ijk) - y_u(ijks))
116 
117  END IF
118 
119 !======================================================================
120 ! du/dz
121 !======================================================================
122 
123  IF(do_k) THEN
124  IF ((fluid_at(ijkt)).AND.(fluid_at(ijkb))) THEN
125 
126  du_dz = (u_g(ijkt) - u_g(ijkb)) / (z_u(ijkt) - z_u(ijkb))
127 
128  ELSE IF ((fluid_at(ijkt)).AND.(.NOT.fluid_at(ijkb))) THEN
129 
130  du_dz = (u_g(ijkt) - u_g(ijk)) / (z_u(ijkt) - z_u(ijk))
131 
132  ELSE IF ((fluid_at(ijkb)).AND.(.NOT.fluid_at(ijkt))) THEN
133 
134  du_dz = (u_g(ijk) - u_g(ijkb)) / (z_u(ijk) - z_u(ijkb))
135 
136  END IF
137  ENDIF
138 
139 !======================================================================
140 ! dv/dx
141 !======================================================================
142 
143  IF ((fluid_at(ijke)).AND.(fluid_at(ijkw))) THEN
144 
145  dv_dx = (v_g(ijke) - v_g(ijkw)) / (x_v(ijke) - x_v(ijkw))
146 
147  ELSE IF ((fluid_at(ijke)).AND.(.NOT.fluid_at(ijkw))) THEN
148 
149  dv_dx = (v_g(ijke) - v_g(ijk)) / (x_v(ijke) - x_v(ijk))
150 
151  ELSE IF ((fluid_at(ijkw)).AND.(.NOT.fluid_at(ijke))) THEN
152 
153  dv_dx = (v_g(ijk) - v_g(ijkw)) / (x_v(ijk) - x_v(ijkw))
154 
155  END IF
156 
157 !======================================================================
158 ! dv/dy
159 !======================================================================
160 
161  IF ((fluid_at(ijkn)).AND.(fluid_at(ijks))) THEN
162 
163  dv_dy = (v_g(ijkn) - v_g(ijks)) / (y_v(ijkn) - y_v(ijks))
164 
165  ELSE IF ((fluid_at(ijkn)).AND.(.NOT.fluid_at(ijks))) THEN
166 
167  dv_dy = (v_g(ijkn) - v_g(ijk)) / (y_v(ijkn) - y_v(ijk))
168 
169  ELSE IF ((fluid_at(ijks)).AND.(.NOT.fluid_at(ijkn))) THEN
170 
171  dv_dy = (v_g(ijk) - v_g(ijks)) / (y_v(ijk) - y_v(ijks))
172 
173  END IF
174 
175 !======================================================================
176 ! dv/dz
177 !======================================================================
178 
179  IF(do_k) THEN
180  IF ((fluid_at(ijkt)).AND.(fluid_at(ijkb))) THEN
181 
182  dv_dz = (v_g(ijkt) - v_g(ijkb)) / (z_v(ijkt) - z_v(ijkb))
183 
184  ELSE IF ((fluid_at(ijkt)).AND.(.NOT.fluid_at(ijkb))) THEN
185 
186  dv_dz = (v_g(ijkt) - v_g(ijk)) / (z_v(ijkt) - z_v(ijk))
187 
188  ELSE IF ((fluid_at(ijkb)).AND.(.NOT.fluid_at(ijkt))) THEN
189 
190  dv_dz = (v_g(ijk) - v_g(ijkb)) / (z_v(ijk) - z_v(ijkb))
191 
192  END IF
193  ENDIF
194 
195 !======================================================================
196 ! dw/dx
197 !======================================================================
198 
199  IF(do_k) THEN
200  IF ((fluid_at(ijke)).AND.(fluid_at(ijkw))) THEN
201 
202  dw_dx = (w_g(ijke) - w_g(ijkw)) / (x_w(ijke) - x_w(ijkw))
203 
204  ELSE IF ((fluid_at(ijke)).AND.(.NOT.fluid_at(ijkw))) THEN
205 
206  dw_dx = (w_g(ijke) - w_g(ijk)) / (x_w(ijke) - x_w(ijk))
207 
208  ELSE IF ((fluid_at(ijkw)).AND.(.NOT.fluid_at(ijke))) THEN
209 
210  dw_dx = (w_g(ijk) - w_g(ijkw)) / (x_w(ijk) - x_w(ijkw))
211 
212  END IF
213 
214 !======================================================================
215 ! dw/dy
216 !======================================================================
217 
218  IF ((fluid_at(ijkn)).AND.(fluid_at(ijks))) THEN
219 
220  dw_dy = (w_g(ijkn) - w_g(ijks)) / (y_w(ijkn) - y_w(ijks))
221 
222  ELSE IF ((fluid_at(ijkn)).AND.(.NOT.fluid_at(ijks))) THEN
223 
224  dw_dy = (w_g(ijkn) - w_g(ijk)) / (y_w(ijkn) - y_w(ijk))
225 
226  ELSE IF ((fluid_at(ijks)).AND.(.NOT.fluid_at(ijkn))) THEN
227 
228  dw_dy = (w_g(ijk) - w_g(ijks)) / (y_w(ijk) - y_w(ijks))
229 
230  END IF
231 
232 !======================================================================
233 ! dw/dz
234 !======================================================================
235 
236  IF ((fluid_at(ijkt)).AND.(fluid_at(ijkb))) THEN
237 
238  dw_dz = (w_g(ijkt) - w_g(ijkb)) / (z_w(ijkt) - z_w(ijkb))
239 
240  ELSE IF ((fluid_at(ijkt)).AND.(.NOT.fluid_at(ijkb))) THEN
241 
242  dw_dz = (w_g(ijkt) - w_g(ijk)) / (z_w(ijkt) - z_w(ijk))
243 
244  ELSE IF ((fluid_at(ijkb)).AND.(.NOT.fluid_at(ijkt))) THEN
245 
246  dw_dz = (w_g(ijk) - w_g(ijkb)) / (z_w(ijk) - z_w(ijkb))
247 
248  END IF
249  ENDIF
250 
251 !======================================================================
252 ! Build Vorticity components and magnitude
253 !======================================================================
254 
255  omega_x = dw_dy - dv_dz
256  omega_y = du_dz - dw_dx
257  omega_z = dv_dx - du_dy
258 
259 
260  vorticity(ijk) = dsqrt(omega_x**2 + omega_y**2 + omega_z**2)
261 
262 !======================================================================
263 ! Build Matrices OMEGA , SS, and AA
264 !
265 ! OMEGA_ij = 1/2 *(dui/dxj + duj/dxi)
266 ! SS_ij = 1/2 *(dui/dxj - duj/dxi)
267 ! AA = OMEGA^2 + SS^2
268 !======================================================================
269 
270  omega(1,1) = du_dx
271  omega(1,2) = half * (du_dy + dv_dx)
272  omega(1,3) = half * (du_dz + dw_dx)
273 
274  omega(2,1) = omega(1,2)
275  omega(2,2) = dv_dy
276  omega(2,3) = half * (dv_dz + dw_dy)
277 
278  omega(3,1) = omega(1,3)
279  omega(3,2) = omega(2,3)
280  omega(3,3) = dw_dz
281 
282 
283  ss(1,1) = zero
284  ss(1,2) = half * (du_dy - dv_dx)
285  ss(1,3) = half * (du_dz - dw_dx)
286 
287  ss(2,1) = - ss(1,2)
288  ss(2,2) = zero
289  ss(2,3) = half * (dv_dz - dw_dy)
290 
291  ss(3,1) = - ss(1,3)
292  ss(3,2) = - ss(2,3)
293  ss(3,3) = zero
294 
295 
296  aa = matmul(omega,omega) + matmul(ss,ss)
297 
298 !======================================================================
299 ! Build Characteristic polynomial of AA
300 !======================================================================
301 
302  poly(1) = - one
303  poly(2) = aa(1,1) + aa(2,2) + aa(3,3)
304  poly(3) = aa(2,1)*aa(1,2) + aa(3,1)*aa(1,3) + aa(3,2)*aa(2,3) &
305  -( aa(1,1)*aa(2,2) + aa(1,1)*aa(3,3) + aa(2,2)*aa(3,3) )
306  poly(4) = aa(1,1)*aa(2,2)*aa(3,3) + aa(1,2)*aa(2,3)*aa(3,1) + aa(2,1)*aa(3,2)*aa(1,3) &
307  -( aa(3,1)*aa(2,2)*aa(1,3) + aa(1,2)*aa(2,1)*aa(3,3) + aa(2,3)*aa(3,2)*aa(1,1) )
308 
309 
310 
311 !======================================================================
312 ! Find roots of characteristic polynomial = eigenvalues
313 ! using Bairstow method
314 !======================================================================
315 
316  CALL bairstow(poly,lambda_1,lambda_2,lambda_3)
317 
318  lambda2(ijk) = lambda_2
319 
320  ENDIF
321 
322  END DO
323 
324  RETURN
325 
326 
327  END SUBROUTINE calc_vorticity
328 
329 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
330 ! C
331 ! Module name: BAIRSTOW C
332 ! Purpose: FIND THE ROOTS OF A POLYNOMIAL USING BAIRSTOW METHOD C
333 ! LIMITED TO POLYNOMIALS OF DEGREE 3 C
334 ! (USED TO FIND THE EIGENVALUES OF A 3x3 MATRIX) C
335 ! C
336 ! Author: Jeff Dietiker Date: 17-Jul-08 C
337 ! Reviewer: Date: C
338 ! C
339 ! Revision Number # Date: ##-###-## C
340 ! Author: # C
341 ! Purpose: # C
342 ! C
343 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
344  SUBROUTINE bairstow(A,X1,X2,X3)
345  USE param1
346 
347  IMPLICIT NONE
348  INTEGER :: N,NMAX,ITERMAX,J,ITER
349 
350  parameter(nmax = 4) ! MAXIMUM NUMBER OF POLYNOMIAL COEFFICIENTS ALLOWED
351 
352  DOUBLE PRECISION, DIMENSION(NMAX):: A,B,C
353  DOUBLE PRECISION:: R,S,DELTA,DELTAR,DELTAS,DENOM
354  DOUBLE PRECISION:: EPS
355  DOUBLE PRECISION:: X1,X2,X3
356  DOUBLE PRECISION:: BUFFER1
357 
358  LOGICAL :: TEST1, TEST2, TEST3
359 
360  n = 3 ! This subroutine is designed for a ploynomial of degree three only
361  itermax = 200 ! It is not anticipated that ITERMAX nor EPS would need to
362  eps = 1.0d-6 ! be modified
363 
364 !=======================================================================
365 ! POLYNOMIAL REDUCTION
366 ! ALWAYS START WITH INITIAL GUESS R = - 0.9 AND S = -1.0
367 ! AND ITERATE
368 !=======================================================================
369 
370  r = -0.9d0
371  s = -1.0d0
372 
373  iter = 0
374 
375 !=======================================================================
376 ! COMPUTING B'S ANS C'S
377 !=======================================================================
378 
379 20 CONTINUE
380 
381  b(1) = a(1)
382  b(2) = a(2) + r * b(1)
383  c(1) = b(1)
384  c(2) = b(2) + r * c(1)
385  b(3) = a(3) + r * b(2) + s * b(1)
386  c(3) = b(3) + r * c(2) + s * c(1)
387  b(4) = a(4) + r * b(3) + s * b(2)
388 
389 !=======================================================================
390 ! UPDATING VALUES OF R AND S
391 !=======================================================================
392 
393  denom = c(n-1) * c(n-1) - c(n)*c(n-2)
394 
395  deltar = (-b(n)*c(n-1) + b(n+1)*c(n-2) ) / denom
396  deltas = (-c(n-1)*b(n+1) + b(n)*c(n) ) / denom
397 
398  r = r + deltar
399  s = s + deltas
400 
401  iter = iter +1
402 
403 !=======================================================================
404 ! CHECKING CONVERGENCE
405 !=======================================================================
406 
407  test1 = (abs( b(n) ) > eps)
408  test2 = (abs( b(n+1)) > eps)
409  test3 = (iter <= itermax)
410 
411  IF(test1.AND.test2.AND.test3) GOTO 20
412 
413  IF (.NOT.test3) THEN
414 ! WRITE(*,*)'ERROR IN BAIRSTOW SUBROUTINE:'
415 ! WRITE(*,*)'DIVERGENCE... THE PROGRAM WILL BE TERMINATED NOW.'
416 ! CALL MFIX_EXIT(myPE)
417  x1 = undefined
418  x2 = undefined
419  x3 = undefined
420  ENDIF
421 
422 !=======================================================================
423 ! NOW, THE ITERATIVE SCHEME HAS CONVERGED
424 ! THE QUADRATIC FACTOR IS X*X-R*X-S
425 ! SOLVING AND DISPLAYING THE TWO ROOTS (LIMITED TO REAL ROOTS)
426 !=======================================================================
427 
428  delta = r*r + 4.0d0 * s
429 
430  IF (delta.LT.0.0d0) THEN
431 
432 ! WRITE(*,*) ' ERROR IN BAIRSTOW SUBROUTINE:TWO COMPLEX ROOTS:'
433 
434  x1 = undefined
435  x2 = undefined
436 
437  ELSE
438 
439  x1 = ( r + dsqrt(delta)) / ( 2.0d0 )
440  x2 = ( r - dsqrt(delta)) / ( 2.0d0 )
441 
442  ENDIF
443 
444 !=======================================================================
445 ! IN PREPARATION FOR THE LAST LINEAR FACTOR,
446 ! SET N = N-2
447 ! RENAME THE B'S AS A'S
448 !=======================================================================
449 
450  n = n - 2
451 
452  DO j=1,n+1
453  a(j) = b(j)
454  END DO
455 
456  IF(dabs(a(1))<1.0d-9) THEN
457  x3 = undefined
458  ELSE
459  x3 = -a(2)/a(1)
460  ENDIF
461 
462 !=======================================================================
463 ! SORING OUT ROOTS IN ASCENDING ORDER
464 !=======================================================================
465 
466  IF(x1>x2) THEN
467  buffer1 = x1
468  x1 = x2
469  x2 = buffer1
470  ENDIF
471 
472  IF(x2>x3) THEN
473  buffer1 = x2
474  x2 = x3
475  x3 = buffer1
476  ENDIF
477 
478  IF(x1>x2) THEN
479  buffer1 = x1
480  x1 = x2
481  x2 = buffer1
482  ENDIF
483 
484  RETURN
485 
486  END SUBROUTINE bairstow
487 
488 
cutcell::y_v
double precision, dimension(:), allocatable y_v
Definition: cutcell_mod.f:55
cutcell::z_u
double precision, dimension(:), allocatable z_u
Definition: cutcell_mod.f:51
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
param1::one
double precision, parameter one
Definition: param1_mod.f:29
cutcell::x_u
double precision, dimension(:), allocatable x_u
Definition: cutcell_mod.f:49
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
cutcell::y_w
double precision, dimension(:), allocatable y_w
Definition: cutcell_mod.f:60
cutcell::z_v
double precision, dimension(:), allocatable z_v
Definition: cutcell_mod.f:56
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
cutcell::x_w
double precision, dimension(:), allocatable x_w
Definition: cutcell_mod.f:59
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
bairstow
subroutine bairstow(A, X1, X2, X3)
Definition: calc_vort_out.f:345
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
cutcell::x_v
double precision, dimension(:), allocatable x_v
Definition: cutcell_mod.f:54
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
param
Definition: param_mod.f:2
cutcell::vorticity
double precision, dimension(:), allocatable vorticity
Definition: cutcell_mod.f:406
cutcell::interior_cell_at
logical, dimension(:), allocatable interior_cell_at
Definition: cutcell_mod.f:40
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
calc_vorticity
subroutine calc_vorticity
Definition: calc_vort_out.f:15
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
cutcell::lambda2
double precision, dimension(:), allocatable lambda2
Definition: cutcell_mod.f:407
cutcell::z_w
double precision, dimension(:), allocatable z_w
Definition: cutcell_mod.f:61
cutcell::y_u
double precision, dimension(:), allocatable y_u
Definition: cutcell_mod.f:50
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27