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