KTH framework for Nek5000 toolboxes; testing version  0.0.1
lpm_comm.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine lpm_comm_setup
3  include 'SIZE'
4  include 'TOTAL'
5  include 'CTIMER'
6 # include "LPM"
7 
8  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
9 
10  nmsh = lpm_rparam(3)
11 
12  call interp_setup(i_fp_hndl,0.0,nmsh,nelt)
13  call fgslib_crystal_setup(i_cr_hndl,nekcomm,np)
14 
15  return
16  end
17 c-----------------------------------------------------------------------
18  subroutine lpm_comm_findpts
19  include 'SIZE'
20  include 'TOTAL'
21 # include "LPM"
22 
23  common /intp_h/ ih_intp(2,1)
24 
25  ih_intp2 = ih_intp(2,i_fp_hndl)
26 
27  ix = 1
28  iy = 2
29  iz = 3
30 
31  call fgslib_findpts(ih_intp2 ! call fgslib_findpts( ihndl,
32  $ , lpm_iprop(1 ,1),lpm_lip ! $ rcode,1,
33  $ , lpm_iprop(3 ,1),lpm_lip ! & proc,1,
34  $ , lpm_iprop(2 ,1),lpm_lip ! & elid,1,
35  $ , lpm_rprop2(1 ,1),lpm_lrp2 ! & rst,ndim,
36  $ , lpm_rprop2(4 ,1),lpm_lrp2 ! & dist,1,
37  $ , lpm_y(ix,1),lpm_lrs ! & pts( 1),1,
38  $ , lpm_y(iy,1),lpm_lrs ! & pts( n+1),1,
39  $ , lpm_y(iz,1),lpm_lrs ,lpm_npart) ! & pts(2*n+1),1,n)
40 
41  do i=1,lpm_npart
42  lpm_iprop(4,i) = lpm_iprop(3,i)
43  enddo
44 
45  ! instead get which bin it is in
46  do i=1,lpm_npart
47  ! check if particles are greater or less than binb bounds....
48  ! add below....
49 
50  ii = floor((lpm_y(ix,i)-lpm_binb(1))/lpm_rdxgp)
51  jj = floor((lpm_y(iy,i)-lpm_binb(3))/lpm_rdygp)
52  kk = floor((lpm_y(iz,i)-lpm_binb(5))/lpm_rdzgp)
53  if (.not. if3d) kk = 0
54  if (ii .eq. lpm_ndxgp) ii = lpm_ndxgp - 1
55  if (jj .eq. lpm_ndygp) jj = lpm_ndygp - 1
56  if (kk .eq. lpm_ndzgp) kk = lpm_ndzgp - 1
57  ndum = ii + lpm_ndxgp*jj + lpm_ndxgp*lpm_ndygp*kk
58  nrank = modulo(ndum, np)
59 
60  lpm_iprop(8,i) = ii
61  lpm_iprop(9,i) = jj
62  lpm_iprop(10,i) = kk
63  lpm_iprop(11,i) = ndum
64 
65  lpm_iprop(4,i) = nrank ! where particle is actually moved
66  enddo
67 
68  n = nid_glcount(lpm_iprop(4,1),lpm_lip,lpm_npart)
69  ierr = 0
70  if (n.gt.lpm_lpart) ierr = 1
71  ierr = iglsum(ierr,1)
72  if (ierr.gt.0) then
73  nmax = iglmax(n,1)
74  call exitti('LPM_LPART too small, require >$',nmax)
75  endif
76 
77  return
78  end
79 c-----------------------------------------------------------------------
80  integer function nid_glcount(a,is,n)
81 c
82 c returns global nid count in distributed integer array
83 c
84  include 'mpif.h'
85  include 'SIZE'
86  include 'PARALLEL'
87 
88  integer a(*)
89 
90  common /nekmpi/ nid_,np_,nekcomm,nekgroup,nekreal
91 
92  integer disp_unit
93  integer*8 wsize, tdisp
94  data tdisp /0/
95 
96  integer win, shared_counter
97  save win, shared_counter
98 
99  integer one
100  parameter(one = 1)
101 
102  integer icalld
103  data icalld /0/
104  save icalld
105 
106 #ifdef MPI
107  if (icalld.eq.0) then
108  disp_unit = isize
109  wsize = disp_unit
110  call mpi_win_create(shared_counter,
111  $ wsize,
112  $ disp_unit,
113  $ mpi_info_null,
114  $ nekcomm,win,ierr)
115  icalld = 1
116  endif
117 
118  shared_counter = 0
119 
120  call mpi_win_fence(mpi_mode_noprecede,win,ierr)
121  do i = 1,n
122  call mpi_accumulate(one,1,mpi_integer,a((i-1)*is+1),
123  $ tdisp,1,mpi_integer,mpi_sum,win,ierr)
124  enddo
125  call mpi_win_fence(mpi_mode_nosucceed,win,ierr)
126 
127  nid_glcount = shared_counter
128 #else
129  nid_glcount = n
130 #endif
131 
132  return
133  end
134 c-----------------------------------------------------------------------
135  subroutine lpm_comm_crystal
136  include 'SIZE'
137  include 'TOTAL'
138  include 'CTIMER'
139 # include "LPM"
140 
141  logical partl
142  integer lpm_ipmap1(1,LPM_LPART)
143  > ,lpm_ipmap2(1,LPM_LPART)
144  > ,lpm_ipmap3(1,LPM_LPART)
145  > ,lpm_ipmap4(1,LPM_LPART)
146 
147  parameter(lrf = lpm_lrs*4 + lpm_lrp + lpm_lrp2)
148  real rwork(lrf,LPM_LPART)
149 
150  do i=1,lpm_npart
151  ic = 1
152  call copy(rwork(ic,i),lpm_y(1,i),lpm_lrs)
153  ic = ic + lpm_lrs
154  call copy(rwork(ic,i),lpm_y1((i-1)*lpm_lrs+1),lpm_lrs)
155  ic = ic + lpm_lrs
156  call copy(rwork(ic,i),lpm_ydot(1,i),lpm_lrs)
157  ic = ic + lpm_lrs
158  call copy(rwork(ic,i),lpm_ydotc(1,i),lpm_lrs)
159  ic = ic + lpm_lrs
160  call copy(rwork(ic,i),lpm_rprop(1,i),lpm_lrp)
161  ic = ic + lpm_lrp
162  call copy(rwork(ic,i),lpm_rprop2(1,i),lpm_lrp2)
163  enddo
164 
165  j0 = 4 ! proc key
166  call fgslib_crystal_tuple_transfer(i_cr_hndl,lpm_npart ,lpm_lpart
167  $ ,lpm_iprop ,lpm_lip,partl,0,rwork,lrf ,j0)
168 
169  do i=1,lpm_npart
170  ic = 1
171  call copy(lpm_y(1,i),rwork(ic,i),lpm_lrs)
172  ic = ic + lpm_lrs
173  call copy(lpm_y1((i-1)*lpm_lrs+1),rwork(ic,i),lpm_lrs)
174  ic = ic + lpm_lrs
175  call copy(lpm_ydot(1,i),rwork(ic,i),lpm_lrs)
176  ic = ic + lpm_lrs
177  call copy(lpm_ydotc(1,i),rwork(ic,i),lpm_lrs)
178  ic = ic + lpm_lrs
179  call copy(lpm_rprop(1,i),rwork(ic,i),lpm_lrp)
180  ic = ic + lpm_lrp
181  call copy(lpm_rprop2(1,i),rwork(ic,i),lpm_lrp2)
182  enddo
183 
184  return
185  end
186 c-----------------------------------------------------------------------
187  subroutine lpm_comm_bin_setup
188  include 'SIZE'
189  include 'TOTAL'
190 # include "LPM"
191 
192  integer el_face_num(18),el_edge_num(36),el_corner_num(24),
193  > nfacegp, nedgegp, ncornergp
194  integer ifac(3), icount(3)
195  real d2new(3)
196  logical partl
197 
198  real lpm_xerange(2,3,lpm_lbmax)
199  common /lpm_elementrange/ lpm_xerange
200 
201 c face, edge, and corner number, x,y,z are all inline, so stride=3
202  el_face_num = (/ -1,0,0, 1,0,0, 0,-1,0, 0,1,0, 0,0,-1, 0,0,1 /)
203  el_edge_num = (/ -1,-1,0 , 1,-1,0, 1,1,0 , -1,1,0 ,
204  > 0,-1,-1, 1,0,-1, 0,1,-1, -1,0,-1,
205  > 0,-1,1 , 1,0,1 , 0,1,1 , -1,0,1 /)
206  el_corner_num = (/
207  > -1,-1,-1, 1,-1,-1, 1,1,-1, -1,1,-1,
208  > -1,-1,1, 1,-1,1, 1,1,1, -1,1,1 /)
209 
210  nfacegp = 4 ! number of faces
211  nedgegp = 4 ! number of edges
212  ncornergp = 0 ! number of corners
213 
214  if (if3d) then
215  nfacegp = 6 ! number of faces
216  nedgegp = 12 ! number of edges
217  ncornergp = 8 ! number of corners
218  endif
219 
220  ix = 1
221  iy = 2
222  iz = 3
223 
224  iperiodicx = int(lpm_rparam(8))
225  iperiodicy = int(lpm_rparam(9))
226  iperiodicz = int(lpm_rparam(10))
227 
228  ! compute binb
229  xmin = 1e8
230  ymin = 1e8
231  zmin = 1e8
232  xmax = 0.
233  ymax = 0.
234  zmax = 0.
235  do i=1,lpm_npart
236  rduml = lpm_y(ix,i) - lpm_d2chk(2)
237  rdumr = lpm_y(ix,i) + lpm_d2chk(2)
238  if (rduml .lt. xmin) xmin = rduml
239  if (rdumr .gt. xmax) xmax = rdumr
240 
241  rduml = lpm_y(iy,i) - lpm_d2chk(2)
242  rdumr = lpm_y(iy,i) + lpm_d2chk(2)
243  if (rduml .lt. ymin) ymin = rduml
244  if (rdumr .gt. ymax) ymax = rdumr
245 
246  rduml = lpm_y(iz,i) - lpm_d2chk(2)
247  rdumr = lpm_y(iz,i) + lpm_d2chk(2)
248  if (rduml .lt. zmin) zmin = rduml
249  if (rdumr .gt. zmax) zmax = rdumr
250  enddo
251 
252  lpm_binb(1) = glmin(xmin,1)
253  lpm_binb(2) = glmax(xmax,1)
254  lpm_binb(3) = glmin(ymin,1)
255  lpm_binb(4) = glmax(ymax,1)
256  lpm_binb(5) = 0.0
257  lpm_binb(6) = 0.0
258  if(if3d) lpm_binb(5) = glmin(zmin,1)
259  if(if3d) lpm_binb(6) = glmax(zmax,1)
260 
261  lpm_binb(1) = max(lpm_binb(1),lpm_xdrange(1,1))
262  lpm_binb(2) = min(lpm_binb(2),lpm_xdrange(2,1))
263  lpm_binb(3) = max(lpm_binb(3),lpm_xdrange(1,2))
264  lpm_binb(4) = min(lpm_binb(4),lpm_xdrange(2,2))
265  if(if3d)lpm_binb(5) = max(lpm_binb(5),lpm_xdrange(1,3))
266  if(if3d)lpm_binb(6) = min(lpm_binb(6),lpm_xdrange(2,3))
267 
268  if (iperiodicx .eq. 1) then
269  lpm_binb(1) = lpm_xdrange(1,1)
270  lpm_binb(2) = lpm_xdrange(2,1)
271  endif
272  if (iperiodicy .eq. 1) then
273  lpm_binb(3) = lpm_xdrange(1,2)
274  lpm_binb(4) = lpm_xdrange(2,2)
275  endif
276  if (iperiodicz .eq. 1 .and. if3d) then
277  lpm_binb(5) = lpm_xdrange(1,3)
278  lpm_binb(6) = lpm_xdrange(2,3)
279  endif
280 
281  ifac(1) = 1
282  ifac(2) = 1
283  ifac(3) = 1
284  icount(1) = 0
285  icount(2) = 0
286  icount(3) = 0
287  d2new(1) = lpm_d2chk(2)
288  d2new(2) = lpm_d2chk(2)
289  d2new(3) = lpm_d2chk(2)
290 
291  lpm_ndxgp = floor( (lpm_binb(2) - lpm_binb(1))/d2new(1))
292  lpm_ndygp = floor( (lpm_binb(4) - lpm_binb(3))/d2new(2))
293  lpm_ndzgp = 1
294  if (if3d) lpm_ndzgp = floor( (lpm_binb(6) - lpm_binb(5))/d2new(3))
295 
296 c if (lpm_ndxgp*lpm_ndygp*lpm_ndzgp .gt. np) then
297  if (lpm_ndxgp*lpm_ndygp*lpm_ndzgp .gt. np .or.
298  > int(lpm_rparam(4)) .eq. 1) then
299  nmax = 1000
300  d2chk_save = lpm_d2chk(2)
301 
302  do i=1,nmax
303  do j=0,ndim-1
304  iflg = 0
305  ifac(j+1) = 1 + i
306  d2new(j+1) = (lpm_binb(2+2*j) - lpm_binb(1+2*j))/ifac(j+1)
307  nbb = ifac(1)*ifac(2)*ifac(3)
308  if (int(lpm_rparam(4)) .eq. 0) then
309  if(d2new(j+1) .lt. d2chk_save .or. nbb .gt. np) iflg = 1
310  elseif (int(lpm_rparam(4)) .eq. 1) then
311  if( nbb .gt. np ) iflg = 1
312  endif
313  if (iflg .eq. 1) then
314  icount(j+1) = icount(j+1) + 1
315  ifac(j+1) = ifac(j+1) - icount(j+1)
316  d2new(j+1) = (lpm_binb(2+2*j) -lpm_binb(1+2*j))/ifac(j+1)
317  endif
318  enddo
319  if (icount(1) .gt. 0) then
320  if (icount(2) .gt. 0) then
321  if (icount(3) .gt. 0) then
322  exit
323  endif
324  endif
325  endif
326  enddo
327  endif
328 
329 ! -------------------------------------------------------
330 c SETUP 3D BACKGROUND GRID PARAMETERS FOR GHOST PARTICLES
331 ! -------------------------------------------------------
332  ! how many spacings in each direction
333  lpm_ndxgp = floor( (lpm_binb(2) - lpm_binb(1))/d2new(1))
334  lpm_ndygp = floor( (lpm_binb(4) - lpm_binb(3))/d2new(2))
335  lpm_ndzgp = 1
336  if (if3d) lpm_ndzgp = floor( (lpm_binb(6) - lpm_binb(5))/d2new(3))
337 
338  ! grid spacing for that many spacings
339  lpm_rdxgp = (lpm_binb(2) - lpm_binb(1))/real(lpm_ndxgp)
340  lpm_rdygp = (lpm_binb(4) - lpm_binb(3))/real(lpm_ndygp)
341  lpm_rdzgp = 1.
342  if (if3d) lpm_rdzgp = (lpm_binb(6) - lpm_binb(5))/real(lpm_ndzgp)
343 
344  ninc = 2
345  rxlbin = lpm_binb(1)
346  rxrbin = lpm_binb(2)
347  rylbin = lpm_binb(3)
348  ryrbin = lpm_binb(4)
349  rzlbin = lpm_binb(5)
350  rzrbin = lpm_binb(6)
351  if (iperiodicx .ne. 1) then
352  rxlbin = rxlbin - ninc/2*lpm_rdxgp
353  rxrbin = rxrbin + ninc/2*lpm_rdxgp
354  rxlbin = max(rxlbin,lpm_xdrange(1,1))
355  rxrbin = min(rxrbin,lpm_xdrange(2,1))
356  endif
357  if (iperiodicy .ne. 1) then
358  rylbin = rylbin - ninc/2*lpm_rdygp
359  ryrbin = ryrbin + ninc/2*lpm_rdygp
360  rylbin = max(rylbin,lpm_xdrange(1,2))
361  ryrbin = min(ryrbin,lpm_xdrange(2,2))
362  endif
363  if (iperiodicz .ne. 1) then
364  if (if3d) then
365  rzlbin = rzlbin - ninc/2*lpm_rdzgp
366  rzrbin = rzrbin + ninc/2*lpm_rdzgp
367  rzlbin = max(rzlbin,lpm_xdrange(1,3))
368  rzrbin = min(rzrbin,lpm_xdrange(2,3))
369  endif
370  endif
371 
372  nbin_now = lpm_ndxgp*lpm_ndygp*lpm_ndzgp
373 
374  if (int(lpm_rparam(4)) .eq. 1) return ! only for projection
375 
376 ! see which bins are in which elements
377  lpm_neltb = 0
378  do ie=1,nelt
379  do k=1,nz1
380  do j=1,ny1
381  do i=1,nx1
382  rxval = xm1(i,j,k,ie)
383  ryval = ym1(i,j,k,ie)
384  rzval = 0.
385  if(if3d) rzval = zm1(i,j,k,ie)
386 
387  if (rxval .gt. lpm_binb(2)) goto 1233
388  if (rxval .lt. lpm_binb(1)) goto 1233
389  if (ryval .gt. lpm_binb(4)) goto 1233
390  if (ryval .lt. lpm_binb(3)) goto 1233
391  if (if3d .and. rzval .gt. lpm_binb(6)) goto 1233
392  if (if3d .and. rzval .lt. lpm_binb(5)) goto 1233
393 
394  ii = floor((rxval-lpm_binb(1))/lpm_rdxgp)
395  jj = floor((ryval-lpm_binb(3))/lpm_rdygp)
396  kk = floor((rzval-lpm_binb(5))/lpm_rdzgp)
397  if (.not. if3d) kk = 0
398  if (ii .eq. lpm_ndxgp) ii = lpm_ndxgp - 1
399  if (jj .eq. lpm_ndygp) jj = lpm_ndygp - 1
400  if (kk .eq. lpm_ndzgp) kk = lpm_ndzgp - 1
401  ndum = ii + lpm_ndxgp*jj + lpm_ndxgp*lpm_ndygp*kk
402  nrank = modulo(ndum,np)
403 
404  lpm_neltb = lpm_neltb + 1
405  if(lpm_neltb .gt. lpm_lbmax) then
406  write(6,*) 'increase lbmax',nid,lpm_neltb,lpm_lbmax
407  call exitt
408  endif
409 
410  lpm_er_map(1,lpm_neltb) = ie
411  lpm_er_map(2,lpm_neltb) = nid
412  lpm_er_map(3,lpm_neltb) = ndum
413  lpm_er_map(4,lpm_neltb) = nrank
414  lpm_er_map(5,lpm_neltb) = nrank
415  lpm_er_map(6,lpm_neltb) = nrank
416 
417  if (lpm_neltb .gt. 1) then
418  do il=1,lpm_neltb-1
419  if (lpm_er_map(1,il) .eq. ie) then
420  if (lpm_er_map(4,il) .eq. nrank) then
421  lpm_neltb = lpm_neltb - 1
422  goto 1233
423  endif
424  endif
425  enddo
426  endif
427  1233 continue
428  enddo
429  enddo
430  enddo
431  enddo
432 
433  nxyz = lx1*ly1*lz1
434  do ie=1,lpm_neltb
435  iee = lpm_er_map(1,ie)
436  call copy(lpm_xm1b(1,1,1,1,ie), xm1(1,1,1,iee),nxyz)
437  call copy(lpm_xm1b(1,1,1,2,ie), ym1(1,1,1,iee),nxyz)
438  call copy(lpm_xm1b(1,1,1,3,ie), zm1(1,1,1,iee),nxyz)
439  enddo
440 
441  lpm_neltbb = lpm_neltb
442 
443  do ie=1,lpm_neltbb
444  call icopy(lpm_er_maps(1,ie),lpm_er_map(1,ie),lpm_lrmax)
445  enddo
446 
447 
448  nl = 0
449  nii = lpm_lrmax
450  njj = 6
451  nrr = nxyz*3
452  nkey = 3
453  call fgslib_crystal_tuple_transfer(i_cr_hndl,lpm_neltb,lpm_lbmax
454  > , lpm_er_map,nii,partl,nl,lpm_xm1b,nrr,njj)
455  call fgslib_crystal_tuple_sort (i_cr_hndl,lpm_neltb
456  $ , lpm_er_map,nii,partl,nl,lpm_xm1b,nrr,nkey,1)
457 
458 
459  do ie=1,lpm_neltb
460  do k=1,nz1
461  do j=1,ny1
462  do i=1,nx1
463  rxval = lpm_xm1b(i,j,k,1,ie)
464  ryval = lpm_xm1b(i,j,k,2,ie)
465  rzval = 0.
466  if(if3d) rzval = lpm_xm1b(i,j,k,3,ie)
467 
468  ii = floor((rxval-lpm_binb(1))/lpm_rdxgp)
469  jj = floor((ryval-lpm_binb(3))/lpm_rdygp)
470  kk = floor((rzval-lpm_binb(5))/lpm_rdzgp)
471  if (.not. if3d) kk = 0
472  if (ii .eq. lpm_ndxgp) ii = lpm_ndxgp - 1
473  if (jj .eq. lpm_ndygp) jj = lpm_ndygp - 1
474  if (kk .eq. lpm_ndzgp) kk = lpm_ndzgp - 1
475  ndum = ii + lpm_ndxgp*jj + lpm_ndxgp*lpm_ndygp*kk
476 
477  lpm_modgp(i,j,k,ie,1) = ii
478  lpm_modgp(i,j,k,ie,2) = jj
479  lpm_modgp(i,j,k,ie,3) = kk
480  lpm_modgp(i,j,k,ie,4) = ndum
481 
482  enddo
483  enddo
484  enddo
485  enddo
486 
487  do ie=1,lpm_neltb
488  lpm_xerange(1,1,ie) = vlmin(lpm_xm1b(1,1,1,1,ie),nxyz)
489  lpm_xerange(2,1,ie) = vlmax(lpm_xm1b(1,1,1,1,ie),nxyz)
490  lpm_xerange(1,2,ie) = vlmin(lpm_xm1b(1,1,1,2,ie),nxyz)
491  lpm_xerange(2,2,ie) = vlmax(lpm_xm1b(1,1,1,2,ie),nxyz)
492  lpm_xerange(1,3,ie) = vlmin(lpm_xm1b(1,1,1,3,ie),nxyz)
493  lpm_xerange(2,3,ie) = vlmax(lpm_xm1b(1,1,1,3,ie),nxyz)
494 
495  ilow = floor((lpm_xerange(1,1,ie) - lpm_binb(1))/lpm_rdxgp)
496  ihigh = floor((lpm_xerange(2,1,ie) - lpm_binb(1))/lpm_rdxgp)
497  jlow = floor((lpm_xerange(1,2,ie) - lpm_binb(3))/lpm_rdygp)
498  jhigh = floor((lpm_xerange(2,2,ie) - lpm_binb(3))/lpm_rdygp)
499  klow = floor((lpm_xerange(1,3,ie) - lpm_binb(5))/lpm_rdzgp)
500  khigh = floor((lpm_xerange(2,3,ie) - lpm_binb(5))/lpm_rdzgp)
501  if (.not. if3d) then
502  klow = 0
503  khigh = 0
504  endif
505 
506  lpm_el_map(1,ie) = ilow + lpm_ndxgp*jlow
507  > + lpm_ndxgp*lpm_ndygp*klow
508  lpm_el_map(2,ie) = ihigh + lpm_ndxgp*jhigh
509  > + lpm_ndxgp*lpm_ndygp*khigh
510  lpm_el_map(3,ie) = ilow
511  lpm_el_map(4,ie) = ihigh
512  lpm_el_map(5,ie) = jlow
513  lpm_el_map(6,ie) = jhigh
514  lpm_el_map(7,ie) = klow
515  lpm_el_map(8,ie) = khigh
516  enddo
517 
518  return
519  end
520 c-----------------------------------------------------------------------
521  subroutine lpm_comm_ghost_create
522  include 'SIZE'
523  include 'TOTAL'
524 # include "LPM"
525 
526  character*132 deathmessage
527  real xdlen,ydlen,zdlen,rxdrng(3),rxnew(3)
528  integer iadd(3),gpsave(27)
529  real map(LPM_LRP_PRO)
530 
531  integer el_face_num(18),el_edge_num(36),el_corner_num(24),
532  > nfacegp, nedgegp, ncornergp
533 
534 c face, edge, and corner number, x,y,z are all inline, so stride=3
535  el_face_num = (/ -1,0,0, 1,0,0, 0,-1,0, 0,1,0, 0,0,-1, 0,0,1 /)
536  el_edge_num = (/ -1,-1,0 , 1,-1,0, 1,1,0 , -1,1,0 ,
537  > 0,-1,-1, 1,0,-1, 0,1,-1, -1,0,-1,
538  > 0,-1,1 , 1,0,1 , 0,1,1 , -1,0,1 /)
539  el_corner_num = (/ -1,-1,-1, 1,-1,-1, 1,1,-1, -1,1,-1,
540  > -1,-1,1, 1,-1,1, 1,1,1, -1,1,1 /)
541 
542  nfacegp = 4 ! number of faces
543  nedgegp = 4 ! number of edges
544  ncornergp = 0 ! number of corners
545 
546  if (if3d) then
547  nfacegp = 6 ! number of faces
548  nedgegp = 12 ! number of edges
549  ncornergp = 8 ! number of corners
550  endif
551 
552  iperiodicx = int(lpm_rparam(8))
553  iperiodicy = int(lpm_rparam(9))
554  iperiodicz = int(lpm_rparam(10))
555 
556 ! ------------------------
557 c CREATING GHOST PARTICLES
558 ! ------------------------
559  jx = 1
560  jy = 2
561  jz = 3
562  jdp = int(lpm_rparam(5))
563 
564  xdlen = lpm_binb(2) - lpm_binb(1)
565  ydlen = lpm_binb(4) - lpm_binb(3)
566  zdlen = -1.
567  if (if3d) zdlen = lpm_binb(6) - lpm_binb(5)
568  if (iperiodicx .ne. 1) xdlen = -1
569  if (iperiodicy .ne. 1) ydlen = -1
570  if (iperiodicz .ne. 1) zdlen = -1
571 
572  rxdrng(1) = xdlen
573  rxdrng(2) = ydlen
574  rxdrng(3) = zdlen
575 
576  lpm_npart_gp = 0
577 
578  rfac = 1.0
579 
580  do ip=1,lpm_npart
581 
582  call lpm_project_map(map,lpm_y(1,ip),lpm_ydot(1,ip)
583  > ,lpm_ydotc(1,ip),lpm_rprop(1,ip))
584  lpm_cp_map(1,ip) = lpm_y(jx,ip) ! x coord
585  lpm_cp_map(2,ip) = lpm_y(jy,ip) ! y coord
586  lpm_cp_map(3,ip) = lpm_y(jz,ip) ! z coord
587  lpm_cp_map(4,ip) = lpm_rprop(jdp,ip) ! filter non-dim scale
588  do j=1,lpm_lrp_pro
589  lpm_cp_map(4+j,ip) = map(j)
590  enddo
591 
592  rxval = lpm_cp_map(1,ip)
593  ryval = lpm_cp_map(2,ip)
594  rzval = 0.
595  if(if3d) rzval = lpm_cp_map(3,ip)
596 
597  iip = lpm_iprop(8,ip)
598  jjp = lpm_iprop(9,ip)
599  kkp = lpm_iprop(10,ip)
600 
601  rxl = lpm_binb(1) + lpm_rdxgp*iip
602  rxr = rxl + lpm_rdxgp
603  ryl = lpm_binb(3) + lpm_rdygp*jjp
604  ryr = ryl + lpm_rdygp
605  rzl = 0.0
606  rzr = 0.0
607  if (if3d) then
608  rzl = lpm_binb(5) + lpm_rdzgp*kkp
609  rzr = rzl + lpm_rdzgp
610  endif
611 
612  isave = 0
613 
614  ! faces
615  do ifc=1,nfacegp
616  ist = (ifc-1)*3
617  ii1 = iip + el_face_num(ist+1)
618  jj1 = jjp + el_face_num(ist+2)
619  kk1 = kkp + el_face_num(ist+3)
620 
621  iig = ii1
622  jjg = jj1
623  kkg = kk1
624 
625  distchk = 0.0
626  dist = 0.0
627  if (ii1-iip .ne. 0) then
628  distchk = distchk + (rfac*lpm_d2chk(2))**2
629  if (ii1-iip .lt. 0) dist = dist +(rxval - rxl)**2
630  if (ii1-iip .gt. 0) dist = dist +(rxval - rxr)**2
631  endif
632  if (jj1-jjp .ne. 0) then
633  distchk = distchk + (rfac*lpm_d2chk(2))**2
634  if (jj1-jjp .lt. 0) dist = dist +(ryval - ryl)**2
635  if (jj1-jjp .gt. 0) dist = dist +(ryval - ryr)**2
636  endif
637  if (if3d) then
638  if (kk1-kkp .ne. 0) then
639  distchk = distchk + (rfac*lpm_d2chk(2))**2
640  if (kk1-kkp .lt. 0) dist = dist +(rzval - rzl)**2
641  if (kk1-kkp .gt. 0) dist = dist +(rzval - rzr)**2
642  endif
643  endif
644  distchk = sqrt(distchk)
645  dist = sqrt(dist)
646  if (dist .gt. distchk) cycle
647 
648  iflgx = 0
649  iflgy = 0
650  iflgz = 0
651  ! periodic if out of domain - add some ifsss
652  if (iig .lt. 0 .or. iig .gt. lpm_ndxgp-1) then
653  iflgx = 1
654  iig =modulo(iig,lpm_ndxgp)
655  if (iperiodicx .ne. 1) cycle
656  endif
657  if (jjg .lt. 0 .or. jjg .gt. lpm_ndygp-1) then
658  iflgy = 1
659  jjg =modulo(jjg,lpm_ndygp)
660  if (iperiodicy .ne. 1) cycle
661  endif
662  if (kkg .lt. 0 .or. kkg .gt. lpm_ndzgp-1) then
663  iflgz = 1
664  kkg =modulo(kkg,lpm_ndzgp)
665  if (iperiodicz .ne. 1) cycle
666  endif
667 
668  iflgsum = iflgx + iflgy + iflgz
669  ndumn = iig + lpm_ndxgp*jjg + lpm_ndxgp*lpm_ndygp*kkg
670  nrank = modulo(ndumn,np)
671 
672  if (nrank .eq. nid .and. iflgsum .eq. 0) cycle
673 
674  do i=1,isave
675  if (gpsave(i) .eq. nrank .and. iflgsum .eq.0) goto 111
676  enddo
677  isave = isave + 1
678  gpsave(isave) = nrank
679 
680  ibctype = iflgx+iflgy+iflgz
681 
682  rxnew(1) = rxval
683  rxnew(2) = ryval
684  rxnew(3) = rzval
685 
686  iadd(1) = ii1
687  iadd(2) = jj1
688  iadd(3) = kk1
689 
690  call lpm_comm_check_periodic_gp(rxnew,rxdrng,iadd)
691 
692  lpm_npart_gp = lpm_npart_gp + 1
693  lpm_iprop_gp(1,lpm_npart_gp) = nrank
694  lpm_iprop_gp(2,lpm_npart_gp) = iig
695  lpm_iprop_gp(3,lpm_npart_gp) = jjg
696  lpm_iprop_gp(4,lpm_npart_gp) = kkg
697  lpm_iprop_gp(5,lpm_npart_gp) = ndumn
698 
699  lpm_rprop_gp(1,lpm_npart_gp) = rxnew(1)
700  lpm_rprop_gp(2,lpm_npart_gp) = rxnew(2)
701  lpm_rprop_gp(3,lpm_npart_gp) = rxnew(3)
702  do k=4,lpm_lrp_gp
703  lpm_rprop_gp(k,lpm_npart_gp) = lpm_cp_map(k,ip)
704  enddo
705  111 continue
706  enddo
707 
708  ! edges
709  do ifc=1,nedgegp
710  ist = (ifc-1)*3
711  ii1 = iip + el_edge_num(ist+1)
712  jj1 = jjp + el_edge_num(ist+2)
713  kk1 = kkp + el_edge_num(ist+3)
714 
715  iig = ii1
716  jjg = jj1
717  kkg = kk1
718 
719  distchk = 0.0
720  dist = 0.0
721  if (ii1-iip .ne. 0) then
722  distchk = distchk + (rfac*lpm_d2chk(2))**2
723  if (ii1-iip .lt. 0) dist = dist +(rxval - rxl)**2
724  if (ii1-iip .gt. 0) dist = dist +(rxval - rxr)**2
725  endif
726  if (jj1-jjp .ne. 0) then
727  distchk = distchk + (rfac*lpm_d2chk(2))**2
728  if (jj1-jjp .lt. 0) dist = dist +(ryval - ryl)**2
729  if (jj1-jjp .gt. 0) dist = dist +(ryval - ryr)**2
730  endif
731  if (if3d) then
732  if (kk1-kkp .ne. 0) then
733  distchk = distchk + (rfac*lpm_d2chk(2))**2
734  if (kk1-kkp .lt. 0) dist = dist +(rzval - rzl)**2
735  if (kk1-kkp .gt. 0) dist = dist +(rzval - rzr)**2
736  endif
737  endif
738  distchk = sqrt(distchk)
739  dist = sqrt(dist)
740  if (dist .gt. distchk) cycle
741 
742  iflgx = 0
743  iflgy = 0
744  iflgz = 0
745  ! periodic if out of domain - add some ifsss
746  if (iig .lt. 0 .or. iig .gt. lpm_ndxgp-1) then
747  iflgx = 1
748  iig =modulo(iig,lpm_ndxgp)
749  if (iperiodicx .ne. 1) cycle
750  endif
751  if (jjg .lt. 0 .or. jjg .gt. lpm_ndygp-1) then
752  iflgy = 1
753  jjg =modulo(jjg,lpm_ndygp)
754  if (iperiodicy .ne. 1) cycle
755  endif
756  if (kkg .lt. 0 .or. kkg .gt. lpm_ndzgp-1) then
757  iflgz = 1
758  kkg =modulo(kkg,lpm_ndzgp)
759  if (iperiodicz .ne. 1) cycle
760  endif
761 
762  iflgsum = iflgx + iflgy + iflgz
763  ndumn = iig + lpm_ndxgp*jjg + lpm_ndxgp*lpm_ndygp*kkg
764  nrank = modulo(ndumn,np)
765 
766  if (nrank .eq. nid .and. iflgsum .eq. 0) cycle
767 
768  do i=1,isave
769  if (gpsave(i) .eq. nrank .and. iflgsum .eq.0) goto 222
770  enddo
771  isave = isave + 1
772  gpsave(isave) = nrank
773 
774  ibctype = iflgx+iflgy+iflgz
775 
776  rxnew(1) = rxval
777  rxnew(2) = ryval
778  rxnew(3) = rzval
779 
780  iadd(1) = ii1
781  iadd(2) = jj1
782  iadd(3) = kk1
783 
784  call lpm_comm_check_periodic_gp(rxnew,rxdrng,iadd)
785 
786  lpm_npart_gp = lpm_npart_gp + 1
787  lpm_iprop_gp(1,lpm_npart_gp) = nrank
788  lpm_iprop_gp(2,lpm_npart_gp) = iig
789  lpm_iprop_gp(3,lpm_npart_gp) = jjg
790  lpm_iprop_gp(4,lpm_npart_gp) = kkg
791  lpm_iprop_gp(5,lpm_npart_gp) = ndumn
792 
793  lpm_rprop_gp(1,lpm_npart_gp) = rxnew(1)
794  lpm_rprop_gp(2,lpm_npart_gp) = rxnew(2)
795  lpm_rprop_gp(3,lpm_npart_gp) = rxnew(3)
796  do k=4,lpm_lrp_gp
797  lpm_rprop_gp(k,lpm_npart_gp) = lpm_cp_map(k,ip)
798  enddo
799  222 continue
800  enddo
801 
802  ! corners
803  do ifc=1,ncornergp
804  ist = (ifc-1)*3
805  ii1 = iip + el_corner_num(ist+1)
806  jj1 = jjp + el_corner_num(ist+2)
807  kk1 = kkp + el_corner_num(ist+3)
808 
809  iig = ii1
810  jjg = jj1
811  kkg = kk1
812 
813  distchk = 0.0
814  dist = 0.0
815  if (ii1-iip .ne. 0) then
816  distchk = distchk + (rfac*lpm_d2chk(2))**2
817  if (ii1-iip .lt. 0) dist = dist +(rxval - rxl)**2
818  if (ii1-iip .gt. 0) dist = dist +(rxval - rxr)**2
819  endif
820  if (jj1-jjp .ne. 0) then
821  distchk = distchk + (rfac*lpm_d2chk(2))**2
822  if (jj1-jjp .lt. 0) dist = dist +(ryval - ryl)**2
823  if (jj1-jjp .gt. 0) dist = dist +(ryval - ryr)**2
824  endif
825  if (if3d) then
826  if (kk1-kkp .ne. 0) then
827  distchk = distchk + (rfac*lpm_d2chk(2))**2
828  if (kk1-kkp .lt. 0) dist = dist +(rzval - rzl)**2
829  if (kk1-kkp .gt. 0) dist = dist +(rzval - rzr)**2
830  endif
831  endif
832  distchk = sqrt(distchk)
833  dist = sqrt(dist)
834  if (dist .gt. distchk) cycle
835 
836  iflgx = 0
837  iflgy = 0
838  iflgz = 0
839  ! periodic if out of domain - add some ifsss
840  if (iig .lt. 0 .or. iig .gt. lpm_ndxgp-1) then
841  iflgx = 1
842  iig =modulo(iig,lpm_ndxgp)
843  if (iperiodicx .ne. 1) cycle
844  endif
845  if (jjg .lt. 0 .or. jjg .gt. lpm_ndygp-1) then
846  iflgy = 1
847  jjg =modulo(jjg,lpm_ndygp)
848  if (iperiodicy .ne. 1) cycle
849  endif
850  if (kkg .lt. 0 .or. kkg .gt. lpm_ndzgp-1) then
851  iflgz = 1
852  kkg =modulo(kkg,lpm_ndzgp)
853  if (iperiodicz .ne. 1) cycle
854  endif
855 
856  iflgsum = iflgx + iflgy + iflgz
857  ndumn = iig + lpm_ndxgp*jjg + lpm_ndxgp*lpm_ndygp*kkg
858  nrank = modulo(ndumn,np)
859 
860  if (nrank .eq. nid .and. iflgsum .eq. 0) cycle
861 
862  do i=1,isave
863  if (gpsave(i) .eq. nrank .and. iflgsum .eq.0) goto 333
864  enddo
865  isave = isave + 1
866  gpsave(isave) = nrank
867 
868  ibctype = iflgx+iflgy+iflgz
869 
870  rxnew(1) = rxval
871  rxnew(2) = ryval
872  rxnew(3) = rzval
873 
874  iadd(1) = ii1
875  iadd(2) = jj1
876  iadd(3) = kk1
877 
878  call lpm_comm_check_periodic_gp(rxnew,rxdrng,iadd)
879 
880  lpm_npart_gp = lpm_npart_gp + 1
881  lpm_iprop_gp(1,lpm_npart_gp) = nrank
882  lpm_iprop_gp(2,lpm_npart_gp) = iig
883  lpm_iprop_gp(3,lpm_npart_gp) = jjg
884  lpm_iprop_gp(4,lpm_npart_gp) = kkg
885  lpm_iprop_gp(5,lpm_npart_gp) = ndumn
886 
887  lpm_rprop_gp(1,lpm_npart_gp) = rxnew(1)
888  lpm_rprop_gp(2,lpm_npart_gp) = rxnew(2)
889  lpm_rprop_gp(3,lpm_npart_gp) = rxnew(3)
890  do k=4,lpm_lrp_gp
891  lpm_rprop_gp(k,lpm_npart_gp) = lpm_cp_map(k,ip)
892  enddo
893  333 continue
894  enddo
895 
896  enddo
897 
898  return
899  end
900 c----------------------------------------------------------------------
901  subroutine lpm_comm_check_periodic_gp(rxnew,rxdrng,iadd)
902  include 'SIZE'
903  include 'TOTAL'
904 # include "LPM"
905 
906  real rxnew(3), rxdrng(3)
907  integer iadd(3), irett(3), ntype, ntypel(7)
908 
909  xloc = rxnew(1)
910  yloc = rxnew(2)
911  zloc = rxnew(3)
912 
913  xdlen = rxdrng(1)
914  ydlen = rxdrng(2)
915  zdlen = rxdrng(3)
916 
917  ii = iadd(1)
918  jj = iadd(2)
919  kk = iadd(3)
920 
921  irett(1) = 0
922  irett(2) = 0
923  irett(3) = 0
924 
925  if (xdlen .gt. 0 ) then
926  if (ii .ge. lpm_ndxgp) then
927  xloc = xloc - xdlen
928  irett(1) = 1
929  goto 123
930  endif
931  endif
932  if (xdlen .gt. 0 ) then
933  if (ii .lt. 0) then
934  xloc = xloc + xdlen
935  irett(1) = 1
936  goto 123
937  endif
938  endif
939 
940  123 continue
941  if (ydlen .gt. 0 ) then
942  if (jj .ge. lpm_ndygp) then
943  yloc = yloc - ydlen
944  irett(2) = 1
945  goto 124
946  endif
947  endif
948  if (ydlen .gt. 0 ) then
949  if (jj .lt. 0) then
950  yloc = yloc + ydlen
951  irett(2) = 1
952  goto 124
953  endif
954  endif
955  124 continue
956 
957  if (if3d) then
958  if (zdlen .gt. 0 ) then
959  if (kk .ge. lpm_ndzgp) then
960  zloc = zloc - zdlen
961  irett(3) = 1
962  goto 125
963  endif
964  endif
965  if (zdlen .gt. 0 ) then
966  if (kk .lt. 0) then
967  zloc = zloc + zdlen
968  irett(3) = 1
969  goto 125
970  endif
971  endif
972  endif
973  125 continue
974 
975  rxnew(1) = xloc
976  rxnew(2) = yloc
977  rxnew(3) = zloc
978 
979  return
980  end
981 c----------------------------------------------------------------------
982  subroutine lpm_comm_ghost_send
983  include 'SIZE'
984 # include "LPM"
985 
986  logical partl
987 
988  ! send ghost particles
989  call fgslib_crystal_tuple_transfer(i_cr_hndl
990  > ,lpm_npart_gp,lpm_lpart_gp
991  > ,lpm_iprop_gp,lpm_lip_gp
992  > ,partl,0
993  > ,lpm_rprop_gp,lpm_lrp_gp
994  > ,1)
995 
996  return
997  end
998 c----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
interp_setup
subroutine interp_setup(ih, tolin, nmsh, nelm)
Definition: interp.f:5
lpm_comm_bin_setup
subroutine lpm_comm_bin_setup
Definition: lpm_comm.f:188
nid_glcount
integer function nid_glcount(a, is, n)
Definition: lpm_comm.f:81
lpm_comm_ghost_create
subroutine lpm_comm_ghost_create
Definition: lpm_comm.f:522
lpm_comm_ghost_send
subroutine lpm_comm_ghost_send
Definition: lpm_comm.f:983
lpm_comm_setup
subroutine lpm_comm_setup
Definition: lpm_comm.f:3
lpm_comm_crystal
subroutine lpm_comm_crystal
Definition: lpm_comm.f:136
lpm_comm_check_periodic_gp
subroutine lpm_comm_check_periodic_gp(rxnew, rxdrng, iadd)
Definition: lpm_comm.f:902
lpm_comm_findpts
subroutine lpm_comm_findpts
Definition: lpm_comm.f:19
glmin
function glmin(a, n)
Definition: math.f:973
icopy
subroutine icopy(a, b, n)
Definition: math.f:289
iglsum
function iglsum(a, n)
Definition: math.f:926
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
iglmax
function iglmax(a, n)
Definition: math.f:913
vlmin
real function vlmin(vec, n)
Definition: math.f:357
glmax
function glmax(a, n)
Definition: math.f:960