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

1     !----------------------------------------------------------------------!
2     !  Module: MPI_INIT_DES                                                !
3     !  Author: Pradeep Gopalakrishnan                                      !
4     !                                                                      !
5     !  Purpose: Contains routines for setting up DES MPI communications.   !
6     !                                                                      !
7     !----------------------------------------------------------------------!
8           module mpi_init_des
9     
10     !-----------------------------------------------
11     ! Modules
12     !-----------------------------------------------
13           use parallel_mpi
14           use mpi_utility
15           use discretelement
16           use desgrid
17           use compar
18           use physprop
19           use sendrecv
20           use des_bc
21           use desmpi_wrapper
22           use sendrecvnode
23           use mfix_pic
24           use des_thermo
25           use run, only: ENERGY_EQ,ANY_SPECIES_EQ
26           use param, only: DIMENSION_N_s
27           use des_rxns
28     
29           use desmpi
30     
31     
32           contains
33     
34     
35     !----------------------------------------------------------------------!
36     !  Module: DESMPI_INIT                                                 !
37     !  Author: Pradeep Gopalakrishnan                                      !
38     !                                                                      !
39     !  Purpose: Allocates and initializes variables used by MPI send/recv  !
40     !  calls. Sets flags related to periodic boundaries and processor      !
41     !  interfaces.                                                         !
42     !                                                                      !
43     !----------------------------------------------------------------------!
44           subroutine desmpi_init
45     
46           use particle_filter, only: DES_INTERP_GARG
47     
48           use desmpi, only: iGhostPacketSize
49           use desmpi, only: iParticlePacketSize
50           use desmpi, only: iPairPacketSize
51     ! Particle orientation
52           use discretelement, only: PARTICLE_ORIENTATION
53           use discretelement, only: DES_USR_VAR_SIZE
54     
55           use error_manager
56     
57     !-----------------------------------------------
58           implicit none
59     !-----------------------------------------------
60     ! local variables
61     !-----------------------------------------------
62           integer :: lfaces
63           integer :: lmaxlen1,lmaxlen2,lmaxarea,lmaxghostpar,ii
64     
65           DOUBLE PRECISION, PARAMETER :: ONEMBo8 = 131072.0
66     
67     !-----------------------------------------------
68     
69     ! Calculate the size of ghost particle packets:
70           iGhostPacketSize = 15 + DES_USR_VAR_SIZE
71           IF(ENERGY_EQ) &
72              iGhostPacketSize = iGhostPacketSize + 1
73     
74     ! Calculate the size of particle packets.
75           iParticlePacketSize = 30 + DES_USR_VAR_SIZE
76           IF(ENERGY_EQ) &
77              iParticlePacketSize = iParticlePacketSize + 1 + DIMENSION_N_s
78           IF(DO_OLD) &
79              iParticlePacketSize = iParticlePacketSize + 15
80           IF(DO_OLD .AND. ENERGY_EQ) &
81              iParticlePacketSize = iParticlePacketSize + 1
82           IF(MPPIC) &
83              iParticlePacketSize = iParticlePacketSize + 1
84           IF(PARTICLE_ORIENTATION) &
85              iParticlePacketSize = iParticlePacketSize + 3
86           IF(DES_EXPLICITLY_COUPLED) THEN
87              iParticlePacketSize = iParticlePacketSize + 3
88              IF(ENERGY_EQ)iParticlePacketSize = iParticlePacketSize + 1
89              IF(ANY_SPECIES_EQ)iParticlePacketSize = iParticlePacketSize + 1
90           ENDIF
91     
92     ! Calculate the size of neighbor data
93           iPairPacketSize = 11
94     
95     ! Calculate the initial size of send and recv buffer based on max_pip,
96     ! total cells max number of boundary cells, and ghost par packet size.
97           lfaces = dimn*2
98           lmaxlen1 = dg_iend2-dg_istart2+1
99           lmaxlen2 = dg_jend2-dg_jstart2+1
100           if (do_K) then
101              lmaxlen1 = max(lmaxlen1,dg_kend2 -dg_kstart2+1)
102              lmaxlen2 = max(lmaxlen2,dg_kend2 -dg_kstart2+1)
103           else
104              lmaxlen1 = max(lmaxlen1,lmaxlen2)
105              lmaxlen2 = 1
106           end if
107     
108     
109     ! Note: 10 is added for buffer and is required for send and recv indices
110           lmaxarea = lmaxlen1*lmaxlen2 + 10
111     
112     ! Random value. This gets resized by DESMPI_CHECK_SENDRECVBUF
113           lmaxghostpar = 100
114           imaxbuf = lmaxghostpar*lmaxarea*iGhostPacketSize
115     
116           WRITE(ERR_MSG, 1000) iMAXBUF/ONEMBo8, ONEMBo8/iGhostPacketSize,  &
117              ONEMBo8/iParticlePacketSize, ONEMBo8/iPairPacketSize
118           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
119     
120      1000 FORMAT(/'DES MPI send/recv buffer: ',F7.1,' MB',/' o  ',F6.0,1X, &
121              'Ghost Particles/MB',/' o  ',F6.0,1X,'Particles/MB',/' o  ',  &
122              F6.0,1X,'Neighbor Pairs/MB')
123     
124     
125           allocate (dsendbuf(2));
126           allocate (drecvbuf(2));
127           do ii=1, size(dsendbuf)
128              allocate (dsendbuf(ii)%facebuf(imaxbuf));
129              allocate (drecvbuf(ii)%facebuf(imaxbuf));
130           end do
131     
132           allocate (isendindices(lmaxarea,lfaces)); isendindices=0
133           allocate (irecvindices(lmaxarea,lfaces)); irecvindices=0
134     
135           allocate (isendreq(lfaces)); isendreq=0
136           allocate (irecvreq(lfaces)); irecvreq=0
137           allocate (isendcnt(lfaces)); isendcnt=0
138     
139           allocate (dcycl_offset(lfaces,dimn)); dcycl_offset=0.0
140           allocate (ineighproc(lfaces)); ineighproc=0
141           allocate (iexchflag(lfaces)); iexchflag=.FALSE.
142     
143     ! allocate variables related to scattergather
144           allocate(iscattercnts(0:numpes-1)); iscattercnts=0
145           allocate(igathercnts(0:numpes-1));  igathercnts=0
146           allocate(idispls(0:numpes-1)); idispls=0
147     
148     ! set the communication flags
149           CALL DESMPI_SETCOMM
150     
151     !      call des_dbgmpi(1)
152     
153           end subroutine desmpi_init
154     
155     
156     !------------------------------------------------------------------------
157     ! Subroutine       : desmpi_setcomm
158     ! Purpose          : sets the flags required for interprocessor communication
159     !------------------------------------------------------------------------
160           subroutine desmpi_setcomm()
161     !-----------------------------------------------
162           implicit none
163     !-----------------------------------------------
164     ! local variables
165     !-----------------------------------------------
166           integer lijkproc,liproc,ljproc,lkproc
167           integer li,lj,lk,lis,lie,ljs,lje,lks,lke,lcount,lface
168           integer listart1,liend1,ljstart1,ljend1,lkstart1,lkend1
169           integer listart2,liend2,ljstart2,ljend2,lkstart2,lkend2
170     !-----------------------------------------------
171     
172     ! set flags for interprocessor boundaries and set the corresponding to proc
173           lijkproc = mype
174           liproc = iofproc(lijkproc)
175           ljproc = jofproc(lijkproc)
176           lkproc = kofproc(lijkproc)
177           iexchflag(:) = .false.
178           ineighproc(:) = 0
179           if(liproc.gt.0) then
180              iexchflag(2) = .true.
181              ineighproc(2)=procijk(liproc-1,ljproc,lkproc)
182           end if
183           if(liproc.lt.nodesi-1) then
184              iexchflag(1) = .true.
185              ineighproc(1)=procijk(liproc+1,ljproc,lkproc)
186           end if
187           if(ljproc.gt.0) then
188              iexchflag(4)= .true.
189              ineighproc(4)=procijk(liproc,ljproc-1,lkproc)
190           end if
191           if(ljproc.lt.nodesj-1) then
192              iexchflag(3) = .true.
193              ineighproc(3)=procijk(liproc,ljproc+1,lkproc)
194           end if
195           if(lkproc.gt.0) then
196              iexchflag(6)=.true.
197              ineighproc(6)=procijk(liproc,ljproc,lkproc-1)
198           end if
199           if(lkproc.lt.nodesk-1) then
200              iexchflag(5) =.true.
201              ineighproc(5)=procijk(liproc,ljproc,lkproc+1)
202           end if
203     
204     !set flags for cyclic boundary conditions and corresponding to proc
205           dcycl_offset(:,:) = 0
206           if (des_periodic_walls_x) then
207              if(liproc.eq.0) then
208                 iexchflag(2)=.true.
209                 ineighproc(2)= procijk(nodesi-1,ljproc,lkproc)
210                 dcycl_offset(2,1)= xlength
211              end if
212              if(liproc.eq.nodesi-1) then
213                 iexchflag(1)=.true.
214                 ineighproc(1)= procijk(0,ljproc,lkproc)
215                 dcycl_offset(1,1)=-xlength
216              end if
217           end if
218           if (des_periodic_walls_y) then
219              if(ljproc.eq.0) then
220                 iexchflag(4)=.true.
221                 ineighproc(4)= procijk(liproc,nodesj-1,lkproc)
222                 dcycl_offset(4,2)= ylength
223              end if
224              if(ljproc.eq.nodesj-1) then
225                 iexchflag(3)=.true.
226                 ineighproc(3)= procijk(liproc,0,lkproc)
227                 dcycl_offset(3,2)=-ylength
228              end if
229           end if
230           if (des_periodic_walls_z) then
231              if(lkproc.eq.0) then
232                 iexchflag(6)=.true.
233                 ineighproc(6)= procijk(liproc,ljproc,nodesk-1)
234                 dcycl_offset(6,3)=zlength
235              end if
236              if(lkproc.eq.nodesk-1) then
237                 iexchflag(5)=.true.
238                 ineighproc(5)= procijk(liproc,ljproc,0)
239                 dcycl_offset(5,3)=-zlength
240              end if
241           end if
242     
243           listart1=dg_istart1; liend1=dg_iend1
244           listart2=dg_istart2; liend2=dg_iend2
245     
246           ljstart1=dg_jstart1; ljend1=dg_jend1
247           ljstart2=dg_jstart2; ljend2=dg_jend2
248     
249           lkstart1=dg_kstart1; lkend1=dg_kend1
250           lkstart2=dg_kstart2; lkend2=dg_kend2
251     
252     ! Extend the domain indices to account for mass inlets and outlets. Do
253     ! not extend the domain for periodic walls becuase 1) they should not
254     ! include inflows or outflows and 2) they are already expanded.
255           IF(.NOT.DES_PERIODIC_WALLS_X) THEN
256              if(listart1.eq.dg_imin1) listart1 = dg_imin1-1
257              if(liend1.eq.dg_imax1) liend1 = dg_imax1+1
258           ENDIF
259     
260           IF(.NOT.DES_PERIODIC_WALLS_Y) THEN
261              if(ljstart1.eq.dg_jmin1) ljstart1 = dg_jmin1-1
262              if(ljend1.eq.dg_jmax1) ljend1 = dg_jmax1+1
263           ENDIF
264     
265           IF(DO_K .AND. .NOT.DES_PERIODIC_WALLS_Z) THEN
266              if(lkstart1.eq.dg_kmin1) lkstart1 = dg_kmin1-1
267              if(lkend1.eq.dg_kmax1) lkend1 = dg_kmax1+1
268           ENDIF
269     
270     ! set the ghost cell indices for e-w, n-s and t-b
271     ! for east and west faces
272           lks = lkstart1
273           lke = lkend1
274           ljs = ljstart1
275           lje = ljend1
276     
277     !east face
278           lface = 1
279           li = liend1
280           lcount = 1
281           do lk = lks,lke
282              do lj = ljs,lje
283                 lcount = lcount + 1
284                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
285                 irecvindices(lcount,lface) = dg_funijk(li+1,lj,lk)
286              end do
287           end do
288           isendindices(1,lface) = lcount - 1
289           irecvindices(1,lface) = lcount - 1
290     
291     !west face
292           lface = 2
293           li = listart1
294           lcount = 1
295           do lk = lks,lke
296              do lj = ljs,lje
297                 lcount = lcount + 1
298                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
299                 irecvindices(lcount,lface) = dg_funijk(li-1,lj,lk)
300              end do
301           end do
302           isendindices(1,lface) = lcount - 1
303           irecvindices(1,lface) = lcount - 1
304     
305     ! for north and south faces
306           lks = lkstart1
307           lke = lkend1
308           lis = listart2
309           lie = liend2
310     
311     !north face
312           lface = 3
313           lj = ljend1
314           lcount = 1
315           do lk = lks,lke
316              do li = lis,lie
317                 lcount = lcount + 1
318                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
319                 irecvindices(lcount,lface) = dg_funijk(li,lj+1,lk)
320              end do
321           end do
322           isendindices(1,lface) = lcount - 1
323           irecvindices(1,lface) = lcount - 1
324     
325     
326     !south face
327           lface = 4
328           lj = ljstart1
329           lcount = 1
330           do lk = lks,lke
331              do li = lis,lie
332                 lcount = lcount + 1
333                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
334                 irecvindices(lcount,lface) = dg_funijk(li,lj-1,lk)
335              end do
336           end do
337           isendindices(1,lface) = lcount - 1
338           irecvindices(1,lface) = lcount - 1
339     
340     ! for top and bottom
341           if (no_k) return
342           lis = listart2
343           lie = liend2
344           ljs = ljstart2
345           lje = ljend2
346     
347     !top face
348           lface = 5
349           lk = lkend1
350           lcount = 1
351           do li = lis,lie
352              do lj = ljs,lje
353                 lcount = lcount + 1
354                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
355                 irecvindices(lcount,lface) = dg_funijk(li,lj,lk+1)
356              end do
357           end do
358           isendindices(1,lface) = lcount - 1
359           irecvindices(1,lface) = lcount - 1
360     
361     
362     !bottom face
363           lface = 6
364           lk = lkstart1
365           lcount = 1
366           do li = lis,lie
367              do lj = ljs,lje
368                 lcount = lcount + 1
369                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
370                 irecvindices(lcount,lface) = dg_funijk(li,lj,lk-1)
371              end do
372           end do
373           isendindices(1,lface) = lcount - 1
374           irecvindices(1,lface) = lcount - 1
375     
376           return
377           end subroutine desmpi_setcomm
378     
379     
380     !------------------------------------------------------------------------
381     ! Subroutine       : des_scatter_particle
382     ! Purpose          : Scatters the particles read from input file or
383     !                    generated based on specified volume fraction
384     ! Comments         : In the main thread (pe_io) particles will be seperated
385     !                    based on the location and it will be packed in drootbuf
386     !                    Each proc receives its particles along with particle count
387     !------------------------------------------------------------------------
388           subroutine des_scatter_particle
389     
390           use mpi_comm_des, only: desmpi_scatterv
391           use des_allocate, only: particle_grow
392     
393     !-----------------------------------------------
394           implicit none
395     !-----------------------------------------------
396     ! local variables
397     !-----------------------------------------------
398           integer lcurpar,lproc,lbuf,lpacketsize
399           integer lproc_parcnt(0:numpes-1),lpar_proc(particles)
400     !-----------------------------------------------
401     
402           integer :: rdimn
403     
404           rdimn = merge(2,3, NO_K)
405     
406     ! set the packet size for transfer
407           lpacketsize = 2*rdimn + 2
408     
409     ! build the send buffer in PE_IO proc
410     ! first pass to get the count of particles
411           lpar_proc(:) =-1
412           lproc_parcnt(:) = 0
413           if(myPE.eq.pe_io) then
414              if (no_k) then
415                 do lcurpar = 1,particles
416                    do lproc= 0,numpes-1
417                       if (   dpar_pos(lcurpar,1).ge.xe(istart1_all(lproc)-1) &
418                        .and. dpar_pos(lcurpar,1).lt.xe(iend1_all(lproc))     &
419                        .and. dpar_pos(lcurpar,2).ge.yn(jstart1_all(lproc)-1) &
420                        .and. dpar_pos(lcurpar,2).lt.yn(jend1_all(lproc))) then
421                          lpar_proc(lcurpar) = lproc
422                          lproc_parcnt(lproc) = lproc_parcnt(lproc) + 1
423                          exit
424                       endif
425                    enddo
426                    if (lpar_proc(lcurpar).eq.-1) then
427                       WRITE(*,500) lcurpar
428                       call des_mpi_stop
429                    endif
430                 enddo
431              else
432                 do lcurpar = 1,particles
433                    do lproc= 0,numpes-1
434                       if (   dpar_pos(lcurpar,1).ge.xe(istart1_all(lproc)-1) &
435                        .and. dpar_pos(lcurpar,1).lt.xe(iend1_all(lproc))     &
436                        .and. dpar_pos(lcurpar,2).ge.yn(jstart1_all(lproc)-1) &
437                        .and. dpar_pos(lcurpar,2).lt.yn(jend1_all(lproc))     &
438                        .and. dpar_pos(lcurpar,3).ge.zt(kstart1_all(lproc)-1) &
439                        .and. dpar_pos(lcurpar,3).lt.zt(kend1_all(lproc))) then
440                          lpar_proc(lcurpar) = lproc
441                          lproc_parcnt(lproc) = lproc_parcnt(lproc) + 1
442                          exit
443                       end if
444                    end do
445                    if (lpar_proc(lcurpar).eq.-1) then
446                       WRITE(*,501) lcurpar
447                       call des_mpi_stop
448                    endif
449                 enddo
450              endif  ! if (no_k)
451           endif ! if (my_pe.eq.pe_io)
452           call bcast(lproc_parcnt(0:numpes-1),pe_io)
453     
454     ! second pass: set and allocate scatter related variables
455           pip = lproc_parcnt(mype)
456           call PARTICLE_GROW(pip)
457           max_pip = max(pip,max_pip)
458           iscr_recvcnt = pip*lpacketsize
459           allocate (dprocbuf(iscr_recvcnt))
460           if (mype.eq.pe_io) then
461              allocate (drootbuf(particles*lpacketsize))
462           else
463              allocate (drootbuf(10))
464           endif
465     
466     ! in the IO processor build the drootbuffer and idispls required
467     ! for mpi communication
468           if(mype.eq.pe_io) then
469              idispls(0) = 0
470              iscattercnts(0) = lproc_parcnt(0)*lpacketsize
471              do lproc = 1,numpes-1
472                 idispls(lproc) = idispls(lproc-1) + iscattercnts(lproc-1)
473                 iscattercnts(lproc) = lproc_parcnt(lproc)*lpacketsize
474              end do
475              lproc_parcnt(:) = 0
476              do lcurpar = 1,particles
477                 lproc = lpar_proc(lcurpar)
478                 lbuf = idispls(lproc)+lproc_parcnt(lproc)*lpacketsize+1
479                 drootbuf(lbuf) = dpar_rad(lcurpar); lbuf = lbuf + 1
480                 drootbuf(lbuf) = dpar_den(lcurpar); lbuf = lbuf + 1
481                 drootbuf(lbuf:lbuf+rdimn-1) = dpar_pos(lcurpar,1:rdimn); lbuf = lbuf + rdimn
482                 drootbuf(lbuf:lbuf+rdimn-1) = dpar_vel(lcurpar,1:rdimn); lbuf = lbuf + rdimn
483                 lproc_parcnt(lproc) = lproc_parcnt(lproc) + 1
484              enddo
485           endif
486           call desmpi_scatterv(ptype=2)
487     
488     ! unpack the particles in each processor and set the pip
489           do lcurpar = 1,pip
490              lbuf = (lcurpar-1)*lpacketsize+1
491              des_radius(lcurpar) = dprocbuf(lbuf); lbuf = lbuf+1
492              ro_sol(lcurpar) = dprocbuf(lbuf); lbuf = lbuf+1
493              des_pos_new(lcurpar,1:rdimn) = dprocbuf(lbuf:lbuf+rdimn-1); lbuf = lbuf+rdimn
494              des_vel_new(lcurpar,1:rdimn) = dprocbuf(lbuf:lbuf+rdimn-1); lbuf = lbuf+rdimn
495              call set_normal(lcurpar)
496           enddo
497           deallocate (dprocbuf,drootbuf)
498     
499      500  FORMAT(/2X,'From: DES_SCATTER_PARTICLE: (0)',/2X,&
500              'ERROR: Unable to locate the particle (no. ',I10,&
501              ') inside the domain')
502      501  FORMAT(/2X,'From: DES_SCATTER_PARTICLE: (1)',/2X,&
503              'ERROR: Unable to locate the particle (no. ',I10,&
504              ') inside the domain')
505     
506           RETURN
507           end subroutine des_scatter_particle
508     
509     
510     
511     !------------------------------------------------------------------------
512     ! Subroutine       : DES_RESTART_GHOST
513     ! Purpose          : restart file contains neighbour information in terms
514     !                    global id. This routine converts the global id into
515     !                    local particle number
516     !                    steps
517     !                    1. Exchange the ghost particles (does not involve any neighbour info)
518     !                    2. loop through particles neighbour and contact list
519     !                       and convert global numbers to local numbers
520     ! Parameters       : none
521     !------------------------------------------------------------------------
522     
523           subroutine DES_RESTART_GHOST
524     
525           use mpi_comm_des, only: desmpi_sendrecv_init
526           use mpi_comm_des, only: desmpi_sendrecv_wait
527     
528           use mpi_funs_des, only: desmpi_check_sendrecvbuf
529           use mpi_pack_des, only: desmpi_pack_ghostpar
530           use mpi_unpack_des, only: desmpi_unpack_ghostpar
531     
532     !-----------------------------------------------
533           implicit none
534     !-----------------------------------------------
535     ! local variables
536     !-----------------------------------------------
537           integer linter,lface
538           integer ii
539     !-----------------------------------------------
540     ! set do_nsearch true so that the ghost cell will be updated
541           do_nsearch = .true.
542           call desgrid_pic(plocate=.true.)
543           call desmpi_check_sendrecvbuf(check_global=.true.)
544     
545     !call ghost particle exchange in E-W, N-S, T-B order
546     
547           do ii=1, size(dsendbuf)
548              dsendbuf(ii)%facebuf(1) = 0
549              drecvbuf(ii)%facebuf(1) = 0
550           end do
551     
552           ighost_updated(:) = .false.
553           ispot = 1
554           do linter = 1,dimn
555              do lface = linter*2-1,linter*2
556                 if(.not.iexchflag(lface))cycle
557                 call desmpi_pack_ghostpar(lface)
558                 call desmpi_sendrecv_init(lface)
559              enddo
560              do lface = linter*2-1,linter*2
561                 if(.not.iexchflag(lface)) cycle
562                 call desmpi_sendrecv_wait(lface)
563                 call desmpi_unpack_ghostpar(lface)
564              enddo
565     ! update pic required as particles in ghost cell can move between ghost cells
566              do lface = linter*2-1,linter*2
567                 if(dsendbuf(1+mod(lface,2))%facebuf(1).gt.0.or.drecvbuf(1+mod(lface,2))%facebuf(1).gt.0) then
568                    call desgrid_pic(plocate=.false.)
569                    exit
570                 endif
571              enddo
572           enddo
573           call des_mpi_barrier
574     
575           end subroutine DES_RESTART_GHOST
576     
577           end module mpi_init_des
578