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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: SET_GEOMETRY                                           !
4     !  Author: M. Syamlal                                 Date: 21-JAN-92  !
5     !                                                                      !
6     !  Purpose: Calculate X, X_E,  oX, oX_E                                !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE SET_GEOMETRY
10     
11     ! Global Variables:
12     !---------------------------------------------------------------------//
13     ! Domain decomposition and dimensions
14           use geometry, only: DX, XLENGTH, oDX, oDX_E
15           use geometry, only: DY, YLENGTH, oDZ, oDZ_T
16           use geometry, only: DZ, ZLENGTH, oDY, oDY_N
17           use geometry, only: X, X_E, oX, oX_E, XMIN
18           use geometry, only: Z, Z_T
19     ! Domain indices.
20           use geometry, only: IJKMAX3
21           use geometry, only: DO_I, IMIN1, IMAX1, IMAX2, IMAX3, IMIN3, IMAX
22           use geometry, only: DO_J, JMIN1, JMAX1, JMAX2, JMAX3, JMIN3
23           use geometry, only: DO_K, KMIN1, KMAX1, KMAX2, KMAX3, KMIN3
24     ! Averaging factors.
25           use geometry, only: FX_E, FX_E_bar, FX, FX_bar
26           use geometry, only: FY_N, FY_N_bar
27           use geometry, only: FZ_T, FZ_T_bar
28     ! Domain flags.
29           use geometry, only: ICBC_FLAG
30           use geometry, only: FLAG, FLAG3
31           use geometry, only: FLAG_E, FLAG_N, FLAG_T
32     ! Domain volumes and areas.
33           use geometry, only: VOL, AYZ, AXZ, AXY          ! Scalar grid
34           use geometry, only: VOL_U, AYZ_U, AXZ_U, AXY_U  ! X-Momentum
35           use geometry, only: VOL_V, AYZ_V, AXZ_V, AXY_V  ! Y-Momentum
36           use geometry, only: VOL_W, AYZ_W, AXZ_W, AXY_W  ! Z-Momentum
37     ! Cyclic domain flags.
38           use geometry, only: CYCLIC
39           use geometry, only: CYCLIC_X, CYCLIC_X_PD, CYCLIC_X_MF
40           use geometry, only: CYCLIC_Y, CYCLIC_Y_PD, CYCLIC_Y_MF
41           use geometry, only: CYCLIC_Z, CYCLIC_Z_PD, CYCLIC_Z_MF
42     ! Flag for cylindrical coordinates.
43           use geometry, only: CYLINDRICAL
44     ! For cylindrical_2d simulations
45           use geometry, only: CYLINDRICAL_2D, I_CYL_NUM, I_CYL_TRANSITION
46           use geometry, only: cyl_X, cyl_X_E
47     ! Axis decomposition
48           USE param, only: DIMENSION_I, DIMENSION_J, DIMENSION_K
49           USE param, only: DIMENSION_3, DIMENSION_4
50           USE param, only: DIMENSION_3G, DIMENSION_3L, DIMENSION_3P
51     ! MPI-Domain decompoint and rank flags.
52           use compar, only: NODESI, NODESJ, NODESK, myPE
53     ! Rank specific decompositions.
54           use compar, only: IJKSIZE3_ALL
55           USE compar, only: iStart3, iEnd3, iStart4, iEnd4
56           USE compar, only: jStart3, jEnd3, jStart4, jEnd4
57           USE compar, only: kStart3, kEnd3, kStart4, kEnd4
58     ! Flag for specificed constant mass flux.
59           use bc, only: Flux_g
60     ! Flag for POST_MFIX
61           use cdist, only: bDoing_postmfix
62     
63     ! Global Parameters:
64     !---------------------------------------------------------------------//
65           use param1, only: ZERO, HALF, ONE, UNDEFINED
66     
67     ! Module procedures
68     !---------------------------------------------------------------------//
69           use mpi_utility, only: GLOBAL_ALL_SUM
70           use error_manager
71           use toleranc
72     
73           IMPLICIT NONE
74     
75     ! Local Variables:
76     !---------------------------------------------------------------------//
77     ! Generic loop indices
78           INTEGER :: I, J, K
79     ! X-direction dimension of U-momentum cell
80           DOUBLE PRECISION :: DX_E
81     ! Y-direction dimension of V-momentum cell
82           DOUBLE PRECISION :: DY_N
83     ! Z-direction dimension of W-momentum cell
84           DOUBLE PRECISION :: DZ_T
85     ! Error Flag
86           INTEGER :: IER
87     ! Local variables for cylindrical_2d simulation
88           integer i_cyl_min, i_cyl_max
89           double precision l_ver, l_ab, rrr, ddy
90     !......................................................................!
91     
92     ! Initialize the error manager.
93           CALL INIT_ERR_MSG("SET_GEOMETRY")
94     
95     ! Allocate geometry arrays.
96           CALL ALLOCATE_ARRAYS_GEOMETRY
97     
98     !  Determine the cyclic direction with a specified mass flux
99           CYCLIC_X_MF = (FLUX_G /= UNDEFINED .AND. CYCLIC_X_PD)
100           CYCLIC_Y_MF = (FLUX_G /= UNDEFINED .AND. CYCLIC_Y_PD)
101           CYCLIC_Z_MF = (FLUX_G /= UNDEFINED .AND. CYCLIC_Z_PD)
102     
103     ! Force the cyclic flag if cyclic with pressure drop.
104           IF (CYCLIC_X_PD) CYCLIC_X = .TRUE.
105           IF (CYCLIC_Y_PD) CYCLIC_Y = .TRUE.
106           IF (CYCLIC_Z_PD) CYCLIC_Z = .TRUE.
107           CYCLIC = CYCLIC_X .OR. CYCLIC_Y .OR. CYCLIC_Z
108     
109           IF(CYLINDRICAL .AND. COMPARE(ZLENGTH,8.D0*ATAN(ONE)) .AND. DO_K) &
110              CYCLIC_Z = .TRUE.
111     
112           IF(CYCLIC_X) THEN
113              DX(1) = DX(IMAX1)
114              DX(IMAX2) = DX(IMIN1)
115              IF(NODESI.NE.1) THEN
116                 DX(IMIN3) = DX(IMAX1-1)
117                 DX(IMAX3) = DX(IMIN1+1)
118              ENDIF
119           ENDIF
120     
121           IF(CYCLIC_Y) THEN
122              DY(1) = DY(JMAX1)
123              DY(JMAX2) = DY(JMIN1)
124              IF(NODESJ.NE.1) THEN
125                 DY(JMIN3) = DY(JMAX1-1)
126                 DY(JMAX3) = DY(JMIN1+1)
127              ENDIF
128           ENDIF
129     
130           IF (CYCLIC_Z) THEN
131              DZ(1) = DZ(KMAX1)
132              DZ(KMAX2) = DZ(KMIN1)
133              IF(NODESK.NE.1) THEN
134                 DZ(KMIN3) = DZ(KMAX1-1)
135                 DZ(KMAX3) = DZ(KMIN1+1)
136              ENDIF
137           ENDIF
138     
139     
140     ! Initialize the X-axis variables for CYLINDRICAL coordinates.
141           IF(CYLINDRICAL) THEN
142              X(IMIN3:IMAX3)    = UNDEFINED
143              X_E(IMIN3:IMAX3)  = UNDEFINED
144              OX(IMIN3:IMAX3)   = UNDEFINED
145              OX_E(IMIN3:IMAX3) = UNDEFINED
146              ODX(IMIN3:IMAX3)  = UNDEFINED
147     
148              IF(XMIN == ZERO) THEN
149                 ODX(1) = ONE/DX(1)
150                 OX(1) = UNDEFINED
151                 OX_E(1) = UNDEFINED
152                 IF (DO_I) THEN
153                    X(1) = -HALF*DX(1)
154                    X_E(1) = 0.0
155                 ELSE
156                    X(1) = HALF*DX(1)
157                    X_E(1) = DX(1)
158                 ENDIF
159              ELSE
160                 IF (DO_I) THEN
161                    X_E(1) = XMIN
162                    X(1) = XMIN - HALF*DX(1)
163                 ELSE
164                    X_E(1) = XMIN + DX(1)
165                    X(1) = XMIN + HALF*DX(1)
166                 ENDIF
167                 OX(1) = ONE/X(1)
168                 OX_E(1) = ONE/X_E(1)
169                 ODX(1) = ONE/DX(1)
170              ENDIF
171     
172              IF (DO_I) THEN
173                 DO I = IMIN1, IMAX2
174                    X(I) = X(I-1) + (DX(I-1)+DX(I))/2.
175                    X_E(I) = X_E(I-1) + DX(I)
176                    OX(I) = ONE/X(I)
177                    OX_E(I) = ONE/X_E(I)
178                    ODX(I) = ONE/DX(I)
179                 END DO
180              ENDIF
181     
182     ! Initialize the X-axis variables for CARTESIAN coordinates.
183           ELSE
184     
185              X(IMIN3:IMAX3) = ONE
186              X_E(IMIN3:IMAX3) = ONE
187              OX(IMIN3:IMAX3) = ONE
188              OX_E(IMIN3:IMAX3) = ONE
189              ODX(IMIN3:IMAX3) = ONE/DX(IMIN3:IMAX3)
190     
191     !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++<<
192     !     Implementation of the 2.5D model by Li et al. CES 123 (2015) 236-246.
193     !     A computational domain of two wedges connected by a thin plate is used.
194     !     The model is invoked by setting CYLINDRICAL_2D to .TRUE.
195     !     The half width of the plate is determined by I_CYL_NUM
196     !     Please refer to the paper for details of the model
197     !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++<<
198              IF(CYLINDRICAL_2D)THEN
199     
200                 IF (XMIN == ZERO) THEN
201                    IF (DO_I) THEN
202                       cyl_X(1) = -HALF*DX(1) -half*xlength
203                       cyl_X_E(1) = 0.D0 -half*xlength
204                    ELSE
205                       cyl_X(1) = HALF*DX(1) -half*xlength
206                       cyl_X_E(1) = DX(1) -half*xlength
207                    ENDIF
208                 ELSE
209                    IF (DO_I) THEN
210                       cyl_X_E(1) = XMIN -half*xlength
211                       cyl_X(1) = XMIN - HALF*DX(1) -half*xlength
212                    ELSE
213                       cyl_X_E(1) = XMIN + DX(1) -half*xlength
214                       cyl_X(1) = XMIN + HALF*DX(1) -half*xlength
215                    ENDIF
216                 ENDIF
217     
218                 IF (DO_I) THEN
219                    DO I = IMIN1, IMAX2
220                       cyl_X(I) = cyl_X(I-1) + (DX(I-1)+DX(I))/2.
221                       cyl_X_E(I) = cyl_X_E(I-1) + DX(I)
222                    END DO
223     
224                    DO I = IMIN3, IMAX3
225                       cyl_X(I) = cyl_X(I) *2.d0*sin(dz(1)/2.d0)
226                       cyl_X_E(I) = cyl_X_E(I) *2.d0*sin(dz(1)/2.d0)
227                       cyl_x(i) = abs(cyl_x(i))
228                       cyl_x_e(i) = abs(cyl_x_e(i))
229                    END DO
230     
231                    if(mod(imax,2).eq.1)then     ! odd
232                       i_cyl_min = (imax+1)/2 + 1 - i_cyl_num
233                       i_cyl_max = (imax+1)/2 + 1 + i_cyl_num
234     
235                       do i=i_cyl_min, i_cyl_max
236                          cyl_x(i)=dx(i)*(i_cyl_num + half) *2.d0*sin(dz(1)/2.d0)
237                       enddo
238                       do i=i_cyl_min ,i_cyl_max -1
239                          cyl_x_e(i)=dx(i)*(i_cyl_num + half) *2.d0*sin(dz(1)/2.d0)
240                       enddo
241                    else                    ! even
242                       i_cyl_min = (imax)/2 + 2 - i_cyl_num
243                       i_cyl_max = (imax)/2 + 1 + i_cyl_num
244     
245                       do i=i_cyl_min, i_cyl_max
246                          cyl_x(i)=dx(i)* i_cyl_num *2.d0*sin(dz(1)/2.d0)
247                       enddo
248                       do i=i_cyl_min, i_cyl_max -1
249                          cyl_x_e(i)=dx(i)* i_cyl_num *2.d0*sin(dz(1)/2.d0)
250                       enddo
251                    endif
252     
253     !  To smooth the transition from cylindrical to 2D domain
254     !  Only one or two cells transition are implemented.
255     
256                    if(i_cyl_transition .eq. 1)then
257                       l_ver=dx(i_cyl_min)*tan(dz(1)/2.0)
258                       l_ab = sqrt(l_ver*l_ver + (dx(i_cyl_min-1)+dx(i_cyl_min))**2)
259                       rrr = half*l_ab/sin(dz(1)/4.0)
260                       ddy = rrr - sqrt(rrr**2 - dx(i_cyl_min)**2)
261     
262                       cyl_x_e(i_cyl_min - 1) = cyl_x_e(i_cyl_min) + 2*ddy
263                       cyl_x_e(i_cyl_max ) = cyl_x_e(i_cyl_max - 1) + 2*ddy
264                    elseif(i_cyl_transition .eq. 2) then
265                       l_ver=2.d0*dx(i_cyl_min)*tan(dz(1)/2.0)
266                       l_ab = sqrt(l_ver*l_ver + (4.d0*dx(i_cyl_min))**2)
267                       rrr = half*l_ab/sin(dz(1)/4.0)
268                       ddy = rrr - sqrt(rrr**2 - dx(i_cyl_min)**2)
269     
270                       cyl_x_e(i_cyl_min ) = cyl_x_e(i_cyl_min + 1) + 2*ddy
271                       cyl_x_e(i_cyl_max - 1) = cyl_x_e(i_cyl_max - 2) + 2*ddy
272     
273                       ddy = rrr - sqrt(rrr**2 - (2.d0*dx(i_cyl_min))**2)
274                       cyl_x_e(i_cyl_min - 1 ) = cyl_x_e(i_cyl_min + 1) + 2*ddy
275                       cyl_x_e(i_cyl_max ) = cyl_x_e(i_cyl_max - 2) + 2*ddy
276     
277                       ddy = rrr - sqrt(rrr**2 - (3.d0*dx(i_cyl_min))**2)
278                       cyl_x_e(i_cyl_min - 2 ) = cyl_x_e(i_cyl_min + 1) + 2*ddy
279                       cyl_x_e(i_cyl_max + 1) = cyl_x_e(i_cyl_max - 2) + 2*ddy
280                    else ! no smooth
281     
282                    endif
283     
284                    do i=i_cyl_min - i_cyl_transition, i_cyl_max + i_cyl_transition
285                       cyl_x(i) = half * (cyl_x_e(i-1) + cyl_x_e(i))
286                    enddo
287     
288     !  To make sure all locations are positive
289                    do i= imin3,imax3
290                       cyl_x(i)=abs(cyl_x(i))
291                       cyl_x_e(i)=abs(cyl_x_e(i))
292                    enddo
293                 ENDIF  !IF (DO_I)
294              ENDIF  !IF(CYLINDRICAL_2D)THEN
295     
296           ENDIF
297     
298     
299     ! Initialize the Y-Axis variables.
300           ODY(JMIN3:JMAX3) = ONE/DY(JMIN3:JMAX3)
301     
302     
303     ! Initialize the Z-Axis variables.
304           DO K = 1, KMAX3
305              IF (K == 1) THEN
306                 Z(K) = ZERO - HALF*DZ(K)
307                 Z_T(K) = ZERO
308                 IF(NODESK.NE.1) THEN
309                    Z(K-1) =Z_T(K) - HALF*DZ(K-1)
310                    Z_T(K-1) = Z_T(K) - DZ(K-1)
311                 ENDIF
312                   ELSE
313                 Z(K) = Z_T(K-1) + HALF*DZ(K)
314                 Z_T(K) = Z_T(K-1) + DZ(K)
315              ENDIF
316              ODZ(K) = ONE/DZ(K)
317     
318              IF(NODESK.NE.1 .AND. K==1) ODZ(K-1) = ONE/DZ(K-1)
319           END DO
320     
321     
322           DX_E = HALF*(DX(1)+DX(IMIN1))
323           DY_N = HALF*(DY(1)+DY(JMIN1))
324           DZ_T = HALF*(DZ(1)+DZ(KMIN1))
325     
326           ODX_E(1) = ONE/DX_E
327           ODY_N(1) = ONE/DY_N
328           ODZ_T(1) = ONE/DZ_T
329     
330           FX(1) = HALF
331           FX_BAR(1) = HALF
332           FX_E(1) = HALF
333           FX_E_BAR(1) = HALF
334           FY_N(1) = HALF
335           FY_N_BAR(1) = HALF
336           FZ_T(1) = HALF
337           FZ_T_BAR(1) = HALF
338     
339           IF(NODESI.NE.1) ODX_E(IMIN3) = ONE/DX_E
340           IF(NODESJ.NE.1) ODY_N(JMIN3) = ONE/DY_N
341           IF(NODESK.NE.1) ODZ_T(KMIN3) = ONE/DZ_T
342     
343           IF(NODESI.NE.1) THEN
344              FX(IMIN3) = HALF
345              FX_BAR(IMIN3) = HALF
346              FX_E(IMIN3) = HALF
347              FX_E_BAR(IMIN3) = HALF
348           ENDIF
349     
350           IF(NODESJ.NE.1) THEN
351              FY_N(JMIN3) = HALF
352              FY_N_BAR(JMIN3) = HALF
353           ENDIF
354     
355           IF(NODESK.NE.1) THEN
356              FZ_T(KMIN3) = HALF
357              FZ_T_BAR(KMIN3) = HALF
358           ENDIF
359     
360     
361     ! Look at 2 through IMAX1 U-momentum cells
362           IF (DO_I) THEN
363              DO I = IMIN1, IMAX1
364                 DX_E = HALF*(DX(I+1)+DX(I))
365                 ODX_E(I) = ONE/DX_E
366                 FX(I) = HALF
367                 FX_BAR(I) = ONE - FX(I)
368                 FX_E(I) = DX(I+1)/(DX(I+1)+DX(I))
369                 FX_E_BAR(I) = ONE - FX_E(I)
370              END DO
371           ENDIF
372     
373     ! Look at 2 through JMAX1 V-momentum cells
374           IF (DO_J) THEN
375              DO J = JMIN1, JMAX1
376                 DY_N = HALF*(DY(J+1)+DY(J))
377                 ODY_N(J) = ONE/DY_N
378                 FY_N(J) = DY(J+1)/(DY(J+1)+DY(J))
379                 FY_N_BAR(J) = ONE - FY_N(J)
380              END DO
381           ENDIF
382     
383     ! Look at 2 through KMAX1 W-momentum cells
384           IF (DO_K) THEN
385              DO K = KMIN1, KMAX1
386                 DZ_T = HALF*(DZ(K+1)+DZ(K))
387                 ODZ_T(K) = ONE/DZ_T
388                 FZ_T(K) = DZ(K+1)/(DZ(K+1)+DZ(K))
389                 FZ_T_BAR(K) = ONE - FZ_T(K)
390              END DO
391           ENDIF
392     
393     ! Look at last U-, V-, and W-momentum cells
394           DX_E = DX(IMAX2)
395           DY_N = DY(JMAX2)
396           DZ_T = DZ(KMAX2)
397           ODX_E(IMAX2) = ONE/DX_E
398           ODY_N(JMAX2) = ONE/DY_N
399           ODZ_T(KMAX2) = ONE/DZ_T
400           FX(IMAX2) = HALF
401           FX_BAR(IMAX2) = HALF
402           FX_E(IMAX2) = HALF
403           FX_E_BAR(IMAX2) = HALF
404     
405           FY_N(JMAX2) = HALF
406           FY_N_BAR(JMAX2) = HALF
407     
408           FZ_T(KMAX2) = HALF
409           FZ_T_BAR(KMAX2) = HALF
410           FZ_T(KMAX3) = HALF
411           FZ_T_BAR(KMAX3) = HALF
412     
413           ODX_E(IMAX3) = ONE/DX_E
414           ODY_N(JMAX3) = ONE/DY_N
415           ODZ_T(KMAX3) = ONE/DZ_T
416           FX(IMAX3) = HALF
417           FX_BAR(IMAX3) = HALF
418           FX_E(IMAX3) = HALF
419           FX_E_BAR(IMAX3) = HALF
420     
421           FY_N(JMAX3) = HALF
422           FY_N_BAR(JMAX3) = HALF
423     
424     
425           IF(CYCLIC_X) THEN
426              FX_E(1) = FX_E(IMAX1)
427              FX_E_BAR(1) = FX_E_BAR(IMAX1)
428               IF(NODESI.NE.1) THEN
429                 FX_E(IMIN3) = FX_E(IMAX1-1)
430                 FX_E_BAR(IMIN3) = FX_E_BAR(IMAX1-1)
431              ENDIF
432           ENDIF
433           IF (CYCLIC_Y) THEN
434              FY_N(1) = FY_N(JMAX1)
435              FY_N_BAR(1) = FY_N_BAR(JMAX1)
436              IF(NODESJ.NE.1) THEN
437                 FY_N(JMIN3) = FY_N(JMAX1-1)
438                 FY_N_BAR(JMIN3) = FY_N_BAR(JMAX1-1)
439              ENDIF
440           ENDIF
441           IF (CYCLIC_Z) THEN
442              FZ_T(1) = FZ_T(KMAX1)
443              FZ_T_BAR(1) = FZ_T_BAR(KMAX1)
444              IF(NODESK.NE.1) THEN
445                 FZ_T(KMIN3) = FZ_T(KMAX1-1)
446                 FZ_T_BAR(KMIN3) = FZ_T_BAR(KMAX1-1)
447              ENDIF
448           ENDIF
449     
450           CALL FINL_ERR_MSG
451     
452           RETURN
453           END SUBROUTINE SET_GEOMETRY
454