File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/calc_vort_out.f

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     
489