KTH framework for Nek5000 toolboxes; testing version  0.0.1
connect2.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine readat
3 C
4 C Read in data from preprocessor input file (.rea)
5 C
6  include 'SIZE'
7  include 'INPUT'
8  include 'GEOM'
9  include 'PARALLEL'
10  include 'CTIMER'
11 
12  logical ifbswap,ifre2,parfound
13  character*132 string
14  integer idum(3*numsts+3)
15 
16  ierr = 0
17  call flush_io
18 
19  ! check if new rea file version exists
20  if(nid.eq.0) inquire(file=parfle, exist=parfound)
21  call bcast(parfound,lsize)
22  if (parfound) then
23  if(nio.eq.0) write(6,'(A,A)') ' Reading ', parfle
24  call readat_par
25  goto 99
26  endif
27 
28  etime0 = dnekclock_sync()
29 
30  if(nid.eq.0) then
31  write(6,'(A,A)') ' Reading ', reafle
32  open (unit=9,file=reafle,status='old', iostat=ierr)
33  endif
34 
35  call bcast(ierr,isize)
36  if (ierr .gt. 0) call exitti('Cannot open rea file!$',1)
37 
38 C Read parameters and logical flags
39  call rdparam
40  meshpartitioner=3 ! HYBRID (RSB+RCB)
41 
42 C Read Mesh Info
43  if(nid.eq.0) then
44  read(9,*) ! xfac,yfac,xzero,yzero
45  read(9,*) ! dummy
46  read(9,*) nelgs,ldimr,nelgv
47  nelgt = abs(nelgs)
48  endif
49  call bcast(ldimr,isize)
50  call bcast(nelgs,isize)
51  call bcast(nelgv,isize)
52  call bcast(nelgt,isize)
53  ifre2 = .false.
54  if (nelgs.lt.0) ifre2 = .true.
55 
56  call usrdat0
57 
58  if (nelgt.gt.350000 .and. .not.ifre2)
59  $ call exitti('Problem size requires .re2!$',1)
60 
61  if (ifre2) call read_re2_hdr(ifbswap, .true.) ! rank0 will open and read
62  call chk_nel ! make certain sufficient array sizes
63 
64  call mapelpr
65 
66  if (ifre2) then
67  call read_re2_data(ifbswap)
68  else
69  maxrd = 32 ! max # procs to read at once
70  mread = (np-1)/maxrd+1 ! mod param
71  iread = 0 ! mod param
72  x = 0
73  do i=0,np-1,maxrd
74  call nekgsync()
75  if (mod(nid,mread).eq.iread) then
76  if (nid.ne.0) then
77  open(unit=9,file=reafle,status='OLD')
78  call cscan(string,'MESH DATA',9)
79  read(9,*) string
80  endif
81  call rdmesh
82  call rdcurve ! Curved side data
83  call rdbdry ! Boundary Conditions
84  if (nid.ne.0) close(unit=9)
85  endif
86  iread = iread + 1
87  enddo
88  endif
89 
90 C Read Restart options / Initial Conditions / Drive Force
91  CALL rdicdf
92 C Read materials property data
93  CALL rdmatp
94 C Read history data
95  CALL rdhist
96 C Read output specs
97  CALL rdout
98 C Read objects
99  CALL rdobj
100 
101  call nekgsync()
102 
103 C End of input data, close read file.
104  if(nid.eq.0) then
105  close(unit=9)
106  call echopar
107  write(6,'(A,g13.5,A,/)') ' done :: read .rea file ',
108  $ dnekclock()-etime0,' sec'
109  endif
110 
111  99 call izero(boundaryid, size(boundaryid))
112  call izero(boundaryidt, size(boundaryidt))
113 
114  ifld = 2
115  if(ifflow) ifld = 1
116  do iel = 1,nelv
117  do ifc = 1,2*ndim
118  boundaryid(ifc,iel) = bc(5,ifc,iel,ifld)
119  enddo
120  enddo
121 
122  ntmsh = 0
123  do i=1,ldimt
124  if(iftmsh(1+i)) ntmsh = ntmsh + 1
125  enddo
126 
127  if (ntmsh.gt.0) then
128  do iel = 1,nelt
129  do ifc = 1,2*ndim
130  boundaryidt(ifc,iel) = bc(5,ifc,iel,2)
131  enddo
132  enddo
133  endif
134 
135  return
136  END
137 c-----------------------------------------------------------------------
138  subroutine vrdsmsh
139 C
140 C=====================================================================
141 C Verify that mesh and dssum are properly defined by performing
142 C a direct stiffness operation on the X,Y and Z coordinates.
143 C Note that periodic faces are not checked here.
144 C=====================================================================
145 C
146  include 'SIZE'
147  include 'TOTAL'
148  COMMON /scrns/ ta(lx1,ly1,lz1,lelt),tb(lx1,ly1,lz1,lelt)
149  $ ,qmask(lx1,ly1,lz1,lelt),tmp(2)
150  CHARACTER*3 CB
151 
152 c call vrdsmshx ! verify mesh topology
153 
154  if(nio.eq.0) write(*,*) 'verify mesh topology'
155 
156  ierr = 0
157  eps = 1.0e-04
158  eps = 1.0e-03
159  ifield = 1
160  if (nelgv.ne.nelgt .or. .not.ifflow) ifield = 2
161  nxyz1 = lx1*ly1*lz1
162  ntot = lx1*ly1*lz1*nelt
163  nfaces = 2*ldim
164 
165  xmx = glmax(xm1,ntot)
166  xmn = glmin(xm1,ntot)
167  ymx = glmax(ym1,ntot)
168  ymn = glmin(ym1,ntot)
169  zmx = glmax(zm1,ntot)
170  zmn = glmin(zm1,ntot)
171  if (nio.eq.0) write(6,*) xmn,xmx,' Xrange'
172  if (nio.eq.0) write(6,*) ymn,ymx,' Yrange'
173  if (nio.eq.0) write(6,*) zmn,zmx,' Zrange'
174 c return
175 C
176 C First check - use 1/Multiplicity
177 C
178  IF (ifheat) THEN
179  CALL copy(ta,tmult,ntot)
180  ELSE
181  CALL copy(ta,vmult,ntot)
182  ENDIF
183 c
184 c write(6,1)
185 c $(nid,'tab4',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
186 c 1 format(i3,a4,i3,16f5.2)
187 c
188  CALL dssum(ta,lx1,ly1,lz1)
189 c
190 c write(6,1)
191 c $(nid,'taaf',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
192 c
193  CALL rone (tb,ntot)
194  CALL sub2 (tb,ta,ntot)
195  DO 1000 ie=1,nelt
196  ieg=lglel(ie)
197  DO 1000 iz=1,lz1
198  DO 1000 iy=1,ly1
199  DO 1000 ix=1,lx1
200  IF (abs(tb(ix,iy,iz,ie)).GT.eps ) THEN
201  WRITE(6,1005) ix,iy,iz,ieg
202  $ ,xm1(ix,iy,iz,ie),ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
203  $ ,ta(ix,iy,iz,ie),eps
204 c WRITE(7,1005) IX,IY,IZ,IEG
205 c $ ,XM1(IX,IY,IZ,IE),TB(IX,IY,IZ,IE),TA(IX,IY,IZ,IE)
206 c $ ,QMASK(IX,IY,IZ,IE)
207  1005 FORMAT(2x,'WARNING: DSSUM problem at:',/
208  $ ,1x,'I,J,K,IE:',3i5,i12,/
209  $ ,2x,'Near X =',3g16.8,', d:',2g16.8)
210  ierr=4
211  ENDIF
212  1000 CONTINUE
213 C
214 C Set up QMASK quickly to annihilate checks on periodic bc's
215 C
216  CALL rone(qmask,ntot)
217  DO 100 iel=1,nelt
218  DO 100 iface=1,nfaces
219  cb =cbc(iface,iel,ifield)
220  IF (cb.EQ.'P '.or.cb.eq.'p ')
221  $ CALL facev(qmask,iel,iface,0.0,lx1,ly1,lz1)
222  100 CONTINUE
223  CALL dsop(qmask,'MUL',lx1,ly1,lz1)
224 
225 c xxmin = glmin(xm1,ntot)
226 c yymin = glmin(ym1,ntot)
227 c zzmin = glmin(zm1,ntot)
228 c xxmax = glmax(xm1,ntot)
229 c yymax = glmax(ym1,ntot)
230 c zzmax = glmax(zm1,ntot)
231 c if (nio.eq.0) write(6,7) xxmin,yymin,zzmin,xxmax,yymax,zzmax
232 c 7 format('xyz minmx2:',6g13.5)
233 
234 
235 
236 C
237 C X-component
238 C
239  CALL copy(ta,xm1,ntot)
240  CALL copy(tb,xm1,ntot)
241  CALL dsop(ta,'MIN',lx1,ly1,lz1)
242  CALL dsop(tb,'MAX',lx1,ly1,lz1)
243  CALL sub2(ta,xm1,ntot)
244  CALL sub2(tb,xm1,ntot)
245  CALL col2(ta,qmask,ntot)
246  CALL col2(tb,qmask,ntot)
247  DO 1100 ie=1,nelt
248  xscmax = vlmax(xm1(1,1,1,ie),nxyz1)
249  xscmin = vlmin(xm1(1,1,1,ie),nxyz1)
250  scal1=abs(xscmax-xscmin)
251  scal2=abs(xscmax)
252  scal3=abs(xscmin)
253  scal1=max(scal1,scal2)
254  scal1=max(scal1,scal3)
255  xscale = 1./scal1
256  ieg=lglel(ie)
257  DO 1100 iz=1,lz1
258  DO 1100 iy=1,ly1
259  DO 1100 ix=1,lx1
260  if (abs(ta(ix,iy,iz,ie)*xscale).gt.eps .or.
261  $ abs(tb(ix,iy,iz,ie)*xscale).gt.eps ) then
262  write(6,1105) ix,iy,iz,ieg
263  $ ,xm1(ix,iy,iz,ie),ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
264  $ ,tb(ix,iy,iz,ie),ta(ix,iy,iz,ie),xscale
265  1105 format(1x,'WARNING1 Element mesh mismatch at:',/
266  $ ,1x,'i,j,k,ie:',3i5,i12,/
267  $ ,1x,'Near X =',3g16.8,', d:',3g16.8)
268  ierr=1
269  endif
270  1100 CONTINUE
271 C
272 C Y-component
273 C
274  CALL copy(ta,ym1,ntot)
275  CALL copy(tb,ym1,ntot)
276  CALL dsop(ta,'MIN',lx1,ly1,lz1)
277  CALL dsop(tb,'MAX',lx1,ly1,lz1)
278  CALL sub2(ta,ym1,ntot)
279  CALL sub2(tb,ym1,ntot)
280  CALL col2(ta,qmask,ntot)
281  CALL col2(tb,qmask,ntot)
282  DO 1200 ie=1,nelt
283  yscmax = vlmax(ym1(1,1,1,ie),nxyz1)
284  yscmin = vlmin(ym1(1,1,1,ie),nxyz1)
285  scal1=abs(yscmax-yscmin)
286  scal2=abs(yscmax)
287  scal3=abs(yscmin)
288  scal1=max(scal1,scal2)
289  scal1=max(scal1,scal3)
290  yscale = 1./scal1
291  ieg=lglel(ie)
292  DO 1200 iz=1,lz1
293  DO 1200 iy=1,ly1
294  DO 1200 ix=1,lx1
295  IF (abs(ta(ix,iy,iz,ie)*yscale).GT.eps .OR.
296  $ abs(tb(ix,iy,iz,ie)*yscale).GT.eps ) THEN
297  WRITE(6,1205) ix,iy,iz,ieg
298  $ ,xm1(ix,iy,iz,ie),ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
299  $ ,tb(ix,iy,iz,ie),ta(ix,iy,iz,ie),yscale
300  1205 FORMAT(1x,'WARNING2 Element mesh mismatch at:',/
301  $ ,1x,'I,J,K,IE:',3i5,i12,/
302  $ ,1x,'Near Y =',3g16.8,', d:',3g16.8)
303  ierr=2
304  ENDIF
305  1200 CONTINUE
306 C
307 C Z-component
308 C
309  IF (if3d) THEN
310  CALL copy(ta,zm1,ntot)
311  CALL copy(tb,zm1,ntot)
312  CALL dsop(ta,'MIN',lx1,ly1,lz1)
313  CALL dsop(tb,'MAX',lx1,ly1,lz1)
314  CALL sub2(ta,zm1,ntot)
315  CALL sub2(tb,zm1,ntot)
316  CALL col2(ta,qmask,ntot)
317  CALL col2(tb,qmask,ntot)
318  DO 1300 ie=1,nelt
319  zscmax = vlmax(zm1(1,1,1,ie),nxyz1)
320  zscmin = vlmin(zm1(1,1,1,ie),nxyz1)
321  scal1=abs(zscmax-zscmin)
322  scal2=abs(zscmax)
323  scal3=abs(zscmin)
324  scal1=max(scal1,scal2)
325  scal1=max(scal1,scal3)
326  zscale = 1./scal1
327  ieg=lglel(ie)
328  DO 1300 iz=1,lz1
329  DO 1300 iy=1,ly1
330  DO 1300 ix=1,lx1
331  IF (abs(ta(ix,iy,iz,ie)*zscale).GT.eps .OR.
332  $ abs(tb(ix,iy,iz,ie)*zscale).GT.eps ) THEN
333  WRITE(6,1305) ix,iy,iz,ieg
334  $ ,xm1(ix,iy,iz,ie),ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
335  $ ,tb(ix,iy,iz,ie),ta(ix,iy,iz,ie),zscale
336  1305 FORMAT(1x,'WARNING3 Element mesh mismatch at:',/
337  $ ,1x,'I,J,K,IE:',3i5,i12,/
338  $ ,1x,'Near Z =',3g16.8,', d:',3g16.8)
339  ierr=3
340  ENDIF
341  1300 CONTINUE
342  ENDIF
343 
344  ierr = iglsum(ierr,1)
345  IF (ierr.gt.0) THEN
346  if(nid.eq.0) WRITE(6,1400)
347  1400 FORMAT
348  $ (' Mesh consistency check failed. EXITING in VRDSMSH.')
349  call exitt
350  ENDIF
351 
352  tmp(1)=ierr
353  CALL gop(tmp,tmp(2),'M ',1)
354  IF (tmp(1).ge.4.0) THEN
355  WRITE(6,1400)
356  $ (' Mesh consistency check failed. EXITING in VRDSMSH.')
357  call exitt
358  ENDIF
359 
360  if(nio.eq.0) then
361  write(6,*) 'done :: verify mesh topology'
362  write(6,*) ' '
363  endif
364 
365  return
366  end
367 c-----------------------------------------------------------------------
368  subroutine vrdsmshx ! verify mesh topology
369 C
370 C=====================================================================
371 C Verify that mesh and dssum are properly defined by performing
372 C a direct stiffness operation on the X,Y and Z coordinates.
373 C Note that periodic faces are not checked here.
374 C=====================================================================
375 C
376  include 'SIZE'
377  include 'TOTAL'
378  common /scrns/ tc(lx1,ly1,lz1,lelt),td(lx1,ly1,lz1,lelt)
379  $ , ta(lx1,ly1,lz1,lelt),tb(lx1,ly1,lz1,lelt)
380  $ , qmask(lx1,ly1,lz1,lelt)
381  CHARACTER*3 CB
382 C
383  ierr = 0
384  eps = 1.0e-04
385  eps = 1.0e-03
386  ifield = 1
387  IF (ifheat) ifield = 2
388  nxyz1 = lx1*ly1*lz1
389  ntot = lx1*ly1*lz1*nelt
390  nfaces = 2*ldim
391 
392  xmx = glmax(xm1,ntot)
393  xmn = glmin(xm1,ntot)
394  ymx = glmax(ym1,ntot)
395  ymn = glmin(ym1,ntot)
396  zmx = glmax(zm1,ntot)
397  zmn = glmin(zm1,ntot)
398  if (nio.eq.0) write(6,*) xmn,xmx,' Xrange'
399  if (nio.eq.0) write(6,*) ymn,ymx,' Yrange'
400  if (nio.eq.0) write(6,*) zmn,zmx,' Zrange'
401 c return
402 C
403 C First check - use 1/Multiplicity
404 C
405  IF (ifheat) THEN
406  CALL copy(ta,tmult,ntot)
407  ELSE
408  CALL copy(ta,vmult,ntot)
409  ENDIF
410 c
411 c write(6,1)
412 c $(nid,'tab4',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
413 c 1 format(i3,a4,i3,16f5.2)
414 c
415  CALL dssum(ta,lx1,ly1,lz1)
416 c
417 c write(6,1)
418 c $(nid,'taaf',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
419 c
420  CALL rone (tb,ntot)
421  CALL sub2 (tb,ta,ntot)
422  DO 1000 ie=1,nelt
423  ieg=lglel(ie)
424  DO 1000 iz=1,lz1
425  DO 1000 iy=1,ly1
426  DO 1000 ix=1,lx1
427  IF (abs(tb(ix,iy,iz,ie)).GT.eps ) THEN
428  WRITE(6,1005) ix,iy,iz,ieg
429  $ ,xm1(ix,iy,iz,ie),ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
430  $ ,ta(ix,iy,iz,ie),eps
431 c WRITE(7,1005) IX,IY,IZ,IEG
432 c $ ,XM1(IX,IY,IZ,IE),TB(IX,IY,IZ,IE),TA(IX,IY,IZ,IE)
433 c $ ,QMASK(IX,IY,IZ,IE)
434  1005 FORMAT(2x,'WARNING: DSSUM problem at:',/
435  $ ,1x,'I,J,K,IE:',3i5,i12,/
436  $ ,2x,'Near X =',3g16.8,', d:',2g16.8)
437  ierr=4
438  ENDIF
439  1000 CONTINUE
440 C
441 C Set up QMASK quickly to annihilate checks on periodic bc's
442 C
443  CALL rone(qmask,ntot)
444  DO 100 iel=1,nelt
445  DO 100 iface=1,nfaces
446  cb =cbc(iface,iel,ifield)
447  IF (cb.EQ.'P '.or.cb.eq.'p ')
448  $ CALL facev(qmask,iel,iface,0.0,lx1,ly1,lz1)
449  100 CONTINUE
450  call dsop(qmask,'MUL',lx1,ly1,lz1)
451 
452  xxmin = glmin(xm1,ntot)
453  yymin = glmin(ym1,ntot)
454  zzmin = glmin(zm1,ntot)
455  xxmax = glmax(xm1,ntot)
456  yymax = glmax(ym1,ntot)
457  zzmax = glmax(zm1,ntot)
458  if (nio.eq.0) write(6,7) xxmin,yymin,zzmin,xxmax,yymax,zzmax
459  7 format('xyz minmx2:',6g13.5)
460 
461 
462 C
463 C X-component
464 C
465  call copy(ta,xm1,ntot)
466  call copy(tb,xm1,ntot)
467  call dsop(ta,'min',lx1,ly1,lz1)
468  call dsop(tb,'max',lx1,ly1,lz1)
469 
470  call copy(tc,xm1,ntot)
471  call copy(td,xm1,ntot)
472  call dsop(tc,'min',lx1,ly1,lz1)
473  call dsop(td,'max',lx1,ly1,lz1)
474 
475  xxmin = glmin(xm1,ntot)
476  xxmax = glmax(xm1,ntot)
477  yymax = glmax(ta ,ntot)
478  yymin = glmin(ta ,ntot)
479  zzmin = glmin(tb ,ntot)
480  zzmax = glmax(tb ,ntot)
481  if (nio.eq.0) write(6,9) xxmin,yymin,zzmin,xxmax,yymax,zzmax
482  9 format('xyz minmx3:',6g13.5)
483 
484  CALL sub2(ta,xm1,ntot)
485  CALL sub2(tb,xm1,ntot)
486 
487  xxmin = glmin(qmask,ntot)
488  xxmax = glmax(qmask,ntot)
489  yymax = glmax(ta ,ntot)
490  yymin = glmin(ta ,ntot)
491  zzmin = glmin(tb ,ntot)
492  zzmax = glmax(tb ,ntot)
493  if (nio.eq.0) write(6,19) xxmin,yymin,zzmin,xxmax,yymax,zzmax
494  19 format('xyz minmx4:',6g13.5)
495 
496  CALL col2(ta,qmask,ntot)
497  CALL col2(tb,qmask,ntot)
498 
499  xxmin = glmin(qmask,ntot)
500  xxmax = glmax(qmask,ntot)
501  yymax = glmax(ta ,ntot)
502  yymin = glmin(ta ,ntot)
503  zzmin = glmin(tb ,ntot)
504  zzmax = glmax(tb ,ntot)
505  if (nio.eq.0) write(6,29) xxmin,yymin,zzmin,xxmax,yymax,zzmax
506  29 format('xyz minmx5:',6g13.5)
507 
508  DO 1100 ie=1,nelt
509  xscmax = vlmax(xm1(1,1,1,ie),nxyz1)
510  xscmin = vlmin(xm1(1,1,1,ie),nxyz1)
511  scal1=abs(xscmax-xscmin)
512  scal2=abs(xscmax)
513  scal3=abs(xscmin)
514  scal1=max(scal1,scal2)
515  scal1=max(scal1,scal3)
516  xscale = 1./scal1
517  ieg=lglel(ie)
518  DO 1100 iz=1,lz1
519  DO 1100 iy=1,ly1
520  DO 1100 ix=1,lx1
521  if (abs(ta(ix,iy,iz,ie)*xscale).gt.eps .or.
522  $ abs(tb(ix,iy,iz,ie)*xscale).gt.eps ) then
523  write(6,1105) nid,ix,iy,iz,ie,ieg
524  $ ,xm1(ix,iy,iz,ie),tc(ix,iy,iz,ie),td(ix,iy,iz,ie)
525  $ ,ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
526  $ ,ta(ix,iy,iz,ie),tb(ix,iy,iz,ie),xscale
527  $ ,qmask(ix,iy,iz,ie)
528  1105 format(i4.4,1x,'ie:',3i3,i10,i10,1p9e11.3)
529 c1105 format(i4.4,1x,'ie:',3i3,i6,1p9e11.3)
530  ierr=1
531  goto 1101
532  endif
533  1100 CONTINUE
534  1101 CONTINUE
535 
536  xxmin = glmin(xm1,ntot)
537  xxmax = glmax(xm1,ntot)
538  yymax = glmax(ta ,ntot)
539  yymin = glmin(ta ,ntot)
540  zzmin = glmin(tb ,ntot)
541  zzmax = glmax(tb ,ntot)
542  if (nio.eq.0) write(6,39) xxmin,yymin,zzmin,xxmax,yymax,zzmax
543  39 format('xyz minmx5:',6g13.5)
544 
545 c ifvo = .true.
546 c ifpo = .false.
547 c ifto = .true.
548 c call outpost(xm1,ta,tb,pr,qmask,' ')
549 c call exitt
550 
551  return
552  end
553 c-----------------------------------------------------------------------
554  subroutine rotat2(xyz,angle,npts)
555 C
556 C Rotate NPTS through ANGLE (in two directions IF3D).
557 C
558  include 'SIZE'
559  include 'INPUT'
560  dimension xyz(3,1)
561  COMMON /ctmp0/ rmtrx(3,3),rx(3,3),rz(3,3),xyzn(3,10)
562 C
563  sina=sin(angle)
564  cosa=cos(angle)
565  CALL rzero(rx,9)
566  CALL rzero(rz,9)
567  rx(1,1)=cosa
568  rx(2,2)=cosa
569  rx(1,2)=sina
570  rx(2,1)=-sina
571  rx(3,3)=1.0
572  IF (if3d) THEN
573  rz(1,1)=cosa
574  rz(3,3)=cosa
575  rz(1,3)=sina
576  rz(3,1)=-sina
577  rz(2,2)=1.0
578  ELSE
579  rz(1,1)=1.0
580  rz(2,2)=1.0
581  rz(3,3)=1.0
582  ENDIF
583  CALL mxm(rx,3,rz,3,rmtrx,3)
584 C
585 C Strip mine mxms in chunks of 10:
586  DO 100 i=1,npts-10,10
587  CALL mxm(rmtrx,3,xyz(1,i),3,xyzn,10)
588  CALL copy(xyz(1,i),xyzn,30)
589  100 CONTINUE
590  n10=mod1(npts,10)
591  i=npts-n10+1
592  CALL rzero(xyzn,30)
593  IF (n10.GT.0) THEN
594  CALL mxm(rmtrx,3,xyz(1,i),3,xyzn,n10)
595  CALL copy(xyz(1,i),xyzn,3*n10)
596  ENDIF
597 C
598  return
599  END
600 c-----------------------------------------------------------------------
601  subroutine scale(xyzl,nl)
602 C
603 C Rescale XYZL such that the mean value of IXX=IYY=IZZ for each element.
604 C
605  include 'SIZE'
606  include 'INPUT'
607  dimension xyzl(3,8,lelt)
608  COMMON /ctmp0/ vo(lelt),xyzi(3,lelt),cg(3,lelt)
609  $ ,ti(6),work(6)
610 C
611 C Compute volumes -
612 C
613  CALL volume2(vo,xyzl,nl)
614  vtot=glsum(vo,nl)
615 C
616 C Compute (weighted) average inertia for each element.
617 C
618  ncrnr=2**ldim
619  CALL rzero(ti,6)
620  DO 100 il=1,nl
621  vo0 = vo(il)/vtot
622  CALL inrtia(xyzi(1,il),cg(1,il),xyzl(1,1,il),ncrnr,1)
623  ti(1)=ti(1)+xyzi(1,il)*vo0
624  ti(2)=ti(2)+xyzi(2,il)*vo0
625  ti(3)=ti(3)+xyzi(3,il)*vo0
626  ti(4)=ti(4)+cg(1,il) *vo0
627  ti(5)=ti(5)+cg(2,il) *vo0
628  ti(6)=ti(6)+cg(3,il) *vo0
629  100 CONTINUE
630  CALL gop(ti,work,'+ ',6)
631  xi =sqrt(ti(1))
632  yi =sqrt(ti(2))
633  zi =1.0
634  IF (if3d) zi=sqrt(ti(3))
635 C
636 C Rescale ( & shift to a nearly mean zero )
637 C
638  DO 200 il=1,nl
639  DO 200 ic=1,ncrnr
640  xyzl(1,ic,il)=(xyzl(1,ic,il)-ti(4))/xi
641  xyzl(2,ic,il)=(xyzl(2,ic,il)-ti(5))/yi
642  xyzl(3,ic,il)=(xyzl(3,ic,il)-ti(6))/zi
643  200 CONTINUE
644 C
645  return
646  END
647 c-----------------------------------------------------------------------
648  subroutine inrtia(xyzi,cg,xyzl,n,itype)
649 C
650 C Compute cg and inertia for a collection of unit point masses.
651 C This is a global (multiprocessor) operation, only IF itype=2.
652 C
653  dimension xyzi(3),cg(3),xyzl(3,1)
654  dimension ti(4),work(4)
655 C
656  ti(1)=0.0
657  ti(2)=0.0
658  ti(3)=0.0
659  ti(4)=n
660  DO 100 i=1,n
661  ti(1)=ti(1)+xyzl(1,i)
662  ti(2)=ti(2)+xyzl(2,i)
663  ti(3)=ti(3)+xyzl(3,i)
664  100 CONTINUE
665  IF (itype.EQ.2) CALL gop(ti,work,'+ ',4)
666  IF (ti(4).EQ.0.0) ti(4)=1.0
667  cg(1)=ti(1)/ti(4)
668  cg(2)=ti(2)/ti(4)
669  cg(3)=ti(3)/ti(4)
670 C
671  ti(1)=0.0
672  ti(2)=0.0
673  ti(3)=0.0
674  DO 200 i=1,n
675  ti(1)=ti(1)+( xyzl(1,i)-cg(1) )**2
676  ti(2)=ti(2)+( xyzl(2,i)-cg(2) )**2
677  ti(3)=ti(3)+( xyzl(3,i)-cg(3) )**2
678  200 CONTINUE
679  IF (itype.EQ.2) CALL gop(ti,work,'+ ',3)
680  ti(1)=ti(1)/ti(4)
681  ti(2)=ti(2)/ti(4)
682  ti(3)=ti(3)/ti(4)
683  IF (itype.EQ.2) THEN
684 C std. def'n of inertia.
685  xyzi(1)=ti(2)+ti(3)
686  xyzi(2)=ti(3)+ti(1)
687  xyzi(3)=ti(1)+ti(2)
688  ELSE
689  xyzi(1)=ti(1)
690  xyzi(2)=ti(2)
691  xyzi(3)=ti(3)
692  ENDIF
693 C
694  return
695  END
696 c-----------------------------------------------------------------------
697  subroutine volume2(vol,xyz,n)
698  include 'SIZE'
699  include 'INPUT'
700  dimension xyz(3,2,2,2,1)
701  dimension vol(1)
702 C
703  DO 1000 ie=1,n
704  vol(ie)=0.0
705  IF (if3d) THEN
706  DO 20 k=1,2
707  DO 20 j=1,2
708  DO 20 i=1,2
709  vol1 = (xyz(1,2,j,k,ie)-xyz(1,1,j,k,ie))
710  $ * (xyz(2,i,2,k,ie)-xyz(2,i,1,k,ie))
711  $ * (xyz(3,i,j,2,ie)-xyz(3,i,j,1,ie))
712  vol2 = (xyz(1,2,j,k,ie)-xyz(1,1,j,k,ie))
713  $ * (xyz(2,i,j,2,ie)-xyz(2,i,j,1,ie))
714  $ * (xyz(3,i,2,k,ie)-xyz(3,i,1,k,ie))
715  vol3 = (xyz(1,i,2,k,ie)-xyz(1,i,1,k,ie))
716  $ * (xyz(2,2,j,k,ie)-xyz(2,1,j,k,ie))
717  $ * (xyz(3,i,j,2,ie)-xyz(3,i,j,1,ie))
718  vol4 = (xyz(1,i,j,2,ie)-xyz(1,i,j,1,ie))
719  $ * (xyz(2,i,2,k,ie)-xyz(2,i,1,k,ie))
720  $ * (xyz(3,i,2,k,ie)-xyz(3,i,1,k,ie))
721  vol5 = (xyz(1,i,2,k,ie)-xyz(1,i,1,k,ie))
722  $ * (xyz(2,i,j,2,ie)-xyz(2,i,j,1,ie))
723  $ * (xyz(3,2,j,k,ie)-xyz(3,1,j,k,ie))
724  vol6 = (xyz(1,i,j,2,ie)-xyz(1,i,j,1,ie))
725  $ * (xyz(2,i,2,k,ie)-xyz(2,i,1,k,ie))
726  $ * (xyz(3,2,j,k,ie)-xyz(3,1,j,k,ie))
727  vol(ie) = vol(ie)+vol1+vol2+vol3+vol4+vol5+vol6
728  20 CONTINUE
729  vol(ie)=vol(ie)/8.0
730  ELSE
731 C 2-D:
732  DO 40 j=1,2
733  DO 40 i=1,2
734  vol1 = (xyz(1,2,j,1,ie)-xyz(1,1,j,1,ie))
735  $ * (xyz(2,i,2,1,ie)-xyz(2,i,1,1,ie))
736  vol3 = (xyz(1,i,2,1,ie)-xyz(1,i,1,1,ie))
737  $ * (xyz(2,2,j,1,ie)-xyz(2,1,j,1,ie))
738  vol(ie)=vol(ie)+vol1+vol3
739  40 CONTINUE
740  vol(ie)=vol(ie)/4.0
741  ENDIF
742  vol(ie)=abs(vol(ie))
743  1000 CONTINUE
744 C
745  return
746  END
747 c-----------------------------------------------------------------------
748  subroutine findcg(cg,xyz,n)
749 C
750 C Compute cg for N elements.
751 C
752  include 'SIZE'
753  dimension cg(3,1),xyz(3,8,1)
754 C
755  ncrnr=2**ldim
756  CALL rzero(cg,3*n)
757  DO 100 i =1,n
758  DO 100 ic=1,ncrnr
759  cg(1,i)=cg(1,i)+xyz(1,ic,i)
760  cg(2,i)=cg(2,i)+xyz(2,ic,i)
761  cg(3,i)=cg(3,i)+xyz(3,ic,i)
762  100 CONTINUE
763  tmp=1.0/(ncrnr)
764  CALL cmult(cg,tmp,3*n)
765  return
766  END
767 c-----------------------------------------------------------------------
768  subroutine divide(list1,list2,nl1,nl2,ifok,list,nl,xyzi,cg,WGT)
769 C
770 C Divide the elements associated with this subdomain according to
771 C the direction having the smallest moment of inertia (the "long"
772 C direction).
773 C
774  include 'SIZE'
775  include 'INPUT'
776  include 'PARALLEL'
777  include 'TSTEP'
778 C
779  dimension list(lelt),list1(lelt),list2(lelt)
780  dimension xyzi(3),cg(3,lelt),wgt(1)
781  COMMON /ctmp0/ xcg(lelt),ycg(lelt),zcg(lelt)
782  REAL IXX,IYY,IZZ
783  INTEGER WORK(2),WRK2(2)
784  LOGICAL IFOK
785 C
786 C Choose "long" direction:
787 C
788  ixx=xyzi(1)
789  iyy=xyzi(2)
790  izz=xyzi(3)
791  IF (if3d) THEN
792  IF (ixx.LE.iyy.AND.ixx.LE.izz) THEN
793  DO 104 ie=1,nl
794  xcg(ie)=cg(1,ie)
795  ycg(ie)=cg(2,ie)
796  zcg(ie)=cg(3,ie)
797  104 CONTINUE
798  ELSEIF (iyy.LE.ixx.AND.iyy.LE.izz) THEN
799  DO 106 ie=1,nl
800  xcg(ie)=cg(2,ie)
801  ycg(ie)=cg(3,ie)
802  zcg(ie)=cg(1,ie)
803  106 CONTINUE
804  ELSEIF (izz.LE.ixx.AND.izz.LE.iyy) THEN
805  DO 108 ie=1,nl
806  xcg(ie)=cg(3,ie)
807  ycg(ie)=cg(1,ie)
808  zcg(ie)=cg(2,ie)
809  108 CONTINUE
810  ENDIF
811  ELSE
812 C 2-D:
813  IF (ixx.LE.iyy) THEN
814  DO 114 ie=1,nl
815  xcg(ie)=cg(1,ie)
816  ycg(ie)=cg(2,ie)
817  114 CONTINUE
818  ELSE
819  DO 116 ie=1,nl
820  xcg(ie)=cg(2,ie)
821  ycg(ie)=cg(1,ie)
822  116 CONTINUE
823  ENDIF
824  ENDIF
825  call col2(xcg,wgt,nl)
826  call col2(ycg,wgt,nl)
827  call col2(zcg,wgt,nl)
828 C
829 C Find median value of CG to determine dividing point:
830 C
831  xm=fmdian(xcg,nl,ifok)
832  ym=fmdian(ycg,nl,ifok)
833  zm=0.0
834  IF (if3d) zm=fmdian(zcg,nl,ifok)
835 C
836 C Diagnostics
837 C
838  IF (.NOT.ifok) THEN
839  WRITE(6,130) nid,nl,xm,ym,zm
840  DO 120 il=1,nl
841  WRITE(6,135) nid,il,xcg(il),ycg(il),zcg(il)
842  120 CONTINUE
843  130 FORMAT(i10,'DIVIDE: NL,XM,YM,ZM',i3,3f12.5)
844  135 FORMAT(i10,'DIVIDE: NID,IL,XC,YC,ZCG',i4,3f12.5)
845  ENDIF
846 C
847 C=============================================================
848 C Divide LIST into LIST1 (XCG < XM) and LIST2 (XCG>XM).
849 C=============================================================
850 C
851  nl1=0
852  nl2=0
853  DO 200 ie=1,nl
854  IF (xcg(ie).LT.xm) THEN
855  nl1=nl1+1
856  list1(nl1)=list(ie)
857  ENDIF
858  IF (xcg(ie).GT.xm) THEN
859  nl2=nl2+1
860  list2(nl2)=list(ie)
861  ENDIF
862  IF (xcg(ie).EQ.xm) THEN
863 C
864 C We have to look at the other directions to arrive at
865 C a unique subdivision algortithm.
866 C
867 C
868 C More Diagnostics
869 C
870  IF (.NOT.ifok) WRITE(6,201) nid,ie,xcg(ie),xm
871  201 FORMAT(i10,'DIVIDE: IE,XCG,XM:',i4,3f12.5)
872 C
873  IF (ycg(ie).LT.ym) THEN
874  nl1=nl1+1
875  list1(nl1)=list(ie)
876  ENDIF
877  IF (ycg(ie).GT.ym) THEN
878  nl2=nl2+1
879  list2(nl2)=list(ie)
880  ENDIF
881  IF (ycg(ie).EQ.ym) THEN
882 C look at 3rd direction.
883  IF (if3d .AND. zcg(ie).LT.zm) THEN
884  nl1=nl1+1
885  list1(nl1)=list(ie)
886  ELSE IF (if3d .AND. zcg(ie).GT.zm) THEN
887  nl2=nl2+1
888  list2(nl2)=list(ie)
889  ELSE
890 C for 2- or 3-D intdeterminate case:
891  nl1=nl1+1
892  list1(nl1)=list(ie)
893  ENDIF
894  ENDIF
895 C
896  ENDIF
897  200 CONTINUE
898 C
899 C Check for an even distribution (i.e. - not different by
900 C more than 1):
901 C
902  ifok=.true.
903  work(1)=nl1
904  work(2)=nl2
905  CALL igop(work,wrk2,'+ ',2)
906  IF (abs(work(1)-work(2)).GT.1) ifok=.false.
907 C
908  return
909  END
910 c-----------------------------------------------------------------------
911  subroutine bufchk(buf,n)
912  integer n
913  real buf(n)
914  do i=1,n
915  write(6,*) buf(i), ' whhhh'
916  enddo
917  return
918  end
919 c-----------------------------------------------------------------------
920  subroutine chk_xyz
921  include 'SIZE'
922  include 'TOTAL'
923  integer e,f,eg
924 
925  do e=1,nelt
926  eg = lglel(e)
927  write(6,1) nid,eg,e,(cbc(f,e,1),f=1,6)
928  enddo
929  1 format(3i12,6(1x,a3),' cbc')
930 
931  return
932  end
933 c-----------------------------------------------------------------------
934  subroutine chk_nel
935  include 'SIZE'
936  include 'TOTAL'
937 
938  neltmx=np*lelt
939  nelvmx=np*lelv
940 
941  neltmx=min(neltmx,lelg)
942  nelvmx=min(nelvmx,lelg)
943 
944  nelgt = iglmax(nelgt,1)
945  nelgv = iglmax(nelgv,1)
946 
947 c write(6,*) nid,' inside chk_nel',nelgt,neltmx,nelvmx
948 
949  if (nelgt.gt.neltmx.or.nelgv.gt.nelvmx) then
950  if (nid.eq.0) then
951  lelt_needed = nelgt/np
952  if (mod(nelgt,np).ne.0) lelt_needed = lelt_needed + 1
953  write(6,12) lelt,lelg,lelt_needed,np,nelgt
954  12 format(//,2x,'ABORT: Problem size too large!'
955  $ ,/,2x
956  $ ,/,2x,'This solver has been compiled for:'
957  $ ,/,2x,' number of elements/proc (lelt):',i12
958  $ ,/,2x,' total number of elements (lelg):',i12
959  $ ,/,2x
960  $ ,/,2x,'Recompile with the following SIZE parameters:'
961  $ ,/,2x,' lelt >= ',i12,' for np = ',i12
962  $ ,/,2x,' lelg >= ',i12,/)
963 c write(6,*)'help:',lp,np,nelvmx,nelgv,neltmx,nelgt
964 c write(6,*)'help:',lelt,lelv,lelgv
965  endif
966  call exitt
967  endif
968 
969  if(nelgt.gt.nelgt_max) then
970  if(nid.eq.0) write(6,*)
971  $ 'ABORT: Total number of elements too large!',
972  $ ' nel_max = ', nelgt_max
973  call exitt
974  endif
975 
976  if (nelt.gt.lelt) then
977  write(6,'(A,3I12)') 'ABORT: nelt>lelt!', nid, nelt, lelt
978  call exitt
979  endif
980 
981  return
982  end
983 c-----------------------------------------------------------------------
984  subroutine cscan(sout,key,nk)
985 
986  character*132 sout,key
987  character*132 string
988  character*1 string1(132)
989  equivalence(string1,string)
990 c
991  do i=1,100000000
992  call blank(string,132)
993  read (nk,80,end=100,err=100) string
994  call chcopy(sout,string,132)
995 c write (6,*) string
996  if (indx1(string,key,nk).ne.0) return
997  enddo
998  100 continue
999 c
1000  80 format(a132)
1001  return
1002 
1003  end
1004 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
igop
subroutine igop(x, w, op, n)
Definition: comm_mpi.f:247
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
nekgsync
subroutine nekgsync()
Definition: comm_mpi.f:502
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
bcast
subroutine bcast(buf, len)
Definition: comm_mpi.f:431
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
dnekclock_sync
real *8 function dnekclock_sync()
Definition: comm_mpi.f:401
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
inrtia
subroutine inrtia(xyzi, cg, xyzl, n, itype)
Definition: connect2.f:649
chk_nel
subroutine chk_nel
Definition: connect2.f:935
chk_xyz
subroutine chk_xyz
Definition: connect2.f:921
bufchk
subroutine bufchk(buf, n)
Definition: connect2.f:912
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
vrdsmsh
subroutine vrdsmsh
Definition: connect2.f:139
volume2
subroutine volume2(vol, xyz, n)
Definition: connect2.f:698
readat
subroutine readat
Definition: connect2.f:3
rotat2
subroutine rotat2(xyz, angle, npts)
Definition: connect2.f:555
cscan
subroutine cscan(sout, key, nk)
Definition: connect2.f:985
divide
subroutine divide(list1, list2, nl1, nl2, ifok, list, nl, xyzi, cg, WGT)
Definition: connect2.f:769
vrdsmshx
subroutine vrdsmshx
Definition: connect2.f:369
findcg
subroutine findcg(cg, xyz, n)
Definition: connect2.f:749
echopar
subroutine echopar
Definition: drive2.f:291
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
indx1
integer function indx1(S1, S2, L2)
Definition: ic.f:1357
mapelpr
subroutine mapelpr()
Definition: map2.f:3
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
glmin
function glmin(a, n)
Definition: math.f:973
mod1
function mod1(i, n)
Definition: math.f:509
iglsum
function iglsum(a, n)
Definition: math.f:926
col2
subroutine col2(a, b, n)
Definition: math.f:564
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
blank
subroutine blank(A, N)
Definition: math.f:19
iglmax
function iglmax(a, n)
Definition: math.f:913
izero
subroutine izero(a, n)
Definition: math.f:215
fmdian
function fmdian(a, n, ifok)
Definition: math.f:1004
glsum
function glsum(x, n)
Definition: math.f:861
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
rone
subroutine rone(a, n)
Definition: math.f:230
vlmin
real function vlmin(vec, n)
Definition: math.f:357
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
readat_par
subroutine readat_par
Definition: reader_par.f:3
read_re2_hdr
subroutine read_re2_hdr(ifbswap, ifverbose)
Definition: reader_re2.f:772
read_re2_data
subroutine read_re2_data(ifbswap)
Definition: reader_re2.f:3
rdout
subroutine rdout
Definition: reader_rea.f:1030
rdmesh
subroutine rdmesh
Definition: reader_rea.f:519
rdparam
subroutine rdparam
Definition: reader_rea.f:3
rdcurve
subroutine rdcurve
Definition: reader_rea.f:609
rdbdry
subroutine rdbdry
Definition: reader_rea.f:702
rdicdf
subroutine rdicdf
Definition: reader_rea.f:881
rdobj
subroutine rdobj
Definition: reader_rea.f:1114
rdhist
subroutine rdhist
Definition: reader_rea.f:997
rdmatp
subroutine rdmatp
Definition: reader_rea.f:948
flush_io
subroutine flush_io
Definition: subs1.f:1561
qmask
subroutine qmask(R1, R2, R3, R1MASK, R2MASK, R3MASK, NEL)
Definition: subs2.f:1836