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