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