File: RELATIVE:/../../../mfix.git/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,lfactor=4
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
78           IF(ANY_SPECIES_EQ) &
79              iParticlePacketSize = iParticlePacketSize + DIMENSION_N_s
80           IF(DES_EXPLICITLY_COUPLED) &
81              iParticlePacketSize = iParticlePacketSize + 3
82           IF(DO_OLD) &
83              iParticlePacketSize = iParticlePacketSize + 15
84           IF(DO_OLD .AND. ENERGY_EQ) &
85              iParticlePacketSize = iParticlePacketSize + 1
86           IF(MPPIC) &
87              iParticlePacketSize = iParticlePacketSize + 1
88           IF(PARTICLE_ORIENTATION) &
89              iParticlePacketSize = iParticlePacketSize + 3
90     
91     ! Calculate the size of neighbor data
92           iPairPacketSize = 11
93     
94     ! Calculate the initial size of send and recv buffer based on max_pip,
95     ! total cells max number of boundary cells, and ghost par packet size.
96           lfaces = dimn*2
97           lmaxlen1 = dg_iend2-dg_istart2+1
98           lmaxlen2 = dg_jend2-dg_jstart2+1
99           if (do_K) then
100              lmaxlen1 = max(lmaxlen1,dg_kend2 -dg_kstart2+1)
101              lmaxlen2 = max(lmaxlen2,dg_kend2 -dg_kstart2+1)
102           else
103              lmaxlen1 = max(lmaxlen1,lmaxlen2)
104              lmaxlen2 = 1
105           end if
106     
107     
108     ! Note: 10 is added for buffer and is required for send and recv indices
109           lmaxarea = lmaxlen1*lmaxlen2 + 10
110     
111     ! Random value. This gets resized by DESMPI_CHECK_SENDRECVBUF
112           lmaxghostpar = 100
113           imaxbuf = lmaxghostpar*lmaxarea*iGhostPacketSize
114     
115           WRITE(ERR_MSG, 1000) iMAXBUF/ONEMBo8, ONEMBo8/iGhostPacketSize,  &
116              ONEMBo8/iParticlePacketSize, ONEMBo8/iPairPacketSize
117           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
118     
119      1000 FORMAT(/'DES MPI send/recv buffer: ',F7.1,' MB',/' o  ',F6.0,1X, &
120              'Ghost Particles/MB',/' o  ',F6.0,1X,'Particles/MB',/' o  ',  &
121              F6.0,1X,'Neighbor Pairs/MB')
122     
123     
124           allocate (dsendbuf(2));
125           allocate (drecvbuf(2));
126           do ii=1, size(dsendbuf)
127              allocate (dsendbuf(ii)%facebuf(imaxbuf));
128              allocate (drecvbuf(ii)%facebuf(imaxbuf));
129           end do
130     
131           allocate (isendindices(lmaxarea,lfaces)); isendindices=0
132           allocate (irecvindices(lmaxarea,lfaces)); irecvindices=0
133     
134           allocate (isendreq(lfaces)); isendreq=0
135           allocate (irecvreq(lfaces)); irecvreq=0
136           allocate (isendcnt(lfaces)); isendcnt=0
137     
138           allocate (dcycl_offset(lfaces,dimn)); dcycl_offset=0.0
139           allocate (ineighproc(lfaces)); ineighproc=0
140           allocate (iexchflag(lfaces)); iexchflag=.FALSE.
141     
142     ! allocate variables related to scattergather
143           allocate(iscattercnts(0:numpes-1)); iscattercnts=0
144           allocate(igathercnts(0:numpes-1));  igathercnts=0
145           allocate(idispls(0:numpes-1)); idispls=0
146     
147     ! set the communication flags
148           CALL DESMPI_SETCOMM
149     
150     !      call des_dbgmpi(1)
151     
152           end subroutine desmpi_init
153     
154     
155     !------------------------------------------------------------------------
156     ! Subroutine       : desmpi_setcomm
157     ! Purpose          : sets the flags required for interprocessor communication
158     !------------------------------------------------------------------------
159           subroutine desmpi_setcomm()
160     !-----------------------------------------------
161           implicit none
162     !-----------------------------------------------
163     ! local variables
164     !-----------------------------------------------
165           integer lijkproc,liproc,ljproc,lkproc
166           integer li,lj,lk,lis,lie,ljs,lje,lks,lke,lcount,lface
167           integer listart1,liend1,ljstart1,ljend1,lkstart1,lkend1
168           integer listart2,liend2,ljstart2,ljend2,lkstart2,lkend2
169     !-----------------------------------------------
170     
171     ! set flags for interprocessor boundaries and set the corresponding to proc
172           lijkproc = mype
173           liproc = iofproc(lijkproc)
174           ljproc = jofproc(lijkproc)
175           lkproc = kofproc(lijkproc)
176           iexchflag(:) = .false.
177           ineighproc(:) = 0
178           if(liproc.gt.0) then
179              iexchflag(2) = .true.
180              ineighproc(2)=procijk(liproc-1,ljproc,lkproc)
181           end if
182           if(liproc.lt.nodesi-1) then
183              iexchflag(1) = .true.
184              ineighproc(1)=procijk(liproc+1,ljproc,lkproc)
185           end if
186           if(ljproc.gt.0) then
187              iexchflag(4)= .true.
188              ineighproc(4)=procijk(liproc,ljproc-1,lkproc)
189           end if
190           if(ljproc.lt.nodesj-1) then
191              iexchflag(3) = .true.
192              ineighproc(3)=procijk(liproc,ljproc+1,lkproc)
193           end if
194           if(lkproc.gt.0) then
195              iexchflag(6)=.true.
196              ineighproc(6)=procijk(liproc,ljproc,lkproc-1)
197           end if
198           if(lkproc.lt.nodesk-1) then
199              iexchflag(5) =.true.
200              ineighproc(5)=procijk(liproc,ljproc,lkproc+1)
201           end if
202     
203     !set flags for cyclic boundary conditions and corresponding to proc
204           dcycl_offset(:,:) = 0
205           if (des_periodic_walls_x) then
206              if(liproc.eq.0) then
207                 iexchflag(2)=.true.
208                 ineighproc(2)= procijk(nodesi-1,ljproc,lkproc)
209                 dcycl_offset(2,1)= xlength
210              end if
211              if(liproc.eq.nodesi-1) then
212                 iexchflag(1)=.true.
213                 ineighproc(1)= procijk(0,ljproc,lkproc)
214                 dcycl_offset(1,1)=-xlength
215              end if
216           end if
217           if (des_periodic_walls_y) then
218              if(ljproc.eq.0) then
219                 iexchflag(4)=.true.
220                 ineighproc(4)= procijk(liproc,nodesj-1,lkproc)
221                 dcycl_offset(4,2)= ylength
222              end if
223              if(ljproc.eq.nodesj-1) then
224                 iexchflag(3)=.true.
225                 ineighproc(3)= procijk(liproc,0,lkproc)
226                 dcycl_offset(3,2)=-ylength
227              end if
228           end if
229           if (des_periodic_walls_z) then
230              if(lkproc.eq.0) then
231                 iexchflag(6)=.true.
232                 ineighproc(6)= procijk(liproc,ljproc,nodesk-1)
233                 dcycl_offset(6,3)=zlength
234              end if
235              if(lkproc.eq.nodesk-1) then
236                 iexchflag(5)=.true.
237                 ineighproc(5)= procijk(liproc,ljproc,0)
238                 dcycl_offset(5,3)=-zlength
239              end if
240           end if
241     
242           listart1=dg_istart1; liend1=dg_iend1
243           listart2=dg_istart2; liend2=dg_iend2
244     
245           ljstart1=dg_jstart1; ljend1=dg_jend1
246           ljstart2=dg_jstart2; ljend2=dg_jend2
247     
248           lkstart1=dg_kstart1; lkend1=dg_kend1
249           lkstart2=dg_kstart2; lkend2=dg_kend2
250     
251     ! Extend the domain indices to account for mass inlets and outlets. Do
252     ! not extend the domain for periodic walls becuase 1) they should not
253     ! include inflows or outflows and 2) they are already expanded.
254           IF(.NOT.DES_PERIODIC_WALLS_X) THEN
255              if(listart1.eq.dg_imin1) listart1 = dg_imin1-1
256              if(liend1.eq.dg_imax1) liend1 = dg_imax1+1
257           ENDIF
258     
259           IF(.NOT.DES_PERIODIC_WALLS_Y) THEN
260              if(ljstart1.eq.dg_jmin1) ljstart1 = dg_jmin1-1
261              if(ljend1.eq.dg_jmax1) ljend1 = dg_jmax1+1
262           ENDIF
263     
264           IF(DO_K .AND. .NOT.DES_PERIODIC_WALLS_Z) THEN
265              if(lkstart1.eq.dg_kmin1) lkstart1 = dg_kmin1-1
266              if(lkend1.eq.dg_kmax1) lkend1 = dg_kmax1+1
267           ENDIF
268     
269     ! set the ghost cell indices for e-w, n-s and t-b
270     ! for east and west faces
271           lks = lkstart1
272           lke = lkend1
273           ljs = ljstart1
274           lje = ljend1
275     
276     !east face
277           lface = 1
278           li = liend1
279           lcount = 1
280           do lk = lks,lke
281              do lj = ljs,lje
282                 lcount = lcount + 1
283                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
284                 irecvindices(lcount,lface) = dg_funijk(li+1,lj,lk)
285              end do
286           end do
287           isendindices(1,lface) = lcount - 1
288           irecvindices(1,lface) = lcount - 1
289     
290     !west face
291           lface = 2
292           li = listart1
293           lcount = 1
294           do lk = lks,lke
295              do lj = ljs,lje
296                 lcount = lcount + 1
297                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
298                 irecvindices(lcount,lface) = dg_funijk(li-1,lj,lk)
299              end do
300           end do
301           isendindices(1,lface) = lcount - 1
302           irecvindices(1,lface) = lcount - 1
303     
304     ! for north and south faces
305           lks = lkstart1
306           lke = lkend1
307           lis = listart2
308           lie = liend2
309     
310     !north face
311           lface = 3
312           lj = ljend1
313           lcount = 1
314           do lk = lks,lke
315              do li = lis,lie
316                 lcount = lcount + 1
317                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
318                 irecvindices(lcount,lface) = dg_funijk(li,lj+1,lk)
319              end do
320           end do
321           isendindices(1,lface) = lcount - 1
322           irecvindices(1,lface) = lcount - 1
323     
324     
325     !south face
326           lface = 4
327           lj = ljstart1
328           lcount = 1
329           do lk = lks,lke
330              do li = lis,lie
331                 lcount = lcount + 1
332                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
333                 irecvindices(lcount,lface) = dg_funijk(li,lj-1,lk)
334              end do
335           end do
336           isendindices(1,lface) = lcount - 1
337           irecvindices(1,lface) = lcount - 1
338     
339     ! for top and bottom
340           if (no_k) return
341           lis = listart2
342           lie = liend2
343           ljs = ljstart2
344           lje = ljend2
345     
346     !top face
347           lface = 5
348           lk = lkend1
349           lcount = 1
350           do li = lis,lie
351              do lj = ljs,lje
352                 lcount = lcount + 1
353                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
354                 irecvindices(lcount,lface) = dg_funijk(li,lj,lk+1)
355              end do
356           end do
357           isendindices(1,lface) = lcount - 1
358           irecvindices(1,lface) = lcount - 1
359     
360     
361     !bottom face
362           lface = 6
363           lk = lkstart1
364           lcount = 1
365           do li = lis,lie
366              do lj = ljs,lje
367                 lcount = lcount + 1
368                 isendindices(lcount,lface) = dg_funijk(li,lj,lk)
369                 irecvindices(lcount,lface) = dg_funijk(li,lj,lk-1)
370              end do
371           end do
372           isendindices(1,lface) = lcount - 1
373           irecvindices(1,lface) = lcount - 1
374     
375           return
376           end subroutine desmpi_setcomm
377     
378     
379     !------------------------------------------------------------------------
380     ! Subroutine       : des_scatter_particle
381     ! Purpose          : Scatters the particles read from input file or
382     !                    generated based on specified volume fraction
383     ! Comments         : In the main thread (pe_io) particles will be seperated
384     !                    based on the location and it will be packed in drootbuf
385     !                    Each proc receives its particles along with particle count
386     !------------------------------------------------------------------------
387           subroutine des_scatter_particle
388     
389           use mpi_comm_des, only: desmpi_scatterv
390           use des_allocate, only: particle_grow
391     
392     !-----------------------------------------------
393           implicit none
394     !-----------------------------------------------
395     ! local variables
396     !-----------------------------------------------
397           integer lcurpar,lproc,lbuf,lpacketsize
398           integer lproc_parcnt(0:numpes-1),lpar_proc(particles)
399     !-----------------------------------------------
400     
401           integer :: rdimn
402     
403           rdimn = merge(2,3, NO_K)
404     
405     ! set the packet size for transfer
406           lpacketsize = 2*rdimn + 2
407     
408     ! build the send buffer in PE_IO proc
409     ! first pass to get the count of particles
410           lpar_proc(:) =-1
411           lproc_parcnt(:) = 0
412           if(myPE.eq.pe_io) then
413              if (no_k) then
414                 do lcurpar = 1,particles
415                    do lproc= 0,numpes-1
416                       if (   dpar_pos(lcurpar,1).ge.xe(istart1_all(lproc)-1) &
417                        .and. dpar_pos(lcurpar,1).lt.xe(iend1_all(lproc))     &
418                        .and. dpar_pos(lcurpar,2).ge.yn(jstart1_all(lproc)-1) &
419                        .and. dpar_pos(lcurpar,2).lt.yn(jend1_all(lproc))) then
420                          lpar_proc(lcurpar) = lproc
421                          lproc_parcnt(lproc) = lproc_parcnt(lproc) + 1
422                          exit
423                       endif
424                    enddo
425                    if (lpar_proc(lcurpar).eq.-1) then
426                       WRITE(*,500) lcurpar
427                       call des_mpi_stop
428                    endif
429                 enddo
430              else
431                 do lcurpar = 1,particles
432                    do lproc= 0,numpes-1
433                       if (   dpar_pos(lcurpar,1).ge.xe(istart1_all(lproc)-1) &
434                        .and. dpar_pos(lcurpar,1).lt.xe(iend1_all(lproc))     &
435                        .and. dpar_pos(lcurpar,2).ge.yn(jstart1_all(lproc)-1) &
436                        .and. dpar_pos(lcurpar,2).lt.yn(jend1_all(lproc))     &
437                        .and. dpar_pos(lcurpar,3).ge.zt(kstart1_all(lproc)-1) &
438                        .and. dpar_pos(lcurpar,3).lt.zt(kend1_all(lproc))) then
439                          lpar_proc(lcurpar) = lproc
440                          lproc_parcnt(lproc) = lproc_parcnt(lproc) + 1
441                          exit
442                       end if
443                    end do
444                    if (lpar_proc(lcurpar).eq.-1) then
445                       WRITE(*,501) lcurpar
446                       call des_mpi_stop
447                    endif
448                 enddo
449              endif  ! if (no_k)
450           endif ! if (my_pe.eq.pe_io)
451           call bcast(lproc_parcnt(0:numpes-1),pe_io)
452     
453     ! second pass: set and allocate scatter related variables
454           pip = lproc_parcnt(mype)
455           if (pip .gt. max_pip) then
456              call PARTICLE_GROW(pip)
457           endif
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(1:rdimn,lcurpar) = dprocbuf(lbuf:lbuf+rdimn-1); lbuf = lbuf+rdimn
494              des_vel_new(1:rdimn,lcurpar) = 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
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      800  FORMAT(/2X,'From: DES_RESTART_GHOST: ',/2X,&
576              'WARNING: Unable to locate neighbor during restart (0)',/)
577      801  FORMAT(/2X,'From: DES_RESTART_GHOST: ',/2X,&
578              'WARNING: Unable to locate neighbor during restart (1)',/)
579     
580           end subroutine DES_RESTART_GHOST
581     
582           end module mpi_init_des
583