source: CIVL/examples/fortran/nek5000/core/map2.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c

  • Property mode set to 100755
File size: 23.8 KB
Line 
1c-----------------------------------------------------------------------
2 subroutine mapelpr()
3 include 'SIZE'
4 include 'INPUT'
5 include 'PARALLEL'
6 include 'SCRCT'
7 include 'SOLN'
8 include 'TSTEP'
9 include 'CTIMER'
10c
11 logical ifverbm
12c
13 if (nio.eq.0) then
14 write(6,12) 'nelgt/nelgv/lelt:',nelgt,nelgv,lelt
15 write(6,12) 'lx1/lx2/lx3/lxd: ',lx1,lx2,lx3,lxd
16 12 format(1X,A,4I12)
17 write(6,*)
18 endif
19
20 etime0 = dnekclock_sync()
21 if(nio.eq.0) write(6,'(A)') ' partioning elements to MPI ranks'
22
23 MFIELD=2
24 IF (IFFLOW) MFIELD=1
25 IF (IFMVBD) MFIELD=0
26
27c Set up TEMPORARY value for NFIELD - NFLDT
28 NFLDT = 1
29 IF (IFHEAT) NFLDT = 2 + NPSCAL
30c
31c Distributed memory processor mapping
32 IF (NP.GT.NELGT) THEN
33 IF(NID.EQ.0) THEN
34 WRITE(6,1000) NP,NELGT
35 1000 FORMAT(2X,'ABORT: Too many processors (',I8
36 $ ,') for to few elements (',I12,').'
37 $ ,/,2X,'ABORTING IN MAPELPR.')
38 ENDIF
39 call exitt
40 ENDIF
41
42 call set_proc_map()
43
44 if (nelt.gt.lelt) then
45 call exitti('nelt > lelt, increase lelt!$',nelt)
46 endif
47
48 DO 1200 IFIELD=MFIELD,NFLDT
49 IF (IFTMSH(IFIELD)) THEN
50 NELG(IFIELD) = NELGT
51 ELSE
52 NELG(IFIELD) = NELGV
53 ENDIF
54 1200 CONTINUE
55
56C Output the processor-element map:
57 ifverbm=.false.
58 if (loglevel .gt. 2) ifverbm=.true.
59
60 if(ifverbm) then
61 idum = 1
62 if(nid.eq.0) then
63 N8 = min(8,nelt)
64 write(6 ,1310) node-1,(lglel(ie),ie=1,n8)
65 if (NELT.GT.8) write(6 ,1315) (lglel(ie),ie=9,NELT)
66 DO inid=1,NP-1
67 mtype = inid
68 call csend(mtype,idum,4,inid,0) ! handshake
69 call crecv(mtype,inelt,4) ! nelt of other cpus
70 N8 = min(8,inelt)
71 ENDDO
72 1310 FORMAT(' RANK',I6,' IEG',8I8)
73 1315 FORMAT(' ',6X,' ',8I8)
74 else
75 mtype = nid
76 call crecv(mtype,idum,4) ! hand-shake
77 call csend(mtype,nelt,4,0,0) ! nelt
78 if (loglevel .gt. 2) then
79 N8 = min(8,nelt)
80 write(6 ,1310) node-1,(lglel(ie),ie=1,n8)
81 if (NELT.GT.8) write(6 ,1315) (lglel(ie),ie=9,NELT)
82 endif
83 endif
84 endif
85
86 dtmp = dnekclock_sync() - etime0
87 if(nio.eq.0) then
88 write(6,*) ' '
89 write(6,'(A,g13.5,A,/)') ' done :: partioning ',dtmp,' sec'
90 endif
91
92 return
93 end
94c-----------------------------------------------------------------------
95 subroutine outmati(u,m,n,name6)
96 integer u(m,n)
97 character*6 name6
98 common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
99c
100c Print out copies of a global matrix
101c
102 do mid=0,np-1
103 call nekgsync
104 if (mid.eq.nid) then
105 n20 = min(n,20)
106 write(6,1) nid,m,n,name6
107 1 format(//,3i6,' Matrix:',2x,a6,/)
108 do i=1,m
109 write(6,2) nid,name6,(u(i,j),j=1,n20)
110 enddo
111 2 format(i3,1x,a6,20i6)
112 endif
113 call nekgsync
114 enddo
115 return
116 end
117c-----------------------------------------------------------------------
118 subroutine get_map
119
120 call get_vert
121
122 return
123 end
124c-----------------------------------------------------------------------
125 subroutine get_vert
126c
127c Distribute and assign partitions using the .map file
128c
129 include 'SIZE'
130 include 'TOTAL'
131
132 parameter(mdw=2+2**ldim)
133 parameter(ndw=7*lx1*ly1*lz1*lelv/mdw)
134 common /scrns/ wk(mdw,ndw)
135 integer wk
136
137 common /ivrtx/ vertex ((2**ldim),lelt)
138 integer vertex
139
140 integer icalld
141 save icalld
142 data icalld /0/
143
144 if(icalld.gt.0) return
145 icalld = 1
146
147 nv = 2**ldim
148 call get_vert_map(vertex,nv,wk,mdw,ndw)
149
150 return
151 end
152c-----------------------------------------------------------------------
153 subroutine get_vert_map(vertex,nlv,wk,mdw,ndw)
154
155 include 'SIZE'
156 include 'TOTAL'
157
158 integer vertex(nlv,1)
159 integer wk(mdw*ndw)
160
161 common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
162
163 integer ibuf(2)
164 integer hrsb
165
166 integer*8 eid8(lelt), vtx8(lelt*2**ldim)
167 integer iwork(lelt)
168 common /SCRCG/ xyz(ldim*lelt*2**ldim)
169 common /ctmp0/ eid8, vtx8, iwork
170
171 integer tt,cnt,nrank,ierr
172
173 integer opt_parrsb(3), opt_parmetis(10)
174 logical ifbswap
175
176 integer start
177
178#if !defined(PARRSB) && !defined(PARMETIS)
179#if defined(DPROCMAP)
180 call exitti('DPROCMAP requires PARRSB or PARMETIS!$',0)
181#else
182 call read_map(vertex,nlv,wk,mdw,ndw)
183 return
184#endif
185#endif
186
187#if defined(PARRSB) || defined(PARMETIS)
188 neli = nelt
189 call get_con(wk,size(wk),neli,nlv)
190
191c fluid elements
192 j = 0
193 ii = 0
194 cnt= 0
195 do i = 1,neli
196 if (wk(ii+1) .le. nelgv) then
197 j = j + 1
198 eid8(j) = wk(ii+1)
199 call icopy48(vtx8((j-1)*nlv+1),wk(ii+2),nlv)
200
201 do tt=1,nlv
202 xyz(cnt+1)=xc(tt,i)
203 xyz(cnt+2)=yc(tt,i)
204 if(ldim.eq.3) then
205 xyz(cnt+3)=zc(tt,i)
206 cnt=cnt+3
207 else
208 cnt=cnt+2
209 endif
210 enddo
211 endif
212 ii = ii + (nlv+1)
213 enddo
214 neliv = j
215
216 nel = neliv
217 call fpartMesh(eid8,vtx8,xyz,lelt,nel,nlv,nekcomm,
218 $ meshPartitioner,ierr)
219 call err_chk(ierr,'partMesh fluid failed!$')
220
221 nelv = nel
222 nelt = nelv
223 ierr = 0
224 if (nelv .gt. lelv) ierr = 1
225 call err_chk(ierr,'nelv > lelv!$')
226
227 do i = 1,nelv
228 lglel(i) = eid8(i)
229 enddo
230 call isort(lglel,iwork,nelv)
231 do i = 1,nelv
232 call icopy84(vertex(1,i),vtx8((iwork(i)-1)*nlv+1),nlv)
233 enddo
234
235 cnt=0
236c solid elements
237 if (nelgt.ne.nelgv) then
238 j = 0
239 ii = 0
240 do i = 1,neli
241 if (wk(ii+1) .gt. nelgv) then
242 j = j + 1
243 eid8(j) = wk(ii+1)
244 call icopy48(vtx8((j-1)*nlv+1),wk(ii+2),nlv)
245
246 do tt=1,nlv
247 xyz(cnt+1)=xc(tt,i)
248 xyz(cnt+2)=yc(tt,i)
249 if(ldim.eq.3) then
250 xyz(cnt+3)=zc(tt,i)
251 cnt=cnt+3
252 else
253 cnt=cnt+2
254 endif
255 enddo
256 endif
257 ii = ii + (nlv+1)
258 enddo
259 nelit = j
260
261 nel = nelit
262 call fpartMesh(eid8,vtx8,xyz,lelt,nel,nlv,nekcomm,
263 $ meshPartitioner,ierr)
264 call err_chk(ierr,'partMesh solid failed!$')
265
266 nelt = nelv + nel
267 ierr = 0
268 if (nelt .gt. lelt) ierr = 1
269 call err_chk(ierr,'nelt > lelt!$')
270
271 do i = 1,nel
272 lglel(nelv+i) = eid8(i)
273 enddo
274 call isort(lglel(nelv+1),iwork,nel) ! sort locally by global element id
275 do i = 1,nel
276 call icopy84(vertex(1,nelv+i),vtx8((iwork(i)-1)*nlv+1),nlv)
277 enddo
278 endif
279
280#ifdef DPROCMAP
281 do i = 1,nelt
282 ieg = lglel(i)
283 if (ieg.lt.1 .or. ieg.gt.nelgt)
284 $ call exitti('invalid ieg!$',ieg)
285 ibuf(1) = i
286 ibuf(2) = nid
287 call dProcmapPut(ibuf,2,0,ieg)
288 enddo
289#else
290 call izero(gllnid,nelgt)
291 do i = 1,nelt
292 ieg = lglel(i)
293 gllnid(ieg) = nid
294 enddo
295 npass = 1 + nelgt/lelt
296 k=1
297 do ipass = 1,npass
298 m = nelgt - k + 1
299 m = min(m,lelt)
300 if (m.gt.0) call igop(gllnid(k),iwork,'+ ',m)
301 k = k+m
302 enddo
303#endif
304
305#endif
306
307 return
308 end
309c-----------------------------------------------------------------------
310 subroutine get_con(wk,nwk,nelr,nv)
311
312 include 'SIZE'
313 include 'INPUT'
314 include 'PARALLEL'
315
316 integer nwk,nelr,nv
317 integer wk(nwk)
318
319 logical ifbswap,if_byte_swap_test
320 logical ifco2, ifcon
321
322 character*132 confle
323 character*1 confle1(132)
324 equivalence (confle,confle1)
325
326 character*132 hdr
327 character*5 version
328 real*4 test
329
330 integer ierr,nvi
331 integer*8 nelgti,nelgvi
332 integer*8 offs, offs0
333
334 ierr = 0
335
336 ifco2 = .false.
337 ifmpiio = .true.
338#ifdef NOMPIIO
339 ifmpiio = .false.
340#endif
341
342 if (nid.eq.0) then
343 lfname = ltrunc(reafle,132) - 4
344 call blank (confle,132)
345 call chcopy(confle,reafle,lfname)
346 call chcopy(confle1(lfname+1),'.con',4)
347 inquire(file=confle, exist=ifcon)
348
349 if (.not.ifcon) then
350 call chcopy(confle1(lfname+1),'.co2',4)
351 inquire(file=confle, exist=ifco2)
352 endif
353
354 if(.not.ifcon .and. .not.ifco2) ierr = 1
355 endif
356
357 call bcast(ierr,sizeof(ierr))
358 if(ierr.ne.0) goto 50
359
360 call bcast(confle,sizeof(confle))
361 if(nid.eq.0) write(6,'(A,A)') ' reading ', confle
362c call err_chk(ierr,' Cannot find con file!$')
363 call bcast(ifco2,lsize)
364 ierr = 0
365
366 ! read header
367 if (nid.eq.0) then
368 if (ifco2) then
369 call byte_open(confle,ierr)
370 if(ierr.ne.0) goto 100
371
372 call blank(hdr,sizeof(hdr))
373 call byte_read(hdr,sizeof(hdr)/4,ierr)
374 if(ierr.ne.0) goto 100
375
376 read (hdr,*) version,nelgti,nelgvi,nvi
377c 1 format(a5,2i12,i2)
378
379 call byte_read(test,1,ierr)
380 if(ierr.ne.0) goto 100
381 ifbswap = if_byte_swap_test(test,ierr)
382 if(ierr.ne.0) goto 100
383 endif
384 endif
385 call bcast(nelgti,sizeof(nelgti))
386 call bcast(nelgvi,sizeof(nelgvi))
387 call bcast(nvi,sizeof(nvi))
388 call bcast(ifbswap,sizeof(ifbswap))
389
390 if (nvi .ne. nv)
391 $ call exitti('Number of vertices do not match!$',0)
392 if (nelgti .ne. nelgt)
393 $ call exitti('nelgt for mesh/con differs!$',0)
394 if (nelgvi .ne. nelgv)
395 $ call exitti('nelgt for mesh/con differs!$',0)
396
397 if (ifco2 .and. ifmpiio) then
398 if (nid.eq.0) call byte_close(ierr)
399 call byte_open_mpi(confle,ifh,.true.,ierr)
400 offs0 = sizeof(hdr) + sizeof(test)
401
402 call lim_chk(nelr*(nvi+1),nwk,'nelr ','nwk ','read_con ')
403
404 nelBr = igl_running_sum(nelr) - nelr
405 offs = offs0 + int(nelBr,8)*(nvi+1)*ISIZE
406
407 call byte_set_view(offs,ifh)
408 call byte_read_mpi(wk,(nvi+1)*nelr,-1,ifh,ierr)
409 call err_chk(ierr,' Error while reading con file!$')
410 call byte_close_mpi(ifh,ierr)
411 if (ifbswap) call byte_reverse(wk,(nvi+1)*nelr,ierr)
412 endif
413
414 return
415
416 50 continue
417
418#if defined(PARRSB)
419 call find_con(wk,nwk,ierr)
420 return
421#endif
422
423 100 continue
424 call err_chk(ierr,'Error opening/reading con file$')
425 return
426
427 end
428c-----------------------------------------------------------------------
429#if defined(PARRSB)
430 subroutine find_con(wk,nwk,ierr)
431
432 include 'SIZE'
433 include 'INPUT'
434 include 'PARALLEL'
435
436 integer nwk,ierr
437 integer wk(nwk)
438
439 common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
440 common /SCRCG/ xyz(ldim*(2**ldim)*lelt)
441
442 integer*8 eid8(4*lelt),vtx8(lelt*(2**ldim+1))
443 common /ctmp0/ eid8, vtx8, iwork
444
445 integer npf,nv,nf,i,j,k,verbose
446 integer*8 start
447 real*8 tol
448
449 nv=2**ndim
450 nf=2*ndim
451
452 k=0
453 if(ndim.eq.3) then
454 do i=1,nelt
455 do j=1,nv
456 xyz(k+1)=xc(j,i)
457 xyz(k+2)=yc(j,i)
458 xyz(k+3)=zc(j,i)
459 k=k+3
460 enddo
461 enddo
462 else
463 do i=1,nelt
464 do j=1,nv
465 xyz(k+1)=xc(j,i)
466 xyz(k+2)=yc(j,i)
467 k=k+2
468 enddo
469 enddo
470 endif
471
472 start=igl_running_sum(nelt)-nelt
473
474 !calculate number of periodic pairs
475 npf=0
476 do i=1,nelt
477 do j=1,nf
478 if(cbc(j,i,1).eq.'P ') then
479 eid8(4*npf+1)=i+start
480 eid8(4*npf+2)=j
481 eid8(4*npf+3)=bc(1,j,i,1)
482 eid8(4*npf+4)=bc(2,j,i,1)
483 npf=npf+1
484 endif
485 enddo
486 enddo
487
488 ierr = 0
489 toli = 5e-3
490 verbose = 1
491 call fparrsb_findConnectivity(vtx8,xyz,nelt,ndim,
492 $ eid8,npf,tol,nekcomm,verbose,ierr)
493 call err_chk(ierr,' findConnectivity failed!$')
494
495 do i=1,nelt*(nv+1)
496 wk(i)=vtx8(i)
497 enddo
498
499 end
500#endif
501c-----------------------------------------------------------------------
502 subroutine set_proc_map()
503C
504C Compute element to processor distribution according to (weighted)
505C physical distribution in an attempt to minimize exposed number of
506C element interfaces.
507C
508 include 'SIZE'
509 include 'PARALLEL'
510 include 'INPUT'
511 include 'DPROCMAP'
512
513 integer ibuf(2)
514 logical ifbswap
515 logical ifre2
516
517 integer iwork(lelt)
518
519 if(nid.eq.0) inquire(file=re2fle, exist=ifre2)
520 call bcast(ifre2,lsize)
521 if(.not.ifre2) then
522 call set_proc_map_legacy
523 return
524 endif
525
526 call read_re2_hdr(ifbswap, .false.)
527 nelt = nelgt/np
528 do i = 1,mod(nelgt,np)
529 if (np-i.eq.nid) nelt = nelt + 1
530 enddo
531
532 if (nelt .gt. lelt)
533 $ call exitti('nelt > lelt!$',nelt)
534
535 ! setup gllnid + gllel
536#if defined(DPROCMAP)
537 call dProcmapInit()
538 dProcmapCache = .false.
539#endif
540 nelB = igl_running_sum(nelt) - nelt
541 do i = 1,nelt
542 ieg = nelB + i
543 lglel(i) = ieg
544 if (ieg.lt.1 .or. ieg.gt.nelgt)
545 $ call exitti('invalid ieg!$',ieg)
546#if defined(DPROCMAP)
547 ibuf(1) = i
548 ibuf(2) = nid
549 call dProcmapPut(ibuf,2,0,ieg)
550#else
551 gllnid(ieg) = nid
552 gllel(ieg) = i
553#endif
554 enddo
555#if !defined(DPROCMAP)
556 npass = 1 + nelgt/lelt
557 k=1
558 do ipass = 1,npass
559 m = nelgt - k + 1
560 m = min(m,lelt)
561 if (m.gt.0) call igop(gllnid(k),iwork,'+ ',m)
562 if (m.gt.0) call igop(gllel(k) ,iwork,'+ ',m)
563 k = k+m
564 enddo
565#endif
566
567 ! read coord for RCB and/or connectivity
568 call read_re2_data(ifbswap, .true., .false., .true.)
569
570 ! get element-proc mapping
571 call get_map()
572
573 itmp = gllnid(0) ! reset last element cache
574 itmp = gllel(0) ! reset last element cache
575 dProcmapCache = .true.
576
577#if !defined(DPROCMAP)
578 IEL=0
579 CALL IZERO(GLLEL,NELGT)
580 DO IEG=1,NELGT
581 IF (GLLNID(IEG).EQ.NID) THEN
582 IEL = IEL + 1
583 GLLEL(IEG)=IEL
584 NELT = IEL
585 if (ieg.le.nelgv) NELV = IEL
586 ENDIF
587c write(6,*) 'map2 ieg:',ieg,nelv,nelt,nelgv,nelgt
588 ENDDO
589 npass = 1 + nelgt/lelt
590 k=1
591 do ipass = 1,npass
592 m = nelgt - k + 1
593 m = min(m,lelt)
594 if (m.gt.0) call igop(gllel(k),iwork,'+ ',m)
595 k = k+m
596 enddo
597
598 do ieg=1,nelgt
599 mid =gllnid(ieg)
600 ie =gllel (ieg)
601 if (mid.eq.nid) lglel(ie)=ieg
602 enddo
603#endif
604
605 return
606 end
607c-----------------------------------------------------------------------
608 subroutine set_proc_map_legacy()
609C
610C Compute element to processor distribution according to (weighted)
611C physical distribution in an attempt to minimize exposed number of
612C element interfaces.
613C
614 include 'SIZE'
615 include 'INPUT'
616 include 'PARALLEL'
617 include 'SOLN'
618 include 'SCRCT'
619 include 'TSTEP'
620 common /ctmp0/ iwork(lelt)
621
622 REAL*8 dnekclock,t0
623
624#if defined(PARRSB) || defined(PARMETIS) || defined(DPROCMAP)
625 call exitti(' DPROCMAP/PARRSB not supported for rea files$',0)
626#else
627 t0 = dnekclock()
628 call get_map
629
630c compute global to local map (no processor info)
631 IEL=0
632 CALL IZERO(GLLEL,NELGT)
633 DO IEG=1,NELGT
634 IF (GLLNID(IEG).EQ.NID) THEN
635 IEL = IEL + 1
636 GLLEL(IEG)=IEL
637 NELT = IEL
638 if (ieg.le.nelgv) NELV = IEL
639 ENDIF
640c write(6,*) 'map2 ieg:',ieg,nelv,nelt,nelgv,nelgt
641 ENDDO
642
643c dist. global to local map to all processors
644 npass = 1 + nelgt/lelt
645 k=1
646 do ipass = 1,npass
647 m = nelgt - k + 1
648 m = min(m,lelt)
649 if (m.gt.0) call igop(gllel(k),iwork,'+ ',m)
650 k = k+m
651 enddo
652
653c compute local to global map
654c (i.e. returns global element number given local index and proc id)
655 do ieg=1,nelgt
656 mid =gllnid(ieg)
657 ie =gllel (ieg)
658 if (mid.eq.nid) lglel(ie)=ieg
659 enddo
660#endif
661
662 return
663 end
664c-----------------------------------------------------------------------
665
666#ifndef DPROCMAP
667
668c-----------------------------------------------------------------------
669 subroutine read_map(vertex,nlv,wk,mdw,ndw)
670
671 include 'SIZE'
672 include 'INPUT'
673 include 'PARALLEL'
674
675 integer vertex(nlv,1)
676 integer wk(mdw,ndw)
677
678 logical ifbswap,if_byte_swap_test
679
680 character*132 mapfle
681 character*1 mapfle1(132)
682 equivalence (mapfle,mapfle1)
683
684 character*132 hdr
685 character*5 version
686 real*4 test
687
688 logical ifma2,ifmap
689 integer e,eg,eg0,eg1
690 integer itmp20(20)
691
692 ierr = 0
693 ifma2 = .false.
694
695 if (nid.eq.0) then
696 lfname = ltrunc(reafle,132) - 4
697 call blank (mapfle,132)
698 call chcopy(mapfle,reafle,lfname)
699 call chcopy(mapfle1(lfname+1),'.map',4)
700 inquire(file=mapfle, exist=ifmap)
701
702 if (.not.ifmap) then
703 call chcopy(mapfle1(lfname+1),'.ma2',4)
704 inquire(file=mapfle, exist=ifma2)
705 endif
706
707 if(.not.ifmap .and. .not.ifma2) ierr = 1
708 endif
709 if(nid.eq.0) write(6,'(A,A)') ' Reading ', mapfle
710 call err_chk(ierr,' Cannot find map file!$')
711 call bcast(ifma2,lsize)
712 ierr = 0
713
714 if (nid.eq.0) then
715 if (ifma2) then
716 call byte_open(mapfle,ierr)
717 if(ierr.ne.0) goto 100
718
719 call blank(hdr,132)
720 call byte_read(hdr,132/4,ierr)
721 if(ierr.ne.0) goto 100
722
723 read (hdr,1) version,neli,nnzi
724 1 format(a5,2i12)
725
726 call byte_read(test,1,ierr)
727 if(ierr.ne.0) goto 100
728 ifbswap = if_byte_swap_test(test,ierr)
729 if(ierr.ne.0) goto 100
730 else
731 open(unit=80,file=mapfle,status='old',err=100)
732 read(80,*,err=100) neli,nnzi
733 endif
734 endif
735
736 call bcast(neli, ISIZE)
737
738 npass = 1 + (neli/ndw)
739 if (npass.gt.np) then
740 if (nid.eq.0) write(6,*) npass,np,neli,ndw,'Error get_vert_map'
741 call exitt
742 endif
743
744 len = 4*mdw*ndw
745 if (nid.gt.0.and.nid.lt.npass) msg_id=irecv(nid,wk,len)
746 call nekgsync
747
748 if (nid.eq.0) then
749 eg0 = 0
750 do ipass=1,npass
751 eg1 = min(eg0+ndw,neli)
752
753 if (ifma2) then
754 nwds = (eg1 - eg0)*(mdw-1)
755 call byte_read(wk,nwds,ierr)
756 if (ierr.ne.0) goto 200
757 if (ifbswap) call byte_reverse(wk,nwds,ierr)
758
759 m = eg1 - eg0
760 do eg=eg1,eg0+1,-1 ! reshuffle array
761 jj = (m-1)*(mdw-1) + 1
762 call icopy(itmp20,wk(jj,1),mdw-1)
763 call icopy(wk(1,m),itmp20 ,mdw-1)
764 m = m - 1
765 enddo
766 else
767 m = 0
768 do eg=eg0+1,eg1
769 m = m + 1
770 read(80,*,err=200) (wk(k,m),k=1,mdw-1)
771 enddo
772 endif
773
774 m = 0
775 do eg=eg0+1,eg1
776 m = m + 1
777 gllnid(eg) = wk(1,m) ! must still be divided
778 wk(mdw,m) = eg
779 enddo
780
781 if (ipass.lt.npass) call csend(ipass,wk,len,ipass,0) !send to ipass
782 eg0 = eg1
783 enddo
784
785 ntuple = m
786
787 if (ifma2) then
788 call byte_close(ierr)
789 else
790 close(80)
791 endif
792 elseif (nid.lt.npass) then
793 call msgwait(msg_id)
794 ntuple = ndw
795 else
796 ntuple = 0
797 endif
798
799 lng = isize*neli
800 call bcast(gllnid,lng)
801 call assign_gllnid(gllnid,gllel,nelgt,nelgv,np) ! gllel is used as scratch
802
803 nelt=0 ! Count number of elements on this processor
804 nelv=0
805 do eg=1,neli
806 if (gllnid(eg).eq.nid) then
807 if (eg.le.nelgv) nelv=nelv+1
808 if (eg.le.nelgt) nelt=nelt+1
809 endif
810 enddo
811
812c if (np.le.64) write(6,*) nid,nelv,nelt,nelgv,nelgt,' NELV'
813
814c NOW: crystal route vertex by processor id
815
816 ntuple_sum = iglsum(ntuple,1)
817 if (ntuple_sum .ne. nelgt) then
818 if (nid.eq.0) write(6,*) 'Error invalid tuple sum!'
819 call exitt
820 endif
821
822 do i=1,ntuple
823 eg=wk(mdw,i)
824 wk(1,i)=gllnid(eg) ! processor id for element eg
825 enddo
826
827 key = 1 ! processor id is in wk(1,:)
828 call fgslib_crystal_ituple_transfer(cr_h,wk,mdw,ntuple,ndw,key)
829
830 key = mdw ! Sort tuple list by eg
831 nkey = 1
832 call fgslib_crystal_ituple_sort(cr_h,wk,mdw,nelt,key,nkey)
833
834 iflag = 0
835 if (ntuple.ne.nelt) then
836 write(6,*) nid,ntuple,nelv,nelt,nelgt,' NELT FAIL'
837 write(6,*) 'Check that .map file and .rea file agree'
838 iflag=1
839 else
840 do e=1,nelt
841 call icopy(vertex(1,e),wk(2,e),nlv)
842 enddo
843 endif
844
845 iflag = iglmax(iflag,1)
846 if (iflag.gt.0) then
847 do mid=0,np-1
848 call nekgsync
849 if (mid.eq.nid)
850 $ write(6,*) nid,ntuple,nelv,nelt,nelgt,' NELT FB'
851 call nekgsync
852 enddo
853 call nekgsync
854 call exitt
855 endif
856
857 return
858
859 100 continue
860 call err_chk(ierr,'Error opening or reading map header$')
861
862 200 continue
863 call err_chk(ierr,'Error while reading map file$')
864
865 return
866 end
867c-----------------------------------------------------------------------
868 subroutine assign_gllnid(gllnid,iunsort,nelgt,nelgv,np)
869c
870 integer gllnid(1),iunsort(1),nelgt,np
871 integer e,eg
872
873
874 log2p = log2(np)
875 np2 = 2**log2p
876 if (np2.eq.np.and.nelgv.eq.nelgt) then ! std power of 2 case
877
878 npstar = ivlmax(gllnid,nelgt)+1
879 nnpstr = npstar/np
880 do eg=1,nelgt
881 gllnid(eg) = gllnid(eg)/nnpstr
882 enddo
883
884 return
885
886 elseif (np2.eq.np) then ! std power of 2 case, conjugate heat xfer
887
888c Assign fluid elements
889 npstar = max(np,ivlmax(gllnid,nelgv)+1)
890 nnpstr = npstar/np
891 do eg=1,nelgv
892 gllnid(eg) = gllnid(eg)/nnpstr
893 enddo
894
895c Assign solid elements
896 nelgs = nelgt-nelgv ! number of solid elements
897 npstar = max(np,ivlmax(gllnid(nelgv+1),nelgs)+1)
898 nnpstr = npstar/np
899 do eg=nelgv+1,nelgt
900 gllnid(eg) = gllnid(eg)/nnpstr
901 enddo
902
903 return
904
905 elseif (nelgv.ne.nelgt) then
906 call exitti
907 $ ('Conjugate heat transfer requires P=power of 2.$',np)
908 endif
909
910
911c Below is the code for P a non-power of two:
912
913c Split the sorted gllnid array (read from .map file)
914c into np contiguous partitions.
915
916c To load balance the partitions in case of mod(nelgt,np)>0
917c add 1 contiguous entry out of the sorted list to NODE_i
918c where i = np-mod(nelgt,np) ... np
919
920
921 nel = nelgt/np ! number of elements per processor
922 nmod = mod(nelgt,np) ! bounded between 1 ... np-1
923 npp = np - nmod ! how many paritions of size nel
924
925 ! sort gllnid
926 call isort(gllnid,iunsort,nelgt)
927
928 ! setup partitions of size nel
929 k = 0
930 do ip = 0,npp-1
931 do e = 1,nel
932 k = k + 1
933 gllnid(k) = ip
934 enddo
935 enddo
936 ! setup partitions of size nel+1
937 if(nmod.gt.0) then
938 do ip = npp,np-1
939 do e = 1,nel+1
940 k = k + 1
941 gllnid(k) = ip
942 enddo
943 enddo
944 endif
945
946 ! unddo sorting to restore initial ordering by
947 ! global element number
948 call iswapt_ip(gllnid,iunsort,nelgt)
949
950 return
951 end
952c-----------------------------------------------------------------------
953
954#endif
Note: See TracBrowser for help on using the repository browser.