MFIX  2016-1
layout_mi_dem.f
Go to the documentation of this file.
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
stl_functions_des::tri_box_overlap
subroutine tri_box_overlap(pCENTER, pHALFSIZE, pVERTS, pOVERLAP)
Definition: stl_functions_des_mod.f:161
mpi_utility
Definition: mpi_utility_mod.f:6
bc::bc_y_n
double precision, dimension(dimension_bc) bc_y_n
Definition: bc_mod.f:42
des_bc::exclude_dem_mi_cell
logical function exclude_dem_mi_cell(lI, lJ, lK)
Definition: des_bc_mod.f:127
funits::dmp_log
logical dmp_log
Definition: funits_mod.f:6
layout_dem_mi_ew
subroutine layout_dem_mi_ew(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
Definition: layout_mi_dem.f:487
error_manager::finl_err_msg
subroutine finl_err_msg
Definition: error_manager_mod.f:216
layout_dem_mi_tb
subroutine layout_dem_mi_tb(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
Definition: layout_mi_dem.f:905
param1::one
double precision, parameter one
Definition: param1_mod.f:29
stl::vertex
double precision, dimension(3, 3, dim_stl) vertex
Definition: stl_mod.f:20
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
stl::stl_start
integer, dimension(4) stl_start
Definition: stl_mod.f:74
bc::bc_x_e
double precision, dimension(dimension_bc) bc_x_e
Definition: bc_mod.f:34
stl_functions_des
Definition: stl_functions_des_mod.f:10
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
geometry::imax
integer imax
Definition: geometry_mod.f:47
bc::bc_y_s
double precision, dimension(dimension_bc) bc_y_s
Definition: bc_mod.f:38
bc::bc_plane
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
error_manager::init_err_msg
subroutine init_err_msg(CALLER)
Definition: error_manager_mod.f:171
stl::default_stl
integer, parameter default_stl
Definition: stl_mod.f:73
cutcell::use_stl
logical use_stl
Definition: cutcell_mod.f:428
mpi_utility::global_all_max
Definition: mpi_utility_mod.f:110
stl
Definition: stl_mod.f:1
des_bc
Definition: des_bc_mod.f:14
run::run_type
character(len=16) run_type
Definition: run_mod.f:33
stl::norm_face
double precision, dimension(3, dim_stl) norm_face
Definition: stl_mod.f:22
geometry::dx
double precision, dimension(0:dim_i) dx
Definition: geometry_mod.f:68
calc_cell_intersect
subroutine calc_cell_intersect(RMIN, LOC, D_DIR, N_DIR, CELL)
Definition: calc_cell.f:127
mpi_utility::global_all_sum
Definition: mpi_utility_mod.f:66
layout_dem_mi_ns
subroutine layout_dem_mi_ns(BCV, BCV_I, MAX_DIA, setDBG, showMAP)
Definition: layout_mi_dem.f:72
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
mpi_utility::global_all_min
Definition: mpi_utility_mod.f:88
geometry::kmax
integer kmax
Definition: geometry_mod.f:51
set_dem_mi_owner
subroutine set_dem_mi_owner(BCV, BCV_I)
Definition: layout_mi_dem.f:1303
stl::n_facets_des
integer n_facets_des
Definition: stl_mod.f:18
bc::bc_z_b
double precision, dimension(dimension_bc) bc_z_b
Definition: bc_mod.f:46
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
compar::mype
integer mype
Definition: compar_mod.f:24
error_manager::err_msg
character(len=line_length), dimension(line_count) err_msg
Definition: error_manager_mod.f:29
geometry::xmin
double precision xmin
Definition: geometry_mod.f:75
geometry::jmax
integer jmax
Definition: geometry_mod.f:49
bc::bc_z_t
double precision, dimension(dimension_bc) bc_z_t
Definition: bc_mod.f:50
error_manager
Definition: error_manager_mod.f:8
des_bc::dem_mi
type(dem_mi_), dimension(:), allocatable, target dem_mi
Definition: des_bc_mod.f:90
layout_mi_dem
subroutine layout_mi_dem(BCV, BCV_I, MAX_DIA)
Definition: layout_mi_dem.f:13
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
error_manager::flush_err_msg
subroutine flush_err_msg(DEBUG, HEADER, FOOTER, ABORT, LOG, CALL_TREE)
Definition: error_manager_mod.f:305
bc
Definition: bc_mod.f:23
bc::bc_area
double precision, dimension(dimension_bc) bc_area
Definition: bc_mod.f:245
bc::bc_x_w
double precision, dimension(dimension_bc) bc_x_w
Definition: bc_mod.f:30