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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: DES_MASS_INLET                                          !
4     !  Author: J.Musser                                   Date: 13-Jul-09  !
5     !                                                                      !
6     !  Purpose:  This routine fills in the necessary information for new   !
7     !  particles entereing the system.                                     !
8     !                                                                      !
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
10           SUBROUTINE MASS_INFLOW_DEM
11     
12           use bc
13           use des_allocate
14           use des_bc
15           use discretelement
16     
17           use mpi_utility, only: BCAST
18     
19           implicit none
20     
21           INTEGER :: IP, LS, M, NP, IJK, LC
22           INTEGER :: BCV, BCV_I
23           LOGICAL :: CHECK_FOR_ERRORS, OWNS
24     
25     ! I/J/K index of fluid cell containing the new particle.
26           INTEGER :: IJKP(3)
27     
28           DOUBLE PRECISION :: DIST, POS(3)
29     ! Random numbers -shared by all ranks-
30           DOUBLE PRECISION :: RAND(3)
31     
32           CHECK_FOR_ERRORS = .FALSE.
33     
34           DO BCV_I = 1, DEM_BCMI
35              BCV = DEM_BCMI_MAP(BCV_I)
36     
37              DO LC=DEM_BCMI_IJKSTART(BCV_I), DEM_BCMI_IJKEND(BCV_I)
38                IJK = DEM_BCMI_IJK(LC)
39     
40                 DO LS= 1,PINC(IJK)
41                    NP = PIC(IJK)%p(LS)
42                    IF(PEA(NP,3)) CYCLE
43                    SELECT CASE (BC_PLANE(BCV))
44                    CASE('N'); DIST = DES_POS_NEW(2,NP) - YN(BC_J_s(BCV))
45                    CASE('S'); DIST = YN(BC_J_s(BCV)-1) - DES_POS_NEW(2,NP)
46                    CASE('E'); DIST = DES_POS_NEW(1,NP) - XE(BC_I_w(BCV))
47                    CASE('W'); DIST = XE(BC_I_w(BCV)-1) - DES_POS_NEW(1,NP)
48                    CASE('T'); DIST = DES_POS_NEW(3,NP) - ZT(BC_K_b(BCV))
49                    CASE('B'); DIST = ZT(BC_K_b(BCV)-1) - DES_POS_NEW(3,NP)
50                    END SELECT
51     ! The particle is still inside the domain
52                    IF(DIST > DES_RADIUS(NP)) PEA(NP,2) = .FALSE.
53                 ENDDO
54              ENDDO
55     
56     ! All ranks generate random numbers but PE_IO BCASTS its values so that
57     ! the calculated position (with random wiggles) is the same on all ranks.
58              CALL RANDOM_NUMBER(RAND)
59              CALL BCAST(RAND)
60     
61     ! Check if any particles need seeded this time step.
62              IF(DEM_MI_TIME(BCV_I) > S_TIME) CYCLE
63     
64              LS = 1
65     
66     ! Loop over the particles being injected
67              PLoop: DO IP = 1, PI_COUNT(BCV_I)
68     
69     ! Increment the global max particle ID (all ranks).
70                 iMAX_GLOBAL_ID = iMAX_GLOBAL_ID + 1
71     
72     ! Determine the location and phase of the incoming particle.
73                 CALL SEED_NEW_PARTICLE(BCV, BCV_I, RAND, M, POS, IJKP, OWNS)
74     
75     ! Only the rank receiving the new particle needs to continue
76                 IF(.NOT.OWNS) CYCLE PLoop
77     
78     ! Increment the number of particle on the processor by one. If the max
79     ! number of particles is exceeded, set the error flag and cycle.
80                 PIP = PIP + 1
81                 CALL PARTICLE_GROW(PIP)
82     
83     ! Find the first free space in the particle existance array.
84                 NP_LP: DO NP = LS, MAX_PIP
85                    IF(.NOT.PEA(NP,1)) THEN
86                       LS = NP
87                       EXIT NP_LP
88                    ENDIF
89                 ENDDO NP_LP
90     
91     ! Set the particle's global ID.
92                 iGLOBAL_ID(LS) = iMAX_GLOBAL_ID
93     
94     ! Set the properties of the new particle.
95                 CALL SET_NEW_PARTICLE_PROPS(BCV, M, LS, POS, IJKP)
96     
97              ENDDO PLoop
98     
99     ! Update the time for seeding the next particle.
100              DEM_MI_TIME(BCV_I) = S_TIME + PI_FACTOR(BCV_I)*DTSOLID
101     ! Set the flag for error checking.
102              CHECK_FOR_ERRORS = .TRUE.
103           ENDDO
104     
105           IF(CHECK_FOR_ERRORS) THEN
106           ENDIF
107     
108           RETURN
109           END SUBROUTINE MASS_INFLOW_DEM
110     
111     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
112     !                                                                      !
113     !  Subroutine: SEED_NEW_PARTICLE                                       !
114     !  Author: J.Musser                                   Date: 14-Aug-09  !
115     !                                                                      !
116     !  Purpose:  This routine uses the classification information to place !
117     !  a new particle in the proper location.                              !
118     !                                                                      !
119     !  Comments:                                                           !
120     !                                                                      !
121     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
122           SUBROUTINE SEED_NEW_PARTICLE(lBCV, lBCV_I, lRAND, lM, lPOS,      &
123              lIJKP, lOWNS)
124     
125     !-----------------------------------------------
126     ! Modules
127     !-----------------------------------------------
128           USE compar
129           USE bc
130           USE des_bc
131           USE desgrid
132           USE discretelement
133           USE funits
134           USE geometry
135           USE param1
136           USE physprop
137           IMPLICIT NONE
138     !-----------------------------------------------
139     ! Dummy arguments
140     !-----------------------------------------------
141     ! The associated bc no.
142           INTEGER, INTENT(IN) :: lBCV, lBCV_I
143     ! Random numbers
144           DOUBLE PRECISION, INTENT(IN) :: lRAND(3)
145     ! Phase of incoming particle.
146           INTEGER, INTENT(OUT) :: lM
147     ! Position of incoming particle.
148           DOUBLE PRECISION, INTENT(OUT) :: lPOS(3)
149     ! I/J/K index of fluid cell containing the new particle.
150           INTEGER, INTENT(OUT) :: lIJKP(3)
151     ! Logical indicating that the current rank is the owner
152           LOGICAL, INTENT(OUT) :: lOWNS
153     
154     !-----------------------------------------------
155     ! Local variables
156     !-----------------------------------------------
157     ! the associated bc no.
158     !      INTEGER :: BCV
159     ! Random position
160           DOUBLE PRECISION RAND1, RAND2
161     ! Array index of vacant position
162           INTEGER :: VACANCY
163     ! Number of array position.
164           INTEGER :: OCCUPANTS
165     
166           INTEGER :: RAND_I
167           INTEGER :: lI, lJ, lK
168           INTEGER :: llI, llJ, llK
169     
170           DOUBLE PRECISION :: WINDOW
171     
172     
173     !      IF(PARTICLE_PLCMNT(lBCV_I) == 'ORDR')THEN
174              VACANCY = DEM_MI(lBCV_I)%VACANCY
175              OCCUPANTS = DEM_MI(lBCV_I)%OCCUPANTS
176              DEM_MI(lBCV_I)%VACANCY = MOD(VACANCY,OCCUPANTS) + 1
177     !      ELSE
178     !         call bcast(lpar_rad)
179     !         call bcast(lpar_pos)
180     !         call bcast(m)
181     !      ENDIF
182     
183     
184     ! Assign the phase of the incoming particle.
185           IF(DEM_MI(lBCV_I)%POLYDISPERSE) THEN
186              RAND_I = ceiling(dble(NUMFRAC_LIMIT)*lRAND(1))
187              lM = DEM_BC_POLY_LAYOUT(lBCV_I, RAND_I)
188           ELSE
189              lM = DEM_BC_POLY_LAYOUT(lBCV_I,1)
190           ENDIF
191     
192           WINDOW = DEM_MI(lBCV_I)%WINDOW
193           RAND1 = HALF*DES_D_p0(lM) + (WINDOW - DES_D_p0(lM))*lRAND(1)
194           RAND2 = HALF*DES_D_p0(lM) + (WINDOW - DES_D_p0(lM))*lRAND(2)
195     
196     
197     ! Set the physical location and I/J/K location of the particle.
198           SELECT CASE(BC_PLANE(lBCV))
199           CASE('N','S')
200     
201              lPOS(1) = DEM_MI(lBCV_I)%P(VACANCY) + RAND1
202              lPOS(3) = DEM_MI(lBCV_I)%Q(VACANCY) + RAND2
203              lPOS(2) = DEM_MI(lBCV_I)%OFFSET
204     
205              lIJKP(1) = DEM_MI(lBCV_I)%W(VACANCY)
206              lIJKP(3) = DEM_MI(lBCV_I)%H(VACANCY)
207              lIJKP(2) = DEM_MI(lBCV_I)%L
208     
209           CASE('E','W')
210     
211              lPOS(2) = DEM_MI(lBCV_I)%P(VACANCY) + RAND1
212              lPOS(3) = DEM_MI(lBCV_I)%Q(VACANCY) + RAND2
213              lPOS(1) = DEM_MI(lBCV_I)%OFFSET
214     
215              lIJKP(2) = DEM_MI(lBCV_I)%W(VACANCY)
216              lIJKP(3) = DEM_MI(lBCV_I)%H(VACANCY)
217              lIJKP(1) = DEM_MI(lBCV_I)%L
218     
219           CASE('T','B')
220     
221              lPOS(1) = DEM_MI(lBCV_I)%P(VACANCY) + RAND1
222              lPOS(2) = DEM_MI(lBCV_I)%Q(VACANCY) + RAND2
223              lPOS(3) = DEM_MI(lBCV_I)%OFFSET
224     
225              lIJKP(1) = DEM_MI(lBCV_I)%W(VACANCY)
226              lIJKP(2) = DEM_MI(lBCV_I)%H(VACANCY)
227              lIJKP(3) = DEM_MI(lBCV_I)%L
228     
229           END SELECT
230     
231     ! Easier and cleaner to clear out the third component at the end.
232           IF(NO_K) lPOS(3) = ZERO
233     
234           lI = IofPOS(lPOS(1))
235           lJ = JofPOS(lPOS(2))
236           lK = KofPOS(lPOS(3))
237     
238           lOWNS = ((DG_ISTART <= lI) .AND. (lI <= DG_IEND) .AND.           &
239              (DG_JSTART <= lJ) .AND. (lJ <= DG_JEND) .AND.                 &
240              (DG_KSTART <= lK) .AND. (lK <= DG_KEND))
241     
242           RETURN
243           END SUBROUTINE SEED_NEW_PARTICLE
244     
245     
246     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
247     !                                                                      !
248     !  Subroutine: SET_NEW_PARTICLE_PROPS                                  !
249     !  Author: J.Musser                                   Date: 14-Aug-09  !
250     !                                                                      !
251     !  Purpose:  Set the properties of the new particle.                   !
252     !                                                                      !
253     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
254           SUBROUTINE SET_NEW_PARTICLE_PROPS(lBCV, lM, lNP, lPOS, lIJKP)
255     
256           USE compar
257           USE bc
258           USE des_bc
259           USE discretelement
260           USE funits
261           USE geometry
262           USE param1
263           USE physprop
264           USE des_thermo
265           use des_rxns
266     
267           use run, only: ENERGY_EQ
268           use run, only: ANY_SPECIES_EQ
269           use constant, only: PI
270     
271           use desgrid
272           use indices
273           use functions
274     
275           IMPLICIT NONE
276     !-----------------------------------------------
277     ! Dummy arguments
278     !-----------------------------------------------
279     ! The associated bc no.
280           INTEGER, INTENT(IN) :: lBCV
281     ! Phase of incoming particle.
282           INTEGER, INTENT(IN) :: lM
283     ! Index of new particle
284           INTEGER, INTENT(IN) :: lNP
285     ! Position of incoming particle.
286           DOUBLE PRECISION, INTENT(IN) :: lPOS(3)
287     ! I/J/K index of fluid cell containing the new particle.
288           INTEGER, INTENT(IN) :: lIJKP(3)
289     ! I/J/K index of DES grid cell
290           INTEGER :: lI, lJ, lK
291     
292     ! Global phase index
293           INTEGER :: BC_M
294     
295     ! Shift the phase index by SMAX to refernece global variables.
296           BC_M = lM + SMAX
297     
298     ! Set the PEA Flags:
299           PEA(lNP,1:2) = .TRUE.  ! The particle exists and is entering
300           PEA(lNP,3:4) = .FALSE. ! It is not exiting nor a ghost particle
301     
302     ! Set the initial position values based on mass inlet class
303           DES_POS_NEW(:,lNP) = lPOS(:)
304           PPOS(:,lNP) = lPOS(:)
305           DES_VEL_NEW(1,lNP) = BC_U_s(lBCV,BC_M)
306           DES_VEL_NEW(2,lNP) = BC_V_s(lBCV,BC_M)
307           DES_VEL_NEW(3,lNP) = BC_W_s(lBCV,BC_M)
308     
309     ! Set the initial velocity values
310           IF (DO_OLD) THEN
311              DES_POS_OLD(:,lNP) = lPOS(:)
312              DES_VEL_OLD(:,lNP) = DES_VEL_NEW(:,lNP)
313              OMEGA_OLD(:,lNP) = 0
314           ENDIF
315     
316     ! Set the initial angular velocity values
317           OMEGA_NEW(:,lNP) = 0
318     
319     ! Set the initial angular position values
320           IF(PARTICLE_ORIENTATION) ORIENTATION(1:3,lNP) = INIT_ORIENTATION
321     
322     ! Set the particle radius value
323           DES_RADIUS(lNP) = HALF * DES_D_P0(lM)
324     
325     ! Set the particle density value
326           RO_Sol(lNP) = DES_RO_S(lM)
327     
328     ! Store the I/J/K indices of the particle.
329           PIJK(lNP,1:3) = lIJKP(:)
330           PIJK(lNP,4) = FUNIJK(lIJKP(1), lIJKP(2), lIJKP(3))
331     
332     ! Set the particle mass phase
333           PIJK(lNP,5) = lM
334     
335     ! Calculate the DES grid cell indices.
336           lI = min(DG_IEND2,max(DG_ISTART2,IOFPOS(DES_POS_NEW(1,lNP))))
337           lJ = min(DG_JEND2,max(DG_JSTART2,JOFPOS(DES_POS_NEW(2,lNP))))
338           IF(NO_K) THEN
339              lK = 1
340           ELSE
341              lK = min(DG_KEND2,max(DG_KSTART2,KOFPOS(DES_POS_NEW(3,lNP))))
342           ENDIF
343     ! Store the triple
344           DG_PIJK(lNP) = DG_FUNIJK(lI,lJ,lK)
345     
346     ! Calculate the new particle's Volume, Mass, OMOI
347           PVOL(lNP) = (4.0d0/3.0d0) * PI * DES_RADIUS(lNP)**3
348           PMASS(lNP) = PVOL(lNP) * RO_Sol(lNP)
349           OMOI(lNP) = 5.d0 / (2.d0 * PMASS(lNP) * DES_RADIUS(lNP)**2)
350     
351     ! Clear the drag force
352           IF(DES_EXPLICITLY_COUPLED) DRAG_FC(:,lNP) = ZERO
353     
354     ! If solving the energy equations, set the temperature
355           IF(ANY_SPECIES_EQ .OR. ENERGY_EQ ) THEN
356              DES_T_s_NEW(lNP) = BC_T_s(lBCV,BC_M)
357              DES_T_s_OLD(lNP) = DES_T_s_NEW(lNP)
358           ENDIF
359     
360     ! Set species mass fractions
361           IF((ENERGY_EQ .AND. C_PS0(BC_M)/=UNDEFINED) .OR. ANY_SPECIES_EQ)&
362              DES_X_s(lNP,1:NMAX(BC_M)) = BC_X_s(lBCV,BC_M,1:NMAX(BC_M))
363     
364     ! Calculate time dependent physical properties
365           CALL DES_PHYSICAL_PROP(lNP, .FALSE.)
366     
367     
368           RETURN
369           END SUBROUTINE SET_NEW_PARTICLE_PROPS
370     
371     
372     
373     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
374     !                                                                      !
375     !  Subroutine:  DES_NEW_PARTICLE_TEST                                  !
376     !                                                                      !
377     !  Purpose:  This routine checks if a new particle placed using the    !
378     !  random inlet was placed in contact with an existing particle.  If   !
379     !  so a flag is set indicating contact, and the new particle is        !
380     !  repositioned within the inlet domain.                               !
381     !                                                                      !
382     !  Author: J.Musser                                   Date: 14-Aug-09  !
383     !                                                                      !
384     !  Purpose: This routine has to be modified for parallel version       !
385     !           the parameter now accepts the lpar_rad and lpar_pos and    !
386     !           tests if it touches any particles                          !
387     !  Comments:                                                           !
388     !                                                                      !
389     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
390     
391           SUBROUTINE DES_NEW_PARTICLE_TEST(BCV_I,ppar_rad,ppar_pos,TOUCHING)
392     
393           USE compar
394           USE constant
395           USE des_bc
396           USE discretelement
397           USE funits
398           USE geometry
399           USE indices
400           USE param1
401           USE physprop
402           USE functions
403     
404           IMPLICIT NONE
405     !-----------------------------------------------
406     ! Dummy arguments
407     !-----------------------------------------------
408     ! index of boundary condition
409           INTEGER, INTENT(IN) :: BCV_I
410           DOUBLE PRECISION, INTENT(IN) :: ppar_pos(DIMN)
411           DOUBLE PRECISION, INTENT(IN) :: ppar_rad
412           LOGICAL, INTENT(INOUT) :: TOUCHING
413     !-----------------------------------------------
414     ! Local variables
415     !-----------------------------------------------
416     ! particle number id of a potential overlapping/contacting particle
417           INTEGER NP2
418     ! total number of particles in current ijk cell and loop counter
419           INTEGER NPG, LL
420     ! i, j, k indices along boundary used for loop counters
421           INTEGER I, J, K, IJK
422     ! for parallel processing
423           integer listart,liend,ljstart,ljend,lkstart,lkend
424     
425           DOUBLE PRECISION  DISTVEC(DIMN), DIST, R_LM
426     !-----------------------------------------------
427     
428           TOUCHING = .FALSE.
429     
430     ! For parallel processing the arrays has to be limited
431     !      select case (des_mi_class(bcv_i))
432     !      case ('XW','XE', 'YZw','YZe')
433     !         listart = gs_array(bcv_i,1)
434     !         liend = gs_array(bcv_i,2)
435     !         ljstart = max(gs_array(bcv_i,3),jstart)
436     !         ljend = min(gs_array(bcv_i,4),jend)
437     !         lkstart = max(gs_array(bcv_i,5),jstart)
438     !         lkend = min(gs_array(bcv_i,6),jend)
439     !      case ('YN','YS', 'XZn','XZs')
440     !         listart = max(gs_array(bcv_i,1),istart)
441     !         liend = min(gs_array(bcv_i,2),iend)
442     !         ljstart = gs_array(bcv_i,3)
443     !         ljend = gs_array(bcv_i,4)
444     !         lkstart = max(gs_array(bcv_i,5),jstart)
445     !         lkend = min(gs_array(bcv_i,6),jend)
446     !      case ('ZT','ZB', 'XYt','XYb')
447     !         listart = max(gs_array(bcv_i,1),istart)
448     !         liend = min(gs_array(bcv_i,2),iend)
449     !         ljstart = max(gs_array(bcv_i,3),jstart)
450     !         ljend = min(gs_array(bcv_i,4),jend)
451     !         lkstart = gs_array(bcv_i,5)
452     !         lkend = gs_array(bcv_i,6)
453     !      end select
454     
455           DO k = lkstart,lkend
456           DO j = ljstart,ljend
457           DO i = listart,liend
458     !      DO K = GS_ARRAY(BCV_I,5), GS_ARRAY(BCV_I,6)
459     !         DO J = GS_ARRAY(BCV_I,3), GS_ARRAY(BCV_I,4)
460     !           DO I =  GS_ARRAY(BCV_I,1), GS_ARRAY(BCV_I,2)
461                  IJK = FUNIJK(I,J,K)
462                  IF(ASSOCIATED(PIC(IJK)%P)) THEN
463                    NPG =  SIZE(PIC(IJK)%P)
464                    DO LL = 1, NPG
465                       NP2 = PIC(IJK)%P(LL)
466                       DISTVEC(:) = ppar_pos(:) - DES_POS_NEW(:,NP2)
467                       DIST = DOT_PRODUCT(DISTVEC,DISTVEC)
468                       R_LM = ppar_rad + DES_RADIUS(NP2)
469                       IF(DIST .LE. R_LM*R_LM) TOUCHING = .TRUE.
470                    ENDDO
471                  ENDIF
472                ENDDO
473              ENDDO
474            ENDDO
475     
476           RETURN
477           END SUBROUTINE DES_NEW_PARTICLE_TEST
478