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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Module name: MAKE_ARRAYS_DES                                        !
3     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
4     !                                                                      !
5     !  Purpose: DES - allocating DES arrays
6     !                                                                      !
7     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
8           SUBROUTINE MAKE_ARRAYS_DES
9     
10           USE calc_collision_wall
11           USE compar
12           USE cutcell
13           USE des_rxns
14           USE des_stl_functions
15           USE des_thermo
16           USE desgrid
17           USE discretelement
18           USE error_manager
19           USE functions
20           USE funits
21           USE geometry
22           USE mpi_utility
23           USE param1
24           USE run
25           USE stl
26           use mfix_pic, only:  MPPIC
27           use mpi_funs_des, only: DES_PAR_EXCHANGE
28           use constant, only: PI
29     
30           IMPLICIT NONE
31     !-----------------------------------------------
32     ! Local variables
33     !-----------------------------------------------
34           INTEGER :: I, J, K, L, IJK, PC, SM_CELL
35           INTEGER :: I1, I2, J1, J2, K1, K2, II, JJ, KK, IJK2
36           INTEGER :: lface, lcurpar, lpip_all(0:numpes-1), lglobal_id  , lparcnt
37           INTEGER :: CELL_ID, I_CELL, J_CELL, K_CELL, COUNT, NF
38           INTEGER :: IMINUS1, IPLUS1, JMINUS1, JPLUS1, KMINUS1, KPLUS1
39     
40           INTEGER :: FACTOR
41     
42     ! MPPIC related quantities
43           DOUBLE PRECISION :: DTPIC_TMPX, DTPIC_TMPY, DTPIC_TMPZ
44           CALL INIT_ERR_MSG("MAKE_ARRAYS_DES")
45     
46     ! Check interpolation input.
47           CALL SET_FILTER_DES
48     
49     ! cfassign and des_init_bc called before reading the particle info
50           CALL CFASSIGN
51     
52           VOL_SURR(:) = ZERO
53     
54           ! initialize VOL_SURR array
55           DO K = KSTART2, KEND1
56              DO J = JSTART2, JEND1
57                 DO I = ISTART2, IEND1
58                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
59                    IJK = funijk(I,J,K)
60                    I1 = I
61                    I2 = I+1
62                    J1 = J
63                    J2 = J+1
64                    K1 = K
65                    K2 = merge(K, K+1, NO_K)
66     
67     ! looping over stencil points (node values)
68                    DO KK = K1, K2
69                       DO JJ = J1, J2
70                          DO II = I1, I2
71                             IF (DEAD_CELL_AT(II,JJ,KK)) CYCLE  ! skip dead cells
72                             IJK2 = funijk_map_c(II, JJ, KK)
73                             IF(FLUID_AT(IJK2)) VOL_SURR(IJK) = &
74                                VOL_SURR(IJK)+VOL(IJK2)
75                          ENDDO
76                       ENDDO
77                    ENDDO
78                 ENDDO
79              ENDDO
80           ENDDO
81     
82           VERTEX(1,:,WEST_FACEID) = (/zero, zero, zero/)
83           VERTEX(2,:,WEST_FACEID) = (/zero, 2*YLENGTH, zero/)
84           VERTEX(3,:,WEST_FACEID) = (/zero, zero, 2*ZLENGTH/)
85     
86           VERTEX(1,:,EAST_FACEID) = (/XLENGTH, zero, zero/)
87           VERTEX(2,:,EAST_FACEID) = (/XLENGTH, 2*YLENGTH, zero/)
88           VERTEX(3,:,EAST_FACEID) = (/XLENGTH, zero, 2*ZLENGTH/)
89     
90           VERTEX(1,:,SOUTH_FACEID) = (/zero, zero, zero/)
91           VERTEX(2,:,SOUTH_FACEID) = (/2*XLENGTH, zero, zero/)
92           VERTEX(3,:,SOUTH_FACEID) = (/zero, zero, 2*ZLENGTH/)
93     
94           VERTEX(1,:,NORTH_FACEID) = (/zero, YLENGTH, zero/)
95           VERTEX(2,:,NORTH_FACEID) = (/2*XLENGTH, YLENGTH, zero/)
96           VERTEX(3,:,NORTH_FACEID) = (/zero, YLENGTH, 2*ZLENGTH/)
97     
98           VERTEX(1,:,BOTTOM_FACEID) = (/zero, zero, zero/)
99           VERTEX(2,:,BOTTOM_FACEID) = (/2*XLENGTH, zero, zero/)
100           VERTEX(3,:,BOTTOM_FACEID) = (/zero, 2*YLENGTH, zero/)
101     
102           VERTEX(1,:,TOP_FACEID) = (/zero, zero, ZLENGTH/)
103           VERTEX(2,:,TOP_FACEID) = (/2*XLENGTH, zero, ZLENGTH/)
104           VERTEX(3,:,TOP_FACEID) = (/zero, 2*YLENGTH, ZLENGTH/)
105     
106     
107           NORM_FACE(:,WEST_FACEID) = (/one, zero, zero/)
108           NORM_FACE(:,EAST_FACEID) = (/-one, zero, zero/)
109           NORM_FACE(:,SOUTH_FACEID) = (/zero, one, zero/)
110           NORM_FACE(:,NORTH_FACEID) = (/zero, -one, zero/)
111           NORM_FACE(:,BOTTOM_FACEID) = (/zero, zero, one/)
112           NORM_FACE(:,TOP_FACEID) = (/zero, zero, -one/)
113     
114     
115           STL_FACET_TYPE(WEST_FACEID) = FACET_TYPE_NORMAL
116           STL_FACET_TYPE(EAST_FACEID) = FACET_TYPE_NORMAL
117           STL_FACET_TYPE(NORTH_FACEID) = FACET_TYPE_NORMAL
118           STL_FACET_TYPE(SOUTH_FACEID) = FACET_TYPE_NORMAL
119           STL_FACET_TYPE(TOP_FACEID) = FACET_TYPE_NORMAL
120           STL_FACET_TYPE(BOTTOM_FACEID) = FACET_TYPE_NORMAL
121     
122     ! initialize CELLNEIGHBOR_FACET array
123           DO K_CELL = DG_KSTART2, DG_KEND2
124           DO J_CELL = DG_JSTART2, DG_JEND2
125           DO I_CELL = DG_ISTART2, DG_IEND2
126     
127              IPLUS1  =  MIN (I_CELL + 1, DG_IEND2)
128              IMINUS1 =  MAX (I_CELL - 1, DG_ISTART2)
129     
130              JPLUS1  =  MIN (J_CELL + 1, DG_JEND2)
131              JMINUS1 =  MAX (J_CELL - 1, DG_JSTART2)
132     
133              KPLUS1  =  MIN (K_CELL + 1, DG_KEND2)
134              KMINUS1 =  MAX (K_CELL - 1, DG_KSTART2)
135     
136              CELL_ID = DG_FUNIJK(I_CELL,J_CELL,K_CELL)
137     
138              IF(.NOT.DES_PERIODIC_WALLS_X)THEN
139                 IF(I_CELL == DG_IMIN1) CALL ADD_FACET(CELL_ID,WEST_FACEID)
140                 IF(I_CELL == DG_IMAX1) CALL ADD_FACET(CELL_ID,EAST_FACEID)
141              ENDIF
142     
143              IF(.NOT.DES_PERIODIC_WALLS_Y)THEN
144                 IF(J_CELL == DG_JMIN1) CALL ADD_FACET(CELL_ID,SOUTH_FACEID)
145                 IF(J_CELL == DG_JMAX1) CALL ADD_FACET(CELL_ID,NORTH_FACEID)
146              ENDIF
147     
148              IF (DO_K .AND. .NOT.DES_PERIODIC_WALLS_Z)THEN
149                 IF(K_CELL == DG_KMIN1) CALL ADD_FACET(CELL_ID,BOTTOM_FACEID)
150                 IF(K_CELL == DG_KMAX1) CALL ADD_FACET(CELL_ID,TOP_FACEID)
151              ENDIF
152     
153              DO KK = KMINUS1, KPLUS1
154              DO JJ = JMINUS1, JPLUS1
155              DO II = IMINUS1, IPLUS1
156                 IJK = DG_FUNIJK(II,JJ,KK)
157     
158                 DO COUNT = 1, LIST_FACET_AT_DES(IJK)%COUNT_FACETS
159                    NF = LIST_FACET_AT_DES(IJK)%FACET_LIST(COUNT)
160                    CALL ADD_FACET(CELL_ID,NF)
161                 ENDDO
162              ENDDO
163              ENDDO
164              ENDDO
165           ENDDO
166           ENDDO
167           ENDDO
168     
169     ! Set the initial particle data.
170           IF(RUN_TYPE == 'NEW') THEN
171     
172              IF(PARTICLES /= 0) THEN
173                 IF(GENER_PART_CONFIG) THEN
174                    CALL COPY_PARTICLE_CONFIG_FROMLISTS
175                 ELSE
176                    CALL READ_PAR_INPUT
177                 ENDIF
178              ENDIF
179     
180     ! Set the global ID for the particles and set the ghost cnt
181              ighost_cnt = 0
182              lpip_all = 0
183              lpip_all(mype) = pip
184              call global_all_sum(lpip_all)
185              lglobal_id = sum(lpip_all(0:mype-1))
186              imax_global_id = 0
187              do lcurpar  = 1,pip
188                 lglobal_id = lglobal_id + 1
189                 iglobal_id(lcurpar) = lglobal_id
190                 imax_global_id = iglobal_id(pip)
191              end do
192              call global_all_max(imax_global_id)
193     
194     ! Initialize old values
195              omega_new(:,:)   = zero
196     
197     ! Particle orientation
198              IF(PARTICLE_ORIENTATION) THEN
199                 ORIENTATION(1,:) = INIT_ORIENTATION(1)
200                 ORIENTATION(2,:) = INIT_ORIENTATION(2)
201                 ORIENTATION(3,:) = INIT_ORIENTATION(3)
202              ENDIF
203     
204     
205              IF (DO_OLD) THEN
206                 omega_old(:,:)   = zero
207                 des_pos_old(:,:) = des_pos_new(:,:)
208                 des_vel_old(:,:) = des_vel_new(:,:)
209              ENDIF
210     
211     ! Read the restart file.
212           ELSEIF(RUN_TYPE == 'RESTART_1' .OR. RUN_TYPE == 'RESTART_2') THEN
213     
214              CALL READ_RES0_DES
215              imax_global_id = maxval(iglobal_id(1:pip))
216              call global_all_max(imax_global_id)
217     
218     ! Initizlie the old values.
219              IF (DO_OLD) THEN
220                 omega_old(:,:)   = omega_new(:,:)
221                 des_pos_old(:,:) = des_pos_new(:,:)
222                 des_vel_old(:,:) = des_vel_new(:,:)
223              ENDIF
224              IF(ENERGY_EQ) DES_T_s_OLD(:) = DES_T_s_NEW(:)
225     
226           ELSE
227     
228              WRITE(ERR_MSG, 1100)
229              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
230      1100 FORMAT('Error 1100: Unsupported RUN_TYPE for DES.')
231     
232           ENDIF
233     
234           IF(RUN_TYPE == 'RESTART_2') VTP_FINDEX=0
235     
236     ! setting additional particle properties now that the particles
237     ! have been identified
238           DO L = 1, MAX_PIP
239     ! Skip 'empty' locations when populating the particle property arrays.
240              IF(.NOT.PEA(L,1)) CYCLE
241              IF(PEA(L,4)) CYCLE
242              PVOL(L) = (4.0D0/3.0D0)*PI*DES_RADIUS(L)**3
243              PMASS(L) = PVOL(L)*RO_SOL(L)
244              OMOI(L) = 2.5D0/(PMASS(L)*DES_RADIUS(L)**2) !ONE OVER MOI
245           ENDDO
246     
247           CALL SET_PHASE_INDEX
248           CALL INIT_PARTICLES_IN_CELL
249     
250     ! do_nsearch should be set before calling particle in cell
251           DO_NSEARCH =.TRUE.
252           CALL DES_PAR_EXCHANGE
253           CALL PARTICLES_IN_CELL
254     
255           IF(DEM_SOLIDS) THEN
256              CALL NEIGHBOUR
257              CALL INIT_SETTLING_DEM
258           ENDIF
259     
260     ! Calculate interpolation weights
261           CALL CALC_INTERP_WEIGHTS
262     ! Calculate mean fields using either interpolation or cell averaging.
263           CALL COMP_MEAN_FIELDS
264     
265           IF(MPPIC) CALL CALC_DTPIC
266     
267           CALL FINL_ERR_MSG
268     
269           RETURN
270           END SUBROUTINE MAKE_ARRAYS_DES
271