File: N:\mfix\model\des\pic\mass_inflow_pic.f
1
2
3
4
5
6
7 SUBROUTINE MASS_INFLOW_PIC
8
9
10
11 USE bc
12 USE constant
13 USE cutcell
14 USE discretelement
15 use error_manager
16 USE functions
17 USE geometry
18 USE mfix_pic
19 USE mpi_utility
20 USE param, only: dimension_m
21 USE param1, only: half, zero
22 use physprop, only: mmax, D_p0, ro_s0
23 USE pic_bc
24 USE randomno
25 USE run
26 IMPLICIT NONE
27
28
29
30 INTEGER :: WDIR, IDIM, IPCOUNT, LAST_EMPTY_SPOT, NEW_SPOT
31 INTEGER :: BCV, BCV_I, L, LC, PIP_ADD_COUNT, IPROC
32 INTEGER :: IFLUID, JFLUID, KFLUID, IJK, M
33 DOUBLE PRECISION :: CORD_START(3), DOML(3),WALL_NORM(3)
34 DOUBLE PRECISION :: AREA_INFLOW, VEL_INFLOW(DIMN), STAT_WT
35
36 LOGICAL :: DELETE_PART
37 DOUBLE PRECISION :: EPS_INFLOW(DIMENSION_M)
38 DOUBLE PRECISION :: REAL_PARTS(DIMENSION_M)
39 DOUBLE PRECISION :: COMP_PARTS(DIMENSION_M)
40 DOUBLE PRECISION :: VEL_NORM_MAG, VOL_INFLOW, VOL_IJK
41
42
43 LOGICAL :: CONST_NPC, CONST_STATWT
44
45 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RANDPOS
46 INTEGER :: CNP_CELL_COUNT
47 DOUBLE PRECISION :: RNP_CELL_COUNT
48 integer :: lglobal_id
49 integer, dimension(0:numpes-1) :: add_count_all
50
51 LOGICAL, parameter :: setDBG = .true.
52 LOGICAL :: dFlag
53
54 type :: ty_spotlist
55 integer spot
56 type(ty_spotlist),pointer :: next => NULL()
57 end type ty_spotlist
58
59 type(ty_spotlist),pointer :: &
60 cur_spotlist => NULL(), &
61 prev_spotlist => NULL(), &
62 temp_spotlist => NULL()
63
64
65
66 CALL INIT_ERR_MSG("PIC_MI_BC")
67
68 dFlag = (DMP_LOG .AND. setDBG)
69 PIP_ADD_COUNT = 0
70 LAST_EMPTY_SPOT = 0
71
72 ALLOCATE(cur_spotlist); cur_spotlist%spot = -1
73
74 DO BCV_I = 1, PIC_BCMI
75 BCV = PIC_BCMI_MAP(BCV_I)
76
77 WALL_NORM(1:3) = PIC_BCMI_NORMDIR(BCV_I,1:3)
78
79
80 = 0
81 DO IDIM = 1, DIMN
82 WDIR = WDIR + ABS(WALL_NORM(IDIM))*IDIM
83 end DO
84
85 DO LC=PIC_BCMI_IJKSTART(BCV_I), PIC_BCMI_IJKEND(BCV_I)
86 IJK = PIC_BCMI_IJK(LC)
87
88 IF(.NOT.FLUID_AT(IJK)) CYCLE
89
90 IFLUID = I_OF(IJK)
91 JFLUID = J_OF(IJK)
92 KFLUID = K_OF(IJK)
93
94 CORD_START(1) = XE(IFLUID) - PIC_BCMI_OFFSET (BCV_I,1)*DX(IFLUID)
95
96 CORD_START(2) = YN(JFLUID) - PIC_BCMI_OFFSET (BCV_I,2)*DY(JFLUID)
97
98
99 CORD_START(3) = merge(zero, ZT(KFLUID) - PIC_BCMI_OFFSET (BCV_I,3)*DZ(KFLUID), no_k)
100
101 DOML(1) = DX(IFLUID)
102 DOML(2) = DY(JFLUID)
103 DOML(3) = MERGE(DZ(1), DZ(KFLUID), NO_K)
104
105 AREA_INFLOW = DOML(1)*DOML(2)*DOML(3)/DOML(WDIR)
106
107 VOL_IJK = DOML(1)*DOML(2)*DOML(3)
108
109 DOML(WDIR) = ZERO
110
111
112
113 DO M = 1, DES_MMAX+MMAX
114
115 IF(SOLIDS_MODEL(M) /= 'PIC') CYCLE
116
117 EPS_INFLOW(M) = BC_ROP_S(BCV, M)/RO_S0(M)
118 VEL_INFLOW(1) = BC_U_S(BCV, M)
119 VEL_INFLOW(2) = BC_V_S(BCV, M)
120 VEL_INFLOW(3) = BC_W_S(BCV, M)
121
122 VEL_NORM_MAG = ABS(DOT_PRODUCT(VEL_INFLOW(1:DIMN), WALL_NORM(1:DIMN)))
123 VOL_INFLOW = AREA_INFLOW*VEL_NORM_MAG*DTSOLID
124
125 REAL_PARTS(M) = 6.d0*EPS_INFLOW(M)*VOL_INFLOW/&
126 (PI*(D_p0(M)**3.d0))
127 COMP_PARTS(M) = zero
128
129 CONST_NPC = (BC_PIC_MI_CONST_NPC (BCV, M) .ne. 0)
130 CONST_STATWT = (BC_PIC_MI_CONST_STATWT(BCV, M) .ne. ZERO)
131 IF(CONST_NPC) THEN
132 IF(EPS_INFLOW(M).GT.ZERO) &
133 COMP_PARTS(M) = REAL(BC_PIC_MI_CONST_NPC(BCV, M))* &
134 VOL_INFLOW/VOL_IJK
135 ELSEIF(CONST_STATWT) THEN
136 COMP_PARTS(M) = REAL_PARTS(M)/ &
137 BC_PIC_MI_CONST_STATWT(BCV, M)
138 ENDIF
139
140 pic_bcmi_rnp(LC,M) = pic_bcmi_rnp(LC,M) + REAL_PARTS(M)
141
142 pic_bcmi_cnp(LC,M) = pic_bcmi_cnp(LC,M) + COMP_PARTS(M)
143
144
145 IF(pic_bcmi_cnp(LC,M).GE.1.d0) THEN
146 CNP_CELL_COUNT = INT(pic_bcmi_cnp(LC,M))
147 pic_bcmi_cnp(LC,M) = pic_bcmi_cnp(LC,M) - CNP_CELL_COUNT
148
149 RNP_CELL_COUNT = pic_bcmi_rnp(LC,M)
150 pic_bcmi_rnp(LC,M) = zero
151
152
153 = RNP_CELL_COUNT/REAL(CNP_CELL_COUNT)
154 ALLOCATE(RANDPOS(CNP_CELL_COUNT*DIMN))
155 CALL UNI_RNO(RANDPOS(:))
156
157 DO IPCOUNT = 1, CNP_CELL_COUNT
158
159
160 CALL PIC_FIND_EMPTY_SPOT(LAST_EMPTY_SPOT, NEW_SPOT)
161
162 DES_POS_OLD( NEW_SPOT,1:DIMN) = CORD_START(1:DIMN) &
163 & + RANDPOS((IPCOUNT-1)*DIMN+1: &
164 & (IPCOUNT-1)*DIMN+DIMN)*DOML(1:DIMN)
165 DES_POS_NEW( NEW_SPOT,:) = DES_POS_OLD(NEW_SPOT,:)
166 DES_VEL_OLD(NEW_SPOT,1:DIMN) = VEL_INFLOW(1:DIMN)
167
168 DES_VEL_NEW( NEW_SPOT,:) = DES_VEL_OLD( NEW_SPOT,:)
169
170 DES_RADIUS(NEW_SPOT) = D_p0(M)*HALF
171
172 RO_Sol(NEW_SPOT) = RO_S0(M)
173
174 DES_STAT_WT(NEW_SPOT) = STAT_WT
175
176 PIJK(NEW_SPOT, 1) = IFLUID
177 PIJK(NEW_SPOT, 2) = JFLUID
178 PIJK(NEW_SPOT, 3) = KFLUID
179 PIJK(NEW_SPOT, 4) = IJK
180 PIJK(NEW_SPOT, 5) = M
181
182 PVOL(NEW_SPOT) = (4.0d0/3.0d0)*Pi*DES_RADIUS(NEW_SPOT)**3
183 PMASS(NEW_SPOT) = PVOL(NEW_SPOT)*RO_SOL(NEW_SPOT)
184
185 FC( NEW_SPOT,:) = zero
186 DELETE_PART = .false.
187 IF(PIC_BCMI_INCL_CUTCELL(BCV_I)) &
188 CALL CHECK_IF_PARCEL_OVERLAPS_STL &
189 (des_pos_new( NEW_SPOT,1:dimn), &
190 DELETE_PART)
191
192 IF(.NOT.DELETE_PART) THEN
193
194 PIP = PIP+1
195 PIP_ADD_COUNT = PIP_ADD_COUNT + 1
196 CALL SET_NORMAL(NEW_SPOT)
197
198 ALLOCATE(temp_spotlist)
199 temp_spotlist%spot = new_spot
200 temp_spotlist%next => cur_spotlist
201 cur_spotlist => temp_spotlist
202 nullify(temp_spotlist)
203
204 ELSE
205 CALL SET_NONEXISTENT(NEW_SPOT)
206 LAST_EMPTY_SPOT = NEW_SPOT - 1
207 ENDIF
208
209
210
211
212
213
214
215 END DO
216 DEALLOCATE(RANDPOS)
217
218 end IF
219 end DO
220 end DO
221 end DO
222
223
224 (:) = 0
225 add_count_all(mype) = pip_add_count
226 call global_all_sum(add_count_all(0:numpes-1))
227 lglobal_id = imax_global_id + sum(add_count_all(0:mype-1))
228
229 do l = 1,pip_add_count
230 lglobal_id = lglobal_id + 1
231 iglobal_id(cur_spotlist%spot)= lglobal_id
232 prev_spotlist=> cur_spotlist
233 cur_spotlist => cur_spotlist%next
234 deallocate(prev_spotlist)
235 end do
236 deallocate(cur_spotlist)
237 imax_global_id = imax_global_id + sum(add_count_all(0:numpes-1))
238
239 IF(PIC_REPORT_SEEDING_STATS) then
240
241 IF(SUM(ADD_COUNT_ALL(:)).GT.0) THEN
242 WRITE(err_msg,'(/,2x,A,2x,i10)') &
243 'TOTAL NUMBER OF PARCELS ADDED GLOBALLY = ',&
244 SUM(ADD_COUNT_ALL(:))
245
246 call flush_err_msg(header = .false., footer = .false.)
247
248 DO IPROC = 0, NUMPES-1
249 IF(DMP_LOG) WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
250 'PARCELS ADDED ON PROC:', IPROC,&
251 ' EQUAL TO', ADD_COUNT_ALL(IPROC)
252 ENDDO
253 ENDIF
254
255 ENDIF
256
257 CALL FINL_ERR_MSG
258 RETURN
259 END SUBROUTINE MASS_INFLOW_PIC
260
261
262
263
264
265
266
267
268
269 SUBROUTINE PIC_FIND_EMPTY_SPOT(LAST_INDEX, EMPTY_SPOT)
270
271
272
273 USE funits
274 USE mpi_utility
275 USE error_manager
276 USE discretelement, only: max_pip
277 USE functions, only: is_nonexistent
278 IMPLICIT NONE
279
280
281
282 INTEGER, INTENT(INOUT) :: LAST_INDEX
283 INTEGER, INTENT(OUT) :: EMPTY_SPOT
284
285
286
287 LOGICAL :: SPOT_FOUND
288 INTEGER :: LL
289
290
291 CALL INIT_ERR_MSG("PIC_FIND_EMPTY_SPOT")
292
293 IF(LAST_INDEX.EQ.MAX_PIP) THEN
294 WRITE(ERR_MSG,2001)
295
296 CALL FLUSH_ERR_MSG(abort = .true.)
297 call mfix_exit(mype)
298 ENDIF
299 SPOT_FOUND = .false.
300
301 DO LL = LAST_INDEX+1, MAX_PIP
302
303 if(IS_NONEXISTENT(LL)) THEN
304 EMPTY_SPOT = LL
305 LAST_INDEX = LL
306 SPOT_FOUND = .true.
307 EXIT
308 ENDIF
309 ENDDO
310
311 IF(.NOT.SPOT_FOUND) THEN
312 WRITE(ERR_MSG,2002)
313 CALL FLUSH_ERR_MSG(abort = .true.)
314 ENDIF
315
316 2001 FORMAT(/,5X, &
317 & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5X, &
318 & 'NO MORE EMPTY SPOT IN THE PARTICLE ARRAY TO ADD A NEW PARTICLE',/5X, &
319 & 'TERMINAL ERROR: STOPPING')
320
321 2002 FORMAT(/,5X, &
322 & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5X, &
323 & 'COULD NOT FIND A SPOT FOR ADDING NEW PARTICLE',/5X, &
324 & 'INCREASE THE SIZE OF THE INITIAL ARRAYS', 5X, &
325 & 'TERMINAL ERROR: STOPPING')
326
327 CALL FINL_ERR_MSG
328 END SUBROUTINE PIC_FIND_EMPTY_SPOT
329
330
331
332
333
334
335
336
337
338
339
340 SUBROUTINE PIC_FIND_NEW_CELL(LL)
341
342
343
344 USE discretelement
345 use mpi_utility
346 USE cutcell
347 USE functions
348 IMPLICIT NONE
349
350
351
352 INTEGER, INTENT(IN) :: LL
353
354
355
356 DOUBLE PRECISION :: XPOS, YPOS, ZPOS
357 INTEGER :: I, J, K
358
359
360 = PIJK(LL,1)
361 J = PIJK(LL,2)
362 K = PIJK(LL,3)
363 XPOS = DES_POS_NEW(LL,1)
364 YPOS = DES_POS_NEW(LL,2)
365 IF (DIMN .EQ. 3) THEN
366 ZPOS = DES_POS_NEW(LL,3)
367 ENDIF
368
369 IF(XPOS >= XE(I-1) .AND. XPOS < XE(I)) THEN
370 PIJK(LL,1) = I
371 ELSEIF(XPOS >= XE(I)) THEN
372 PIJK(LL,1) = I+1
373 ELSE
374 PIJK(LL,1) = I-1
375 END IF
376
377 IF(YPOS >= YN(J-1) .AND. YPOS < YN(J)) THEN
378 PIJK(LL,2) = J
379 ELSEIF(YPOS >= YN(J))THEN
380 PIJK(LL,2) = J+1
381 ELSE
382 PIJK(LL,2) = J-1
383 END IF
384
385 IF(DIMN.EQ.2) THEN
386 PIJK(LL,3) = 1
387 ELSE
388 IF(ZPOS >= ZT(K-1) .AND. ZPOS < ZT(K)) THEN
389 PIJK(LL,3) = K
390 ELSEIF(ZPOS >= ZT(K)) THEN
391 PIJK(LL,3) = K+1
392 ELSE
393 PIJK(LL,3) = K-1
394 END IF
395 ENDIF
396
397 I = PIJK(LL,1)
398 J = PIJK(LL,2)
399 K = PIJK(LL,3)
400
401 IF(I.EQ.IEND1+1) then
402 IF(XPOS>=XE(IEND1-1) .AND. XPOS<=XE(IEND1)) PIJK(LL,1) = IEND1
403 ENDIF
404
405 IF(J.EQ.JEND1+1) then
406 IF(YPOS>=YN(JEND1-1) .AND. YPOS<=YN(JEND1)) PIJK(LL,2) = JEND1
407 ENDIF
408
409 IF(DIMN.EQ.3.AND.K.EQ.KEND1+1) THEN
410 IF(ZPOS>=ZT(KEND1-1) .AND. ZPOS<=ZT(KEND1)) PIJK(LL,3) = KEND1
411 ENDIF
412
413 PIJK(LL,4) = FUNIJK(PIJK(LL,1), PIJK(LL,2),PIJK(LL,3))
414
415 END SUBROUTINE PIC_FIND_NEW_CELL
416
417