File: N:\mfix\model\des\layout_mi_dem.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: SET_BC_DEM                                              !
4     !  Author: J.Musser                                   Date: 13-Jul-09  !
5     !                                                                      !
6     !  Purpose: Check the data provided for the des mass inflow boundary   !
7     !  condition and flag errors if the data is improper.  This module is  !
8     !  also used to convert the proveded information into the format       !
9     !  necessary for the dependent subrountines to function properly.      !
10     !                                                                      !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
12           SUBROUTINE LAYOUT_MI_DEM(BCV, BCV_I, MAX_DIA)
13     
14           use bc, only: BC_PLANE
15           USE run, only: RUN_TYPE
16     
17           use error_manager
18     
19           IMPLICIT NONE
20     
21           INTEGER, INTENT(IN) :: BCV
22           INTEGER, INTENT(IN) :: BCV_I      ! BC loop counter
23           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
24     ! Local debug flags.
25           LOGICAL, parameter :: setDBG = .FALSE.
26           LOGICAL, parameter :: showMAP = .FALSE.
27     
28           CALL INIT_ERR_MSG("LAYOUT_MI_DEM")
29     
30     ! This subroutine determines the pattern that the particles will need to
31     ! enter the system, if any. This routine only needs to be called if a
32     ! run is new.  If a run is a RESTART_1, all of the setup information
33     ! provided by this subroutine is will be obtained from the *_DES.RES file.
34     ! This is done due to this routine's strong dependence on the
35     ! RANDOM_NUMBER() subroutine.
36           IF(RUN_TYPE == 'RESTART_1') THEN
37              CALL SET_DEM_MI_OWNER(BCV, BCV_I)
38           ELSE
39              SELECT CASE (BC_PLANE(BCV))
40              CASE('N','S')
41                 CALL LAYOUT_DEM_MI_NS(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
42              CASE('E','W')
43                 CALL LAYOUT_DEM_MI_EW(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
44              CASE('T','B')
45                 CALL LAYOUT_DEM_MI_TB(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
46              END SELECT
47           ENDIF
48     
49           CALL FINL_ERR_MSG
50     
51           RETURN
52           END SUBROUTINE LAYOUT_MI_DEM
53     
54     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
55     !                                                                      !
56     !  Subroutine: LAYOUT_DEM_MI_NS                                        !
57     !                                                                      !
58     !  Purpose:  This routine determines the layout of the mass inlet as   !
59     !  either ordered or random based upon the inlet conditions.  This     !
60     !  routine also verifies that the specified inlet conditions for mass  !
61     !  or volumetric flow rates along with inlet size (length or area) and !
62     !  particle inlet velocity will work.  If not an error is flagged and  !
63     !  the program is exited.                                              !
64     !                                                                      !
65     !                                                                      !
66     !  Author: J.Musser                                   Date: 14-Aug-09  !
67     !                                                                      !
68     !  Comments:                                                           !
69     !                                                                      !
70     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
71           SUBROUTINE LAYOUT_DEM_MI_NS(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
72     
73     ! Global variables:
74     !---------------------------------------------------------------------//
75           use bc, only: BC_PLANE, BC_Y_s, BC_AREA
76           use bc, only: BC_X_w, BC_X_e
77           use bc, only: BC_Z_b, BC_Z_t
78     
79           use des_bc, only: DEM_MI
80     
81           use stl, only: N_FACETS_DES
82           use stl, only: STL_START, DEFAULT_STL
83           use stl, only: VERTEX, NORM_FACE
84           use cutcell, only: USE_STL
85     
86           use compar, only: myPE
87           use geometry, only: IMAX, JMAX, KMAX
88           use geometry, only: DX, DY, DZ
89           use geometry, only: XMIN, DO_K
90     
91           use funits, only: DMP_LOG
92     
93     ! Global parameters:
94     !---------------------------------------------------------------------//
95           use param1, only: ZERO, HALF, ONE
96     
97     ! Module procedures
98     !---------------------------------------------------------------------//
99           use des_bc, only: EXCLUDE_DEM_MI_CELL
100           use mpi_utility, only: GLOBAL_ALL_SUM
101           use mpi_utility, only: GLOBAL_ALL_MAX
102           use mpi_utility, only: GLOBAL_ALL_MIN
103           use stl_functions_des, only: TRI_BOX_OVERLAP
104           use functions, only: IS_ON_myPE_OWNS
105     
106           use error_manager
107     
108           IMPLICIT NONE
109     
110     ! Dummy arguments
111     !---------------------------------------------------------------------//
112     ! passed arguments giving index information of the boundary
113           INTEGER, INTENT(IN) :: BCV
114           INTEGER, INTENT(IN) :: BCV_I
115     ! Max diameter of incoming particles at bc
116           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
117     ! Debug flas.
118           LOGICAL, INTENT(IN) :: setDBG, showMAP
119     
120     ! Local variables
121     !---------------------------------------------------------------------//
122     ! Loop counters.
123           INTEGER LL, LC
124     ! Indices for mapping to fluid grid.
125           INTEGER I, J, K
126     ! Local MESH indices
127           INTEGER H, W
128     ! Temp variable for double precision
129           DOUBLE PRECISION :: TMP_DP
130     ! Temporary variable for integers
131           INTEGER :: TMP_INT
132     ! Window indicies.
133           INTEGER, allocatable :: MESH_H(:)
134           INTEGER, allocatable :: MESH_W(:)
135     ! Window positions.
136           DOUBLE PRECISION, allocatable :: MESH_P(:)
137           DOUBLE PRECISION, allocatable :: MESH_Q(:)
138     ! Map of array index to MI cell
139           INTEGER, allocatable :: RAND_MAP(:)
140     ! Map of MI cells to owner processes
141           INTEGER, allocatable :: FULL_MAP(:,:)
142     ! max number of partitions along length of inlet
143           INTEGER :: WMAX, HMAX
144           INTEGER :: maxEXT(2), minEXT(2)
145     ! the length of each side of the inlet boundary
146           DOUBLE PRECISION :: PLEN, QLEN
147     ! Number of occupied mesh cells
148           INTEGER :: OCCUPANTS
149     ! Offset and window size.
150           DOUBLE PRECISION :: SHIFT, WINDOW
151     ! The origin and dimension of MI cells. (STL intersection tests)
152           DOUBLE PRECISION :: CENTER(3), HALFSIZE(3)
153     ! Indicates that a separating axis exists
154           LOGICAL :: OVERLAP
155     ! Debug flag.
156           LOGICAL :: dFlag
157     !......................................................................!
158     
159     ! Initialize the error manager.
160           CALL INIT_ERR_MSG('LAYOUT_DEM_MI_NS')
161     
162           dFlag = (DMP_LOG .AND. setDBG)
163           if(dFlag) write(*,"(2/,'Building NS DEM_MI: ',I3)") BCV_I
164     
165     ! Store the index that maps back to the user input.
166     
167           OCCUPANTS = 0
168     
169     ! Calculate number of partitions in first direction.
170           PLEN = BC_X_e(BCV) - BC_X_w(BCV)
171           WMAX = FLOOR(real(PLEN/MAX_DIA))
172           allocate( MESH_W(WMAX) )
173           allocate( MESH_P(0:WMAX) )
174     
175     ! Calculate number of partitions in the second direction.
176           QLEN = merge(BC_Z_t(BCV) - BC_Z_b(BCV), MAX_DIA, DO_K)
177           HMAX = FLOOR(real(QLEN/MAX_DIA))
178           allocate( MESH_H(HMAX) )
179           allocate( MESH_Q(0:HMAX) )
180     
181     ! Allocate the full map.
182           allocate( FULL_MAP(WMAX, HMAX))
183     
184     ! Set the value of the boundary condtion offset value used in the
185     ! placement of new particles.
186           CALL CALC_CELL_INTERSECT(ZERO, BC_Y_s(BCV), DY, JMAX, J)
187           SHIFT = merge(-ONE, ONE, BC_PLANE(BCV) == 'N')
188           DEM_MI(BCV_I)%OFFSET = BC_Y_s(BCV) + MAX_DIA*SHIFT
189           DEM_MI(BCV_I)%L = J + int(SHIFT)
190           if(dFlag) write(*,"(2x,'Offset: ',3x,I4,3x,g12.5)") &
191              DEM_MI(BCV_I)%L, DEM_MI(BCV_I)%OFFSET
192     
193     
194     ! Dimension of grid cell for seeding particles; this may be larger than
195     ! than the particle diameter but not smaller:
196           DEM_MI(BCV_I)%WINDOW = MIN(PLEN/WMAX, QLEN/HMAX)
197           WINDOW = DEM_MI(BCV_I)%WINDOW
198           if(dFlag) write(*,"(2x,'Windows size: ',g12.5)") WINDOW
199     
200     ! Setup the first direction.
201           SHIFT = HALF*(PLEN - WMAX*WINDOW)
202           MESH_P(0) = BC_X_w(BCV) + SHIFT
203           if(dFlag) write(*,8005) 'P', SHIFT, 'P', MESH_P(0)
204           DO LC=1,WMAX
205              MESH_P(LC) = MESH_P(0) + dble(LC-1)*WINDOW
206              SHIFT = MESH_P(LC) + HALF*WINDOW
207              CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, MESH_W(LC))
208              IF(dFlag)WRITE(*,8006) LC, 'W', MESH_W(LC), 'P', MESH_P(LC)
209           ENDDO
210     
211     ! Setup the second direction.
212           IF(DO_K) THEN
213              SHIFT = HALF*(QLEN - HMAX*WINDOW)
214              MESH_Q(0) = BC_Z_b(BCV) + SHIFT
215              if(dFlag) write(*,8005) 'Q',SHIFT, 'Q',MESH_Q(0)
216              DO LC=1,HMAX
217                 MESH_Q(LC) = MESH_Q(0) + dble(LC-1)*WINDOW
218                 SHIFT = MESH_Q(LC) + HALF*WINDOW
219                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, MESH_H(LC))
220                 IF(dFlag)WRITE(*,8006) LC, 'H', MESH_H(LC), 'Q', MESH_Q(LC)
221              ENDDO
222           ELSE
223              MESH_H(1) = 1
224              MESH_Q(1) = 0.0d0
225           ENDIF
226     
227     ! Get the Jth index of the fluid cell
228           CALL CALC_CELL_INTERSECT(ZERO, BC_Y_s(BCV), DY, JMAX, J)
229     
230     
231     ! If the computationsl cell adjacent to the DEM_MI mesh cell is a
232     ! fluid cell and has not been cut, store the ID of the cell owner.
233           DO H=1,HMAX
234           DO W=1,WMAX
235     
236              I = MESH_W(W)
237              K = MESH_H(H)
238              FULL_MAP(W,H) = 0
239     
240              IF(.NOT.IS_ON_myPE_owns(I,J,K)) CYCLE
241     
242              IF(DO_K) THEN
243     
244                 CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
245                 CALL CALC_CELL_INTERSECT(ZERO, MESH_Q(H), DZ, KMAX, K)
246                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
247     
248                 SHIFT = MESH_P(W)+WINDOW
249                 CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, I)
250                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
251     
252                 SHIFT = MESH_Q(H)+WINDOW
253                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, K)
254                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
255     
256                 CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
257                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
258     
259              ELSE
260     
261                 K = MESH_H(1)
262                 CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
263                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
264     
265                 SHIFT = MESH_P(W)+WINDOW
266                 CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, I)
267                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
268              ENDIF
269     
270              FULL_MAP(W,H) = myPE+1
271           ENDDO
272           ENDDO
273     
274     
275     ! For complex boundaries defined by STLs, exclude any mass inflow cells
276     ! that intersect the boundary and all cells opposite the normal. The MI
277     ! cell sizes are increased by 10% to provide a small buffer.
278           IF(USE_STL) THEN
279     
280              HALFSIZE(2) = 1.10d0*MAX_DIA
281              HALFSIZE(1) = HALF*(WINDOW * 1.10d0)
282              HALFSIZE(3) = HALF*(WINDOW * 1.10d0)
283     
284              minEXT(1) = HMAX+1; maxEXT(1) = 0
285              minEXT(2) = WMAX+1; maxEXT(2) = 0
286              DO H=1,HMAX
287              DO W=1,WMAX
288     
289                 CENTER(2) = BC_Y_s(BCV)
290                 CENTER(1) = MESH_P(W) + HALF*WINDOW
291                 CENTER(3) = MESH_Q(H) + HALF*WINDOW
292     
293                 FACET_LP: DO LC=1, STL_START(DEFAULT_STL)-1
294     
295                    IF(BC_Y_s(BCV) > maxval(VERTEX(:,2,LC))) CYCLE FACET_LP
296                    IF(BC_Y_s(BCV) < minval(VERTEX(:,2,LC))) CYCLE FACET_LP
297     
298                    IF(BC_X_w(BCV) > maxval(VERTEX(:,1,LC))) CYCLE FACET_LP
299                    IF(BC_X_e(BCV) < minval(VERTEX(:,1,LC))) CYCLE FACET_LP
300     
301                    IF(BC_Z_b(BCV) > maxval(VERTEX(:,3,LC))) CYCLE FACET_LP
302                    IF(BC_Z_t(BCV) < minval(VERTEX(:,3,LC))) CYCLE FACET_LP
303     
304                    CALL TRI_BOX_OVERLAP(CENTER, HALFSIZE, &
305                       VERTEX(:,:,LC), OVERLAP)
306     
307                    IF(OVERLAP) THEN
308                       IF(NORM_FACE(1,LC) >= 0) THEN
309                          FULL_MAP(1:W,H) = 0
310                       ELSE
311                          FULL_MAP(W:WMAX,H) = 0
312                       ENDIF
313                       IF(NORM_FACE(3,LC) >= 0) THEN
314                          FULL_MAP(W,1:H) = 0
315                       ELSE
316                          FULL_MAP(W,H:HMAX) = 0
317                       ENDIF
318                       minEXT(1) = min(minEXT(1),H)
319                       minEXT(2) = min(minEXT(2),W)
320     
321                       maxEXT(1) = max(maxEXT(1),H)
322                       maxEXT(2) = max(maxEXT(2),W)
323                    ENDIF
324                 ENDDO FACET_LP
325              ENDDO
326              ENDDO
327              CALL GLOBAL_ALL_MIN(minEXT)
328              CALL GLOBAL_ALL_MAX(maxEXT)
329     
330              if(minEXT(1) < HMAX+1) FULL_MAP(:,:minEXT(1)) = 0
331              if(maxEXT(1) > 0) FULL_MAP(:,maxEXT(1):) = 0
332     
333              if(minEXT(2) /= WMAX+1) FULL_MAP(:minEXT(2),:) = 0
334              if(maxEXT(2) > 0) FULL_MAP(maxEXT(2):,:) = 0
335           ENDIF
336     
337     ! Add up the total number of available positions in the seeding map.
338           DO H=1,HMAX
339           DO W=1,WMAX
340              IF(FULL_MAP(W,H) /= 0) OCCUPANTS = OCCUPANTS + 1
341           ENDDO
342           ENDDO
343     
344     ! Sync the full map across all ranks.
345           CALL GLOBAL_ALL_SUM(OCCUPANTS)
346           CALL GLOBAL_ALL_SUM(FULL_MAP)
347     
348     ! Throw an error and exit if there are no occupants.
349           IF(OCCUPANTS == 0) THEN
350              WRITE(ERR_MSG, 1100) BCV_I
351              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
352           ENDIF
353     
354      1100 FORMAT('Error 1100: No un-cut fluid cells adjacent to DEM_MI ',  &
355              'staging area.',/'Unable to setup the discrete solids mass ', &
356              'inlet for BC:',I3)
357     
358     ! Store the number of occupants.
359           DEM_MI(BCV_I)%OCCUPANTS = OCCUPANTS
360     
361     ! Display the fill map if debugging
362           IF(dFlag .OR. (DMP_LOG .AND. showMAP)) THEN
363              WRITE(*,"(2/,2x,'Displaying Fill Map:')")
364              DO H=HMAX,1,-1
365                 WRITE(*,"(2x,'H =',I3)",advance='no')H
366                 DO W=1,WMAX
367                    IF(FULL_MAP(W,H) == 0) then
368                       WRITE(*,"(' *')",advance='no')
369                    ELSE
370                       WRITE(*,"(' .')",advance='no')
371                    ENDIF
372                 ENDDO
373                 WRITE(*,*)' '
374              ENDDO
375           ENDIF
376     
377     ! Construct an array of integers randomly ordered.
378           if(dFLAG) write(*,"(2/,2x,'Building RAND_MAP:')")
379           allocate( RAND_MAP(OCCUPANTS) )
380           RAND_MAP = 0
381     
382     ! Only Rank 0 will randomize the layout and broadcast it to the other
383     ! ranks. This will ensure that all ranks see the same layout.
384           IF(myPE == 0) THEN
385              LL = 1
386              DO WHILE (RAND_MAP(OCCUPANTS) .EQ. 0)
387                 CALL RANDOM_NUMBER(TMP_DP)
388                 TMP_INT = CEILING(real(TMP_DP*dble(OCCUPANTS)))
389                 DO LC = 1, LL
390                   IF(TMP_INT .EQ. RAND_MAP(LC) )EXIT
391                   IF(LC .EQ. LL)THEN
392                      if(dFlag) WRITE(*,"(4x,'LC:',I9,' : ',I9)") LC, TMP_INT
393                      RAND_MAP(LC) = TMP_INT
394                      LL = LL + 1
395                   ENDIF
396                 ENDDO
397              ENDDO
398           ENDIF
399     
400           CALL GLOBAL_ALL_SUM(RAND_MAP)
401     
402     ! Initialize the vacancy pointer.
403           DEM_MI(BCV_I)%VACANCY = 1
404     
405     ! Allocate the mass inlet storage arrays.
406           allocate( DEM_MI(BCV_I)%W(OCCUPANTS) )
407           allocate( DEM_MI(BCV_I)%P(OCCUPANTS) )
408           allocate( DEM_MI(BCV_I)%H(OCCUPANTS) )
409           allocate( DEM_MI(BCV_I)%Q(OCCUPANTS) )
410           allocate( DEM_MI(BCV_I)%OWNER(OCCUPANTS) )
411     
412           if(dFlag) write(*,8010)
413     ! Store the MI layout in the randomized order.
414           LC = 0
415           DO H=1,HMAX
416           DO W=1,WMAX
417              IF(FULL_MAP(W,H) == 0) CYCLE
418              LC = LC + 1
419              LL = RAND_MAP(LC)
420              DEM_MI(BCV_I)%OWNER(LL) = FULL_MAP(W,H) - 1
421     
422              DEM_MI(BCV_I)%W(LL) = MESH_W(W)
423              DEM_MI(BCV_I)%H(LL) = MESH_H(H)
424     
425              DEM_MI(BCV_I)%P(LL) = MESH_P(W)
426              DEM_MI(BCV_I)%Q(LL) = MESH_Q(H)
427     
428              if(dFlag) write(*,8011) DEM_MI(BCV_I)%OWNER(LL), &
429                 DEM_MI(BCV_I)%W(LL), DEM_MI(BCV_I)%H(LL), DEM_MI(BCV_I)%L, &
430                 DEM_MI(BCV_I)%P(LL), DEM_MI(BCV_I)%Q(LL), DEM_MI(BCV_I)%OFFSET
431     
432           ENDDO
433           ENDDO
434     
435     
436      8010 FORMAT(2/,2x,'Storing DEM_MI data:',/4X,'OWNER',5X,'W',5X,'H',   &
437              5X,'L',7X,'P',12X,'Q',12X,'R')
438      8011 FORMAT(4x,I5,3(2X,I4),3(2x,g12.5))
439     
440     
441           if(dFlag .OR. (DMP_LOG .AND. showMAP)) THEN
442              write(*,"(2/,2x,'Inlet area sizes:')")
443              write(*,9000) 'mfix.dat: ', PLEN * QLEN
444              write(*,9000) 'BC_AREA:  ', BC_AREA(BCV)
445              write(*,9000) 'DEM_MI:   ', OCCUPANTS * (WINDOW**2)
446           endif
447      9000 FORMAT(2x,A,g12.5)
448     
449     ! House keeping.
450           IF( allocated(MESH_H)) deallocate(MESH_H)
451           IF( allocated(MESH_W)) deallocate(MESH_W)
452           IF( allocated(MESH_P)) deallocate(MESH_P)
453           IF( allocated(MESH_Q)) deallocate(MESH_Q)
454     
455           IF( allocated(RAND_MAP)) deallocate(RAND_MAP)
456           IF( allocated(FULL_MAP)) deallocate(FULL_MAP)
457     
458           CALL FINL_ERR_MSG
459     
460           RETURN
461     
462      8005 FORMAT(2/,2x,'Building MESH_',A1,':',/4x,'Shift:',f8.4,/4x,      &
463              'MESH_',A1,'(0) = ',f8.4,/)
464     
465      8006 FORMAT(4x,'LC = ',I4,3x,A1,' =',I3,3x,A1,' =',f8.4)
466     
467           END SUBROUTINE LAYOUT_DEM_MI_NS
468     
469     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
470     !                                                                      !
471     !  Subroutine: LAYOUT_DEM_MI_EW                                        !
472     !                                                                      !
473     !  Purpose:  This routine determines the layout of the mass inlet as   !
474     !  either ordered or random based upon the inlet conditions.  This     !
475     !  routine also verifies that the specified inlet conditions for mass  !
476     !  or volumetric flow rates along with inlet size (length or area) and !
477     !  particle inlet velocity will work.  If not an error is flagged and  !
478     !  the program is exited.                                              !
479     !                                                                      !
480     !                                                                      !
481     !  Author: J.Musser                                   Date: 14-Aug-09  !
482     !                                                                      !
483     !  Comments:                                                           !
484     !                                                                      !
485     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
486           SUBROUTINE LAYOUT_DEM_MI_EW(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
487     
488           use bc, only: BC_PLANE, BC_X_w, BC_AREA
489           use bc, only: BC_Y_s, BC_Y_n
490           use bc, only: BC_Z_b, BC_Z_t
491     
492           use des_bc, only: DEM_MI
493     
494           use stl, only: STL_START, DEFAULT_STL
495           use stl, only: N_FACETS_DES
496           use stl, only: VERTEX, NORM_FACE
497           use cutcell, only: USE_STL
498     
499           use compar, only: myPE
500           use geometry, only: IMAX, JMAX, KMAX
501           use geometry, only: DX, DY, DZ
502           use geometry, only: XMIN, DO_K
503     
504           use funits, only: DMP_LOG
505     
506     ! Global parameters:
507     !---------------------------------------------------------------------//
508           use param1, only: ZERO, HALF, ONE
509     
510     ! Module procedures
511     !---------------------------------------------------------------------//
512           use des_bc, only: EXCLUDE_DEM_MI_CELL
513           use stl_functions_des, only: TRI_BOX_OVERLAP
514           use mpi_utility, only: GLOBAL_ALL_SUM
515           use mpi_utility, only: GLOBAL_ALL_MAX
516           use mpi_utility, only: GLOBAL_ALL_MIN
517           use functions, only: IS_ON_myPE_OWNS
518     
519           use error_manager
520     
521           IMPLICIT NONE
522     
523     ! Dummy arguments
524     !---------------------------------------------------------------------//
525     ! passed arguments giving index information of the boundary
526           INTEGER, INTENT(IN) :: BCV
527           INTEGER, INTENT(IN) :: BCV_I
528     ! Max diameter of incoming particles at bc
529           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
530     ! Debug flas.
531           LOGICAL, INTENT(IN) :: setDBG, showMAP
532     
533     ! Local variables
534     !---------------------------------------------------------------------//
535     ! Loop counters.
536           INTEGER LL, LC
537     ! Indices for mapping to fluid grid.
538           INTEGER I, J, K
539     ! Local MESH indices
540           INTEGER H, W
541     ! Temp variable for double precision
542           DOUBLE PRECISION :: TMP_DP
543     ! Temporary variable for integers
544           INTEGER :: TMP_INT
545     ! Window indicies.
546           INTEGER, allocatable :: MESH_H(:)
547           INTEGER, allocatable :: MESH_W(:)
548     ! Window positions.
549           DOUBLE PRECISION, allocatable :: MESH_P(:)
550           DOUBLE PRECISION, allocatable :: MESH_Q(:)
551     ! Map of array index to MI cell
552           INTEGER, allocatable :: RAND_MAP(:)
553     ! Map of MI cells to owner processes
554           INTEGER, allocatable :: FULL_MAP(:,:)
555     ! max number of partitions along length of inlet
556           INTEGER :: WMAX, HMAX
557           INTEGER :: maxEXT(2), minEXT(2)
558     ! the length of each side of the inlet boundary
559           DOUBLE PRECISION :: PLEN, QLEN
560     ! Number of occupied mesh cells
561           INTEGER :: OCCUPANTS
562     ! Offset and window size.
563           DOUBLE PRECISION :: SHIFT, WINDOW
564     ! The origin and dimension of MI cells. (STL intersection tests)
565           DOUBLE PRECISION :: CENTER(3), HALFSIZE(3)
566     ! Separating axis test dummy variable
567           INTEGER :: SEP_AXIS
568     ! Indicates that a separating axis exists
569           LOGICAL :: OVERLAP
570     ! Local debug flag.
571           LOGICAL :: dFlag
572     !......................................................................!
573     
574     
575     ! Initialize the error manager.
576           CALL INIT_ERR_MSG('LAYOUT_DEM_MI_EW')
577     
578           dFlag = (DMP_LOG .AND. setDBG)
579           if(dFlag) write(*,"(2/,'Building EW DEM_MI: ',I3)") BCV_I
580     
581           OCCUPANTS = 0
582     
583     ! Calculate number of partitions in first direction.
584           PLEN = BC_Y_n(BCV) - BC_Y_s(BCV)
585           WMAX = FLOOR(real(PLEN/MAX_DIA))
586           allocate( MESH_W(WMAX) )
587           allocate( MESH_P(0:WMAX) )
588     
589     ! Calculate number of partitions in the second direction.
590           QLEN = merge(BC_Z_t(BCV) - BC_Z_b(BCV), MAX_DIA, DO_K)
591           HMAX = FLOOR(real(QLEN/MAX_DIA))
592           allocate( MESH_H(HMAX) )
593           allocate( MESH_Q(0:HMAX) )
594     
595     ! Allocate the full map.
596           allocate( FULL_MAP(WMAX, HMAX))
597     
598     ! Set the value of the boundary condtion offset value used in the
599     ! placement of new particles.
600           CALL CALC_CELL_INTERSECT(XMIN, BC_X_w(BCV), DX, IMAX, I)
601           SHIFT = merge(-ONE, ONE, BC_PLANE(BCV) == 'E')
602           DEM_MI(BCV_I)%OFFSET = BC_X_w(BCV) + MAX_DIA*SHIFT
603           DEM_MI(BCV_I)%L = I + int(SHIFT)
604           if(dFlag) write(*,"(2x,'Offset: ',3x,I4,3x,g12.5)") &
605              DEM_MI(BCV_I)%L, DEM_MI(BCV_I)%OFFSET
606     
607     
608     ! Dimension of grid cell for seeding particles; this may be larger than
609     ! than the particle diameter but not smaller:
610           DEM_MI(BCV_I)%WINDOW = MIN(PLEN/WMAX, QLEN/HMAX)
611           WINDOW = DEM_MI(BCV_I)%WINDOW
612           if(dFlag) write(*,"(2x,'Windows size: ',g12.5)") WINDOW
613     
614     ! Setup the first direction.
615           SHIFT = HALF*(PLEN - WMAX*WINDOW)
616           MESH_P(0) = BC_Y_s(BCV) + SHIFT
617           if(dFlag) write(*,8005) 'P', SHIFT, 'P', MESH_P(0)
618           DO LC=1,WMAX
619              MESH_P(LC) = MESH_P(0) + dble(LC-1)*WINDOW
620              SHIFT = MESH_P(LC) + HALF*WINDOW
621              CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, MESH_W(LC))
622              IF(dFlag)WRITE(*,8006) LC, 'W', MESH_W(LC), 'P', MESH_P(LC)
623           ENDDO
624     
625     ! Setup the second direction.
626           IF(DO_K) THEN
627              SHIFT = HALF*(QLEN - HMAX*WINDOW)
628              MESH_Q(0) = BC_Z_b(BCV) + SHIFT
629              if(dFlag) write(*,8005) 'Q',SHIFT, 'Q',MESH_Q(0)
630              DO LC=1,HMAX
631                 MESH_Q(LC) = MESH_Q(0) + dble(LC-1)*WINDOW
632                 SHIFT = MESH_Q(LC) + HALF*WINDOW
633                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, MESH_H(LC))
634                 IF(dFlag)WRITE(*,8006) LC, 'H', MESH_H(LC), 'Q', MESH_Q(LC)
635              ENDDO
636           ELSE
637              MESH_H(1) = 1
638              MESH_Q(1) = 0.0d0
639           ENDIF
640     
641     ! Get the Jth index of the fluid cell
642           CALL CALC_CELL_INTERSECT(XMIN, BC_X_w(BCV), DX, IMAX, I)
643     
644     
645     ! If the computationsl cell adjacent to the DEM_MI mesh cell is a
646     ! fluid cell, store the rank of the cell owner.
647           DO H=1,HMAX
648           DO W=1,WMAX
649     
650              J = MESH_W(W)
651              K = MESH_H(H)
652              FULL_MAP(W,H) = 0
653     
654              IF(.NOT.IS_ON_myPE_owns(I,J,K)) CYCLE
655     
656              IF(DO_K) THEN
657     
658                 CALL CALC_CELL_INTERSECT(ZERO, MESH_P(W), DY, JMAX, J)
659                 CALL CALC_CELL_INTERSECT(ZERO, MESH_Q(H), DZ, KMAX, K)
660                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
661     
662                 SHIFT = MESH_P(W)+WINDOW
663                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, J)
664                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
665     
666                 SHIFT = MESH_Q(H)+WINDOW
667                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DZ, KMAX, K)
668                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
669     
670                 CALL CALC_CELL_INTERSECT(ZERO, MESH_P(W), DY, JMAX, J)
671                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
672     
673              ELSE
674     
675                 K = MESH_H(1)
676                 CALL CALC_CELL_INTERSECT(ZERO, MESH_P(W), DY, JMAX, J)
677                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
678     
679                 SHIFT = MESH_P(W)+WINDOW
680                 CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, J)
681                 IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
682              ENDIF
683     
684              FULL_MAP(W,H) = myPE+1
685           ENDDO
686           ENDDO
687     
688     
689     ! For complex boundaries defined by STLs, exclude any mass inflow cells
690     ! that intersect the boundary and all cells opposite the normal. The MI
691     ! cell sizes are increased by 10% to provide a small buffer.
692           IF(USE_STL) THEN
693     
694              HALFSIZE(1) = 1.10d0*MAX_DIA
695              HALFSIZE(2) = HALF*(WINDOW * 1.10d0)
696              HALFSIZE(3) = HALF*(WINDOW * 1.10d0)
697     
698              minEXT(1) = HMAX+1; maxEXT(1) = 0
699              minEXT(2) = WMAX+1; maxEXT(2) = 0
700     
701              DO H=1,HMAX
702              DO W=1,WMAX
703     
704                 CENTER(1) = BC_X_w(BCV)
705                 CENTER(2) = MESH_P(W) + HALF*WINDOW
706                 CENTER(3) = MESH_Q(H) + HALF*WINDOW
707     
708                 FACET_LP: DO LC=1, STL_START(DEFAULT_STL)-1
709     
710                    IF(BC_X_w(BCV) > maxval(VERTEX(:,1,LC))) CYCLE FACET_LP
711                    IF(BC_X_w(BCV) < minval(VERTEX(:,1,LC))) CYCLE FACET_LP
712     
713                    IF(BC_Y_s(BCV) > maxval(VERTEX(:,2,LC))) CYCLE FACET_LP
714                    IF(BC_Y_n(BCV) < minval(VERTEX(:,2,LC))) CYCLE FACET_LP
715     
716                    IF(BC_Z_b(BCV) > maxval(VERTEX(:,3,LC))) CYCLE FACET_LP
717                    IF(BC_Z_t(BCV) < minval(VERTEX(:,3,LC))) CYCLE FACET_LP
718     
719                    CALL TRI_BOX_OVERLAP(CENTER, HALFSIZE, &
720                       VERTEX(:,:,LC), OVERLAP)
721     
722                    IF(OVERLAP) THEN
723                       IF(NORM_FACE(2,LC) >= 0) THEN
724                          FULL_MAP(1:W,H) = 0
725                       ELSE
726                          FULL_MAP(W:WMAX,H) = 0
727                       ENDIF
728                       IF(NORM_FACE(3,LC) >= 0) THEN
729                          FULL_MAP(W,1:H) = 0
730                       ELSE
731                          FULL_MAP(W,H:HMAX) = 0
732                       ENDIF
733     
734                       minEXT(1) = min(minEXT(1),H)
735                       minEXT(2) = min(minEXT(2),W)
736     
737                       maxEXT(1) = max(maxEXT(1),H)
738                       maxEXT(2) = max(maxEXT(2),W)
739     
740                    ENDIF
741                 ENDDO FACET_LP
742              ENDDO
743              ENDDO
744     
745              CALL GLOBAL_ALL_MIN(minEXT)
746              CALL GLOBAL_ALL_MAX(maxEXT)
747     
748              if(minEXT(1) < HMAX+1) FULL_MAP(:,:minEXT(1)) = 0
749              if(maxEXT(1) > 0) FULL_MAP(:,maxEXT(1):) = 0
750     
751              if(minEXT(2) /= WMAX+1) FULL_MAP(:minEXT(2),:) = 0
752              if(maxEXT(2) > 0) FULL_MAP(maxEXT(2):,:) = 0
753     
754           ENDIF
755     
756     ! Add up the total number of available positions in the seeding map.
757           DO H=1,HMAX
758           DO W=1,WMAX
759              IF(FULL_MAP(W,H) /= 0) OCCUPANTS = OCCUPANTS + 1
760           ENDDO
761           ENDDO
762     
763     ! Sync the full map across all ranks.
764           CALL GLOBAL_ALL_SUM(OCCUPANTS)
765           CALL GLOBAL_ALL_SUM(FULL_MAP)
766     
767     ! Throw an error and exit if there are no occupants.
768           IF(OCCUPANTS == 0) THEN
769              WRITE(ERR_MSG, 1100) BCV_I
770              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
771           ENDIF
772     
773      1100 FORMAT('Error 1100: No un-cut fluid cells adjacent to DEM_MI ',  &
774              'staging area.',/'Unable to setup the discrete solids mass ', &
775              'inlet for BC:',I3)
776     
777     ! Store the number of occupants.
778           DEM_MI(BCV_I)%OCCUPANTS = OCCUPANTS
779     
780     ! Display the fill map if debugging
781           IF(dFlag .OR. (DMP_LOG .AND. showMAP)) THEN
782              WRITE(*,"(2/,2x,'Displaying Fill Map:')")
783              DO H=HMAX,1,-1
784                 WRITE(*,"(2x,'H =',I3)",advance='no')H
785                 DO W=1,WMAX
786                    IF(FULL_MAP(W,H) == 0) then
787                       WRITE(*,"(' *')",advance='no')
788                    ELSE
789                       WRITE(*,"(' .')",advance='no')
790                    ENDIF
791                 ENDDO
792                 WRITE(*,*)' '
793              ENDDO
794           ENDIF
795     
796     ! Construct an array of integers randomly ordered.
797           if(dFLAG) write(*,"(2/,2x,'Building RAND_MAP:')")
798           allocate( RAND_MAP(OCCUPANTS) )
799           RAND_MAP = 0
800     
801     ! Only Rank 0 will randomize the layout and broadcast it to the other
802     ! ranks. This will ensure that all ranks see the same layout.
803           IF(myPE == 0) THEN
804              LL = 1
805              DO WHILE (RAND_MAP(OCCUPANTS) .EQ. 0)
806                 CALL RANDOM_NUMBER(TMP_DP)
807                 TMP_INT = CEILING(real(TMP_DP*dble(OCCUPANTS)))
808                 DO LC = 1, LL
809                   IF(TMP_INT .EQ. RAND_MAP(LC) )EXIT
810                   IF(LC .EQ. LL)THEN
811                      if(dFlag) WRITE(*,"(4x,'LC:',I6,' : ',I6)") LC, TMP_INT
812                      RAND_MAP(LC) = TMP_INT
813                      LL = LL + 1
814                   ENDIF
815                 ENDDO
816              ENDDO
817           ENDIF
818     
819           CALL GLOBAL_ALL_SUM(RAND_MAP)
820     
821     ! Initialize the vacancy pointer.
822           DEM_MI(BCV_I)%VACANCY = 1
823     
824     ! Allocate the mass inlet storage arrays.
825           allocate( DEM_MI(BCV_I)%W(OCCUPANTS) )
826           allocate( DEM_MI(BCV_I)%P(OCCUPANTS) )
827           allocate( DEM_MI(BCV_I)%H(OCCUPANTS) )
828           allocate( DEM_MI(BCV_I)%Q(OCCUPANTS) )
829           allocate( DEM_MI(BCV_I)%OWNER(OCCUPANTS) )
830     
831           if(dFlag) write(*,8010)
832     ! Store the MI layout in the randomized order.
833           LC = 0
834           DO H=1,HMAX
835           DO W=1,WMAX
836              IF(FULL_MAP(W,H) == 0) CYCLE
837              LC = LC + 1
838              LL = RAND_MAP(LC)
839              DEM_MI(BCV_I)%OWNER(LL) = FULL_MAP(W,H) - 1
840     
841              DEM_MI(BCV_I)%W(LL) = MESH_W(W)
842              DEM_MI(BCV_I)%H(LL) = MESH_H(H)
843     
844              DEM_MI(BCV_I)%P(LL) = MESH_P(W)
845              DEM_MI(BCV_I)%Q(LL) = MESH_Q(H)
846     
847              if(dFlag) write(*,8011) DEM_MI(BCV_I)%OWNER(LL), &
848                 DEM_MI(BCV_I)%W(LL), DEM_MI(BCV_I)%H(LL), DEM_MI(BCV_I)%L, &
849                 DEM_MI(BCV_I)%P(LL), DEM_MI(BCV_I)%Q(LL), DEM_MI(BCV_I)%OFFSET
850     
851           ENDDO
852           ENDDO
853     
854     
855      8010 FORMAT(2/,2x,'Storing DEM_MI data:',/4X,'OWNER',5X,'W',5X,'H',   &
856              5X,'L',7X,'P',12X,'Q',12X,'R')
857      8011 FORMAT(4x,I5,3(2X,I4),3(2x,g12.5))
858     
859           if(dFlag .OR. (DMP_LOG .AND. showMAP)) THEN
860              write(*,"(2/,2x,'Inlet area sizes:')")
861              write(*,9000) 'mfix.dat: ', PLEN * QLEN
862              write(*,9000) 'BC_AREA:  ', BC_AREA(BCV)
863              write(*,9000) 'DEM_MI:   ', OCCUPANTS * (WINDOW**2)
864           endif
865      9000 FORMAT(2x,A,g12.5)
866     
867     ! House keeping.
868           IF( allocated(MESH_H)) deallocate(MESH_H)
869           IF( allocated(MESH_W)) deallocate(MESH_W)
870           IF( allocated(MESH_P)) deallocate(MESH_P)
871           IF( allocated(MESH_Q)) deallocate(MESH_Q)
872     
873           IF( allocated(RAND_MAP)) deallocate(RAND_MAP)
874           IF( allocated(FULL_MAP)) deallocate(FULL_MAP)
875     
876           CALL FINL_ERR_MSG
877     
878           RETURN
879     
880      8005 FORMAT(2/,2x,'Building MESH_',A1,':',/4x,'Shift:',f8.4,/4x,      &
881              'MESH_',A1,'(0) = ',f8.4,/)
882     
883      8006 FORMAT(4x,'LC = ',I4,3x,A1,' =',I3,3x,A1,' =',f8.4)
884     
885           END SUBROUTINE LAYOUT_DEM_MI_EW
886     
887     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
888     !                                                                      !
889     !  Subroutine: LAYOUT_DEM_MI_TB                                        !
890     !                                                                      !
891     !  Purpose:  This routine determines the layout of the mass inlet as   !
892     !  either ordered or random based upon the inlet conditions.  This     !
893     !  routine also verifies that the specified inlet conditions for mass  !
894     !  or volumetric flow rates along with inlet size (length or area) and !
895     !  particle inlet velocity will work.  If not an error is flagged and  !
896     !  the program is exited.                                              !
897     !                                                                      !
898     !                                                                      !
899     !  Author: J.Musser                                   Date: 14-Aug-09  !
900     !                                                                      !
901     !  Comments:                                                           !
902     !                                                                      !
903     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
904           SUBROUTINE LAYOUT_DEM_MI_TB(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
905     
906     ! Global variables:
907     !---------------------------------------------------------------------//
908           use bc, only: BC_PLANE, BC_Z_b, BC_AREA
909           use bc, only: BC_X_w, BC_X_e
910           use bc, only: BC_Y_s, BC_Y_n
911     
912           use des_bc, only: DEM_MI
913     
914           use stl, only: STL_START, DEFAULT_STL
915           use stl, only: N_FACETS_DES
916           use stl, only: VERTEX, NORM_FACE
917           use cutcell, only: USE_STL
918           use compar, only: myPE
919           use geometry, only: IMAX, JMAX, KMAX
920           use geometry, only: DX, DY, DZ
921           use geometry, only: XMIN
922     
923           use funits, only: DMP_LOG
924     
925     ! Global parameters:
926     !---------------------------------------------------------------------//
927           use param1, only: ZERO, HALF, ONE
928     
929     ! Module procedures
930     !---------------------------------------------------------------------//
931           use des_bc, only: EXCLUDE_DEM_MI_CELL
932           use mpi_utility, only: GLOBAL_ALL_SUM
933           use mpi_utility, only: GLOBAL_ALL_MAX
934           use mpi_utility, only: GLOBAL_ALL_MIN
935           use stl_functions_des, only: TRI_BOX_OVERLAP
936           use functions, only: IS_ON_myPE_OWNS
937     
938           use error_manager
939     
940           IMPLICIT NONE
941     
942     ! Dummy arguments
943     !---------------------------------------------------------------------//
944     ! passed arguments giving index information of the boundary
945           INTEGER, INTENT(IN) :: BCV
946           INTEGER, INTENT(IN) :: BCV_I
947     ! Max diameter of incoming particles at bc
948           DOUBLE PRECISION, INTENT(IN) :: MAX_DIA
949     ! Debug flas.
950           LOGICAL, INTENT(IN) :: setDBG, showMAP
951     
952     ! Local variables
953     !---------------------------------------------------------------------//
954     ! Loop counters.
955           INTEGER LL, LC
956     ! Indices for mapping to fluid grid.
957           INTEGER I, J, K
958     ! Local MESH indices
959           INTEGER H, W
960     ! Temp variable for double precision
961           DOUBLE PRECISION :: TMP_DP
962     ! Temporary variable for integers
963           INTEGER :: TMP_INT
964     ! Window indicies.
965           INTEGER, allocatable :: MESH_H(:)
966           INTEGER, allocatable :: MESH_W(:)
967     ! Window positions.
968           DOUBLE PRECISION, allocatable :: MESH_P(:)
969           DOUBLE PRECISION, allocatable :: MESH_Q(:)
970     ! Map of array index to MI cell
971           INTEGER, allocatable :: RAND_MAP(:)
972     ! Map of MI cells to owner processes
973           INTEGER, allocatable :: FULL_MAP(:,:)
974     ! max number of partitions along length of inlet
975           INTEGER :: WMAX, HMAX
976           INTEGER :: maxEXT(2), minEXT(2)
977     ! the length of each side of the inlet boundary
978           DOUBLE PRECISION :: PLEN, QLEN
979     ! Number of occupied mesh cells
980           INTEGER :: OCCUPANTS
981     ! Offset and and window size.
982           DOUBLE PRECISION :: SHIFT, WINDOW
983     ! The origin and dimension of MI cells. (STL intersection tests)
984           DOUBLE PRECISION :: CENTER(3), HALFSIZE(3)
985     ! Separating axis test dummy variable
986           INTEGER :: SEP_AXIS
987     ! Indicates that a separating axis exists
988           LOGICAL :: OVERLAP
989     ! Local Debug flag.
990           LOGICAL :: dFlag
991     
992     !......................................................................!
993     
994     ! Initialize the error manager.
995           CALL INIT_ERR_MSG('LAYOUT_DEM_MI_TB')
996     
997           dFlag = (DMP_LOG .AND. setDBG)
998           if(dFlag) write(*,"(2/,'Building TB DEM_MI: ',I3)") BCV_I
999     
1000     ! Store the index that maps back to the user input.
1001     
1002           OCCUPANTS = 0
1003     
1004     ! Calculate number of partitions in first direction.
1005           PLEN = BC_X_e(BCV) - BC_X_w(BCV)
1006           WMAX = FLOOR(real(PLEN/MAX_DIA))
1007           allocate( MESH_W(WMAX) )
1008           allocate( MESH_P(0:WMAX) )
1009     
1010     ! Calculate number of partitions in the second direction.
1011           QLEN = BC_Y_n(BCV) - BC_Y_s(BCV)
1012           HMAX = FLOOR(real(QLEN/MAX_DIA))
1013           allocate( MESH_H(HMAX) )
1014           allocate( MESH_Q(0:HMAX) )
1015     
1016     ! Allocate the full map.
1017           allocate( FULL_MAP(WMAX, HMAX))
1018     
1019     ! Set the value of the boundary condtion offset value used in the
1020     ! placement of new particles.
1021           CALL CALC_CELL_INTERSECT(ZERO, BC_Z_b(BCV), DZ, KMAX, K)
1022           SHIFT = merge(-ONE, ONE, BC_PLANE(BCV) == 'T')
1023           DEM_MI(BCV_I)%OFFSET = BC_Z_b(BCV) + MAX_DIA*SHIFT
1024           DEM_MI(BCV_I)%L = K + int(SHIFT)
1025           if(dFlag) write(*,"(2x,'Offset: ',3x,I4,3x,g12.5)") &
1026              DEM_MI(BCV_I)%L, DEM_MI(BCV_I)%OFFSET
1027     
1028     
1029     ! Dimension of grid cell for seeding particles; this may be larger than
1030     ! than the particle diameter but not smaller:
1031           DEM_MI(BCV_I)%WINDOW = MIN(PLEN/WMAX, QLEN/HMAX)
1032           WINDOW = DEM_MI(BCV_I)%WINDOW
1033           if(dFlag) write(*,"(2x,'Windows size: ',g12.5)") WINDOW
1034     
1035     ! Setup the first direction.
1036           SHIFT = HALF*(PLEN - WMAX*WINDOW)
1037           MESH_P(0) = BC_X_w(BCV) + SHIFT
1038           if(dFlag) write(*,8005) 'P', SHIFT, 'P', MESH_P(0)
1039           DO LC=1,WMAX
1040              MESH_P(LC) = MESH_P(0) + dble(LC-1)*WINDOW
1041              SHIFT = MESH_P(LC) + HALF*WINDOW
1042              CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, MESH_W(LC))
1043              IF(dFlag)WRITE(*,8006) LC, 'W', MESH_W(LC), 'P', MESH_P(LC)
1044           ENDDO
1045     
1046     ! Setup the second direction.
1047           SHIFT = HALF*(QLEN - HMAX*WINDOW)
1048           MESH_Q(0) = BC_Y_s(BCV) + SHIFT
1049           if(dFlag) write(*,8005) 'Q',SHIFT, 'Q',MESH_Q(0)
1050           DO LC=1,HMAX
1051              MESH_Q(LC) = MESH_Q(0) + dble(LC-1)*WINDOW
1052              SHIFT = MESH_Q(LC) + HALF*WINDOW
1053              CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, MESH_H(LC))
1054              IF(dFlag)WRITE(*,8006) LC, 'H', MESH_H(LC), 'Q', MESH_Q(LC)
1055           ENDDO
1056     
1057     ! Get the Jth index of the fluid cell
1058           CALL CALC_CELL_INTERSECT(ZERO, BC_Z_b(BCV), DZ, KMAX, K)
1059     
1060     ! If the computationsl cell adjacent to the DEM_MI mesh cell is a
1061     ! fluid cell and has not been cut, store the ID of the cell owner.
1062           DO H=1,HMAX
1063           DO W=1,WMAX
1064     
1065              I = MESH_W(W)
1066              J = MESH_H(H)
1067              FULL_MAP(W,H) = 0
1068     
1069              IF(.NOT.IS_ON_myPE_owns(I,J,K)) CYCLE
1070     
1071              CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
1072              CALL CALC_CELL_INTERSECT(ZERO, MESH_Q(H), DY, JMAX, J)
1073              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1074     
1075              SHIFT = MESH_P(W)+WINDOW
1076              CALL CALC_CELL_INTERSECT(XMIN, SHIFT, DX, IMAX, I)
1077              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1078     
1079              SHIFT = MESH_Q(H)+WINDOW
1080              CALL CALC_CELL_INTERSECT(ZERO, SHIFT, DY, JMAX, J)
1081              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1082     
1083              CALL CALC_CELL_INTERSECT(XMIN, MESH_P(W), DX, IMAX, I)
1084              IF(EXCLUDE_DEM_MI_CELL(I, J, K)) CYCLE
1085     
1086              FULL_MAP(W,H) = myPE+1
1087           ENDDO
1088           ENDDO
1089     
1090     ! For complex boundaries defined by STLs, exclude any mass inflow cells
1091     ! that intersect the boundary and all cells opposite the normal. The MI
1092     ! cell sizes are increased by 10% to provide a small buffer.
1093           IF(USE_STL) THEN
1094     
1095              HALFSIZE(3) = 1.10d0*MAX_DIA
1096              HALFSIZE(1) = HALF*(WINDOW * 1.10d0)
1097              HALFSIZE(2) = HALF*(WINDOW * 1.10d0)
1098     
1099              minEXT(1) = HMAX+1; maxEXT(1) = 0
1100              minEXT(2) = WMAX+1; maxEXT(2) = 0
1101     
1102              DO H=1,HMAX
1103              DO W=1,WMAX
1104     
1105                 CENTER(3) = BC_Z_b(BCV)
1106                 CENTER(1) = MESH_P(W) + HALF*WINDOW
1107                 CENTER(2) = MESH_Q(H) + HALF*WINDOW
1108     
1109                 FACET_LP: DO LC=1, STL_START(DEFAULT_STL)-1
1110     
1111                    IF(BC_Z_b(BCV) > maxval(VERTEX(:,3,LC))) CYCLE FACET_LP
1112                    IF(BC_Z_b(BCV) < minval(VERTEX(:,3,LC))) CYCLE FACET_LP
1113     
1114                    IF(BC_X_w(BCV) > maxval(VERTEX(:,1,LC))) CYCLE FACET_LP
1115                    IF(BC_X_e(BCV) < minval(VERTEX(:,1,LC))) CYCLE FACET_LP
1116     
1117                    IF(BC_Y_s(BCV) > maxval(VERTEX(:,2,LC))) CYCLE FACET_LP
1118                    IF(BC_Y_n(BCV) < minval(VERTEX(:,2,LC))) CYCLE FACET_LP
1119     
1120                    CALL TRI_BOX_OVERLAP(CENTER, HALFSIZE, &
1121                       VERTEX(:,:,LC), OVERLAP)
1122     
1123                    IF(OVERLAP) THEN
1124                       IF(NORM_FACE(1,LC) >= 0) THEN
1125                          FULL_MAP(1:W,H) = 0
1126                       ELSE
1127                          FULL_MAP(W:WMAX,H) = 0
1128                       ENDIF
1129                       IF(NORM_FACE(2,LC) >= 0) THEN
1130                          FULL_MAP(W,1:H) = 0
1131                       ELSE
1132                          FULL_MAP(W,H:HMAX) = 0
1133                       ENDIF
1134     
1135                       minEXT(1) = min(minEXT(1),H)
1136                       minEXT(2) = min(minEXT(2),W)
1137     
1138                       maxEXT(1) = max(maxEXT(1),H)
1139                       maxEXT(2) = max(maxEXT(2),W)
1140     
1141                    ENDIF
1142                 ENDDO FACET_LP
1143              ENDDO
1144              ENDDO
1145     
1146              CALL GLOBAL_ALL_MIN(minEXT)
1147              CALL GLOBAL_ALL_MAX(maxEXT)
1148     
1149              if(minEXT(1) < HMAX+1) FULL_MAP(:,:minEXT(1)) = 0
1150              if(maxEXT(1) > 0) FULL_MAP(:,maxEXT(1):) = 0
1151     
1152              if(minEXT(2) /= WMAX+1) FULL_MAP(:minEXT(2),:) = 0
1153              if(maxEXT(2) > 0) FULL_MAP(maxEXT(2):,:) = 0
1154     
1155           ENDIF
1156     
1157     ! Add up the total number of available positions in the seeding map.
1158           DO H=1,HMAX
1159           DO W=1,WMAX
1160              IF(FULL_MAP(W,H) /= 0) OCCUPANTS = OCCUPANTS + 1
1161           ENDDO
1162           ENDDO
1163     
1164     ! Sync the full map across all ranks.
1165           CALL GLOBAL_ALL_SUM(OCCUPANTS)
1166           CALL GLOBAL_ALL_SUM(FULL_MAP)
1167     
1168     ! Throw an error and exit if there are no occupants.
1169           IF(OCCUPANTS == 0) THEN
1170              WRITE(ERR_MSG, 1100) BCV_I
1171              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
1172           ENDIF
1173     
1174      1100 FORMAT('Error 1100: No un-cut fluid cells adjacent to DEM_MI ',  &
1175              'staging area.',/'Unable to setup the discrete solids mass ', &
1176              'inlet for BC:',I3)
1177     
1178     ! Store the number of occupants.
1179           DEM_MI(BCV_I)%OCCUPANTS = OCCUPANTS
1180     
1181     ! Display the fill map if debugging
1182           IF(dFlag .OR. (DMP_LOG .AND. showMAP)) THEN
1183              WRITE(*,"(2/,2x,'Displaying Fill Map:')")
1184              DO H=HMAX,1,-1
1185                 WRITE(*,"(2x,'H =',I3)",advance='no')H
1186                 DO W=1,WMAX
1187                    IF(FULL_MAP(W,H) == 0) then
1188                       WRITE(*,"(' *')",advance='no')
1189                    ELSE
1190                       WRITE(*,"(' .')",advance='no')
1191                    ENDIF
1192                 ENDDO
1193                 WRITE(*,*)' '
1194              ENDDO
1195           ENDIF
1196     
1197     ! Construct an array of integers randomly ordered.
1198           if(dFLAG) write(*,"(2/,2x,'Building RAND_MAP:')")
1199           allocate( RAND_MAP(OCCUPANTS) )
1200           RAND_MAP = 0
1201     
1202     ! Only Rank 0 will randomize the layout and broadcast it to the other
1203     ! ranks. This will ensure that all ranks see the same layout.
1204           IF(myPE == 0) THEN
1205              LL = 1
1206              DO WHILE (RAND_MAP(OCCUPANTS) .EQ. 0)
1207                 CALL RANDOM_NUMBER(TMP_DP)
1208                 TMP_INT = CEILING(real(TMP_DP*dble(OCCUPANTS)))
1209                 DO LC = 1, LL
1210                   IF(TMP_INT .EQ. RAND_MAP(LC) )EXIT
1211                   IF(LC .EQ. LL)THEN
1212                      if(dFlag) WRITE(*,"(4x,'LC:',I3,' : ',I3)") LC, TMP_INT
1213                      RAND_MAP(LC) = TMP_INT
1214                      LL = LL + 1
1215                   ENDIF
1216                 ENDDO
1217              ENDDO
1218           ENDIF
1219     
1220           CALL GLOBAL_ALL_SUM(RAND_MAP)
1221     
1222     ! Initialize the vacancy pointer.
1223           DEM_MI(BCV_I)%VACANCY = 1
1224     
1225     ! Allocate the mass inlet storage arrays.
1226           allocate( DEM_MI(BCV_I)%W(OCCUPANTS) )
1227           allocate( DEM_MI(BCV_I)%P(OCCUPANTS) )
1228           allocate( DEM_MI(BCV_I)%H(OCCUPANTS) )
1229           allocate( DEM_MI(BCV_I)%Q(OCCUPANTS) )
1230           allocate( DEM_MI(BCV_I)%OWNER(OCCUPANTS) )
1231     
1232           if(dFlag) write(*,8010)
1233     ! Store the MI layout in the randomized order.
1234           LC = 0
1235           DO H=1,HMAX
1236           DO W=1,WMAX
1237              IF(FULL_MAP(W,H) == 0) CYCLE
1238              LC = LC + 1
1239              LL = RAND_MAP(LC)
1240              DEM_MI(BCV_I)%OWNER(LL) = FULL_MAP(W,H) - 1
1241     
1242              DEM_MI(BCV_I)%W(LL) = MESH_W(W)
1243              DEM_MI(BCV_I)%H(LL) = MESH_H(H)
1244     
1245              DEM_MI(BCV_I)%P(LL) = MESH_P(W)
1246              DEM_MI(BCV_I)%Q(LL) = MESH_Q(H)
1247     
1248              if(dFlag) write(*,8011) DEM_MI(BCV_I)%OWNER(LL), &
1249                 DEM_MI(BCV_I)%W(LL), DEM_MI(BCV_I)%H(LL), DEM_MI(BCV_I)%L, &
1250                 DEM_MI(BCV_I)%P(LL), DEM_MI(BCV_I)%Q(LL), DEM_MI(BCV_I)%OFFSET
1251     
1252           ENDDO
1253           ENDDO
1254     
1255     
1256      8010 FORMAT(2/,2x,'Storing DEM_MI data:',/4X,'OWNER',5X,'W',5X,'H',   &
1257              5X,'L',7X,'P',12X,'Q',12X,'R')
1258      8011 FORMAT(4x,I5,3(2X,I4),3(2x,g12.5))
1259     
1260           if(dFlag .OR. (DMP_LOG .AND. showMAP)) THEN
1261              write(*,"(2/,2x,'Inlet area sizes:')")
1262              write(*,9000) 'mfix.dat: ', PLEN * QLEN
1263              write(*,9000) 'BC_AREA:  ', BC_AREA(BCV)
1264              write(*,9000) 'DEM_MI:   ', OCCUPANTS * (WINDOW**2)
1265           endif
1266      9000 FORMAT(2x,A,g12.5)
1267     
1268     ! House keeping.
1269           IF( allocated(MESH_H)) deallocate(MESH_H)
1270           IF( allocated(MESH_W)) deallocate(MESH_W)
1271           IF( allocated(MESH_P)) deallocate(MESH_P)
1272           IF( allocated(MESH_Q)) deallocate(MESH_Q)
1273     
1274           IF( allocated(RAND_MAP)) deallocate(RAND_MAP)
1275           IF( allocated(FULL_MAP)) deallocate(FULL_MAP)
1276     
1277           CALL FINL_ERR_MSG
1278     
1279           RETURN
1280     
1281      8005 FORMAT(2/,2x,'Building MESH_',A1,':',/4x,'Shift:',f8.4,/4x,      &
1282              'MESH_',A1,'(0) = ',f8.4,/)
1283     
1284      8006 FORMAT(4x,'LC = ',I4,3x,A1,' =',I3,3x,A1,' =',f8.4)
1285     
1286           END SUBROUTINE LAYOUT_DEM_MI_TB
1287     
1288     
1289     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1290     !                                                                      !
1291     !  Subroutine: SET_DEM_MI_OWNER                                        !
1292     !  Author: J.Musser                                   Date: 08-Aug-14  !
1293     !                                                                      !
1294     !  Purpose: Set the owner of the background mesh used for seeding new  !
1295     !  particles into the domain. The mesh is stored and read from the RES !
1296     !  file, however, the owers need to be set per-run in case the DMP     !
1297     !  partitions changed.                                                 !
1298     !                                                                      !
1299     !  Comments:                                                           !
1300     !                                                                      !
1301     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1302           SUBROUTINE SET_DEM_MI_OWNER(BCV, BCV_I)
1303     
1304           use bc, only: BC_PLANE
1305           use des_bc, only: DEM_MI
1306     
1307           use compar
1308           use geometry
1309           use indices
1310     
1311     ! Module procedures
1312     !---------------------------------------------------------------------//
1313           use mpi_utility, only: GLOBAL_ALL_SUM
1314           use error_manager
1315           use functions
1316     
1317           IMPLICIT NONE
1318     !-----------------------------------------------
1319     ! Dummy arguments
1320     !-----------------------------------------------
1321     ! passed arguments giving index information of the boundary
1322           INTEGER, INTENT(IN) :: BCV
1323           INTEGER, INTENT(IN) :: BCV_I
1324     
1325     ! Generic loop counter
1326           INTEGER :: LC1, OCCUPANTS
1327           INTEGER :: I, J, K
1328     
1329           OCCUPANTS = DEM_MI(BCV_I)%OCCUPANTS
1330           allocate(DEM_MI(BCV_I)%OWNER(OCCUPANTS))
1331     
1332           DEM_MI(BCV_I)%OWNER = 0
1333     
1334           SELECT CASE (BC_PLANE(BCV))
1335           CASE('N','S')
1336     
1337              J = DEM_MI(BCV_I)%L
1338              DO LC1=1,OCCUPANTS
1339                 I = DEM_MI(BCV_I)%W(LC1)
1340                 K = DEM_MI(BCV_I)%H(LC1)
1341                 IF(IS_ON_myPE_owns(I,J,K)) &
1342                    DEM_MI(BCV_I)%OWNER(LC1) = myPE
1343              ENDDO
1344     
1345           CASE('E','W')
1346     
1347              I = DEM_MI(BCV_I)%L
1348              DO LC1=1,OCCUPANTS
1349                 J = DEM_MI(BCV_I)%W(LC1)
1350                 K = DEM_MI(BCV_I)%H(LC1)
1351                 IF(IS_ON_myPE_owns(I,J,K)) &
1352                    DEM_MI(BCV_I)%OWNER(LC1) = myPE
1353              ENDDO
1354     
1355           CASE('T','B')
1356     
1357              K = DEM_MI(BCV_I)%L
1358              DO LC1=1,OCCUPANTS
1359                 I = DEM_MI(BCV_I)%W(LC1)
1360                 J = DEM_MI(BCV_I)%H(LC1)
1361                 IF(IS_ON_myPE_owns(I,J,K)) &
1362                    DEM_MI(BCV_I)%OWNER(LC1) = myPE
1363              ENDDO
1364     
1365           END SELECT
1366     
1367           CALL GLOBAL_ALL_SUM(DEM_MI(BCV_I)%OWNER(:))
1368     
1369           RETURN
1370           END SUBROUTINE SET_DEM_MI_OWNER
1371