source: CIVL/examples/fortran/nek5000/core/ic.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: 80.8 KB
Line 
1c-----------------------------------------------------------------------
2 subroutine setics
3C-----------------------------------------------------------------------
4C
5C Set initial conditions.
6C
7C-----------------------------------------------------------------------
8 INCLUDE 'SIZE'
9 INCLUDE 'DEALIAS'
10 INCLUDE 'INPUT'
11 INCLUDE 'IXYZ'
12 INCLUDE 'GEOM'
13 INCLUDE 'SOLN'
14 INCLUDE 'MASS'
15 INCLUDE 'MVGEOM'
16 INCLUDE 'PARALLEL'
17 INCLUDE 'TSTEP'
18
19 logical iffort( ldimt1,0:lpert)
20 $ , ifrest(0:ldimt1,0:lpert)
21 $ , ifprsl( ldimt1,0:lpert)
22
23 LOGICAL IFANYP
24 common /inelr/ nelrr
25 common /ctmp1/ work(lx1,ly1,lz1,lelv)
26 $ , ta1 (lx2,ly1,lz1)
27 $ , ta2 (lx2,ly2,lz1)
28 integer*8 ntotg,nn
29
30 real psmax(ldimt)
31
32 if(nio.eq.0) write(6,*) 'set initial conditions'
33
34
35 nxyz2=lx2*ly2*lz2 ! Initialize all fields:
36 ntot2=nxyz2*nelv
37 nxyz1=lx1*ly1*lz1
38 ntott=nelt*nxyz1
39 ntotv=nelv*nxyz1
40 ltott=lelt*nxyz1
41
42 call rzero(vx,ntott)
43 call rzero(vy,ntott)
44 call rzero(vz,ntott)
45 call rzero(pr,nxyz2*nelt)
46 do 10 ifld=1,ldimt
47 call rzero(t(1,1,1,1,ifld),ntott)
48 10 continue
49
50 jp = 0 ! Set counter for perturbation analysis
51
52 irst = param(46) ! for lee's restart (rarely used)
53 if (irst.gt.0) call setup_convect(2)
54
55c If moving geometry then add a perturbation to the
56c mesh coordinates (see Subroutine INIGEOM)
57
58 if (ifmvbd) call ptbgeom
59
60C Find out what type of i.c. is requested
61C Current options:
62C
63C (1) - User specified fortran function (default is zero i.c.)
64C (2) - Restart from file(s)
65C (3) - Activate pre-solver => steady diffusion / steady Stokes
66C
67C If option (2) is requested, also return with the name of the
68C restart file(s) together with the associated dump number
69
70 call slogic (iffort,ifrest,ifprsl,nfiles)
71
72C ***** TEMPERATURE AND PASSIVE SCALARS ******
73C
74C Check if any pre-solv necessary for temperature/passive scalars
75
76 IFANYP = .FALSE.
77 DO 100 IFIELD=2,NFIELD
78 IF (IFPRSL(IFIELD,jp)) THEN
79 IF (NIO.EQ.0) WRITE(6,101) IFIELD
80 IFANYP = .TRUE.
81 ENDIF
82 100 CONTINUE
83 101 FORMAT(2X,'Using PRESOLVE option for field',I2,'.')
84
85C Fortran function initial conditions for temp/pass. scalars.
86 maxfld = nfield
87 if (ifmhd) maxfld = npscal+3
88
89c Always call nekuic (pff, 12/7/11)
90 do ifield=1,maxfld
91 if (nio.eq.0) write(6,*) 'nekuic (1) for ifld ', ifield
92 call nekuic
93 enddo
94
95C If any pre-solv, do pre-solv for all temperatur/passive scalar fields
96 if (ifanyp) call prsolvt
97
98 jp = 0 ! jp=0 --> base field, not perturbation field
99 do 200 ifield=2,maxfld
100 if (iffort(ifield,jp)) then
101 if (nio.eq.0) write(6,*) 'call nekuic for ifld ', ifield
102 call nekuic
103 endif
104 200 continue
105
106 if (ifpert) then
107 ifield=2
108 do jp=1,npert
109 if (nio.eq.0) write(6,*) 'nekuicP',ifield,jp,iffort(ifield,jp)
110 if (iffort(ifield,jp)) call nekuic
111 enddo
112 endif
113 jp = 0
114
115
116 call nekgsync()
117 call restart(nfiles) ! Check restart files
118 call nekgsync()
119
120
121C ***** VELOCITY ******
122C (If restarting for V, we're done,
123C ...else, do pre-solv for fluid if requested.)
124
125 ifield = 1
126 if (ifprsl(ifield,jp)) call prsolvv
127
128
129C Fortran function initial conditions for velocity.
130 ifield = 1
131 if (iffort(ifield,jp)) then
132 if (nio.eq.0) write(6,*) 'call nekuic for vel '
133 call nekuic
134 endif
135c
136 if (ifpert) then
137 ifield=1
138 do jp=1,npert
139 if (iffort(ifield,jp)) call nekuic
140 if (nio.eq.0) write(6,*) 'ic vel pert:',iffort(1,jp),jp
141 enddo
142 endif
143 jp = 0
144
145 ntotv = lx1*ly1*lz1*nelv
146
147C Initial mesh velocities
148 if (ifmvbd) call opcopy (wx,wy,wz,vx,vy,vz)
149 if (ifmvbd.and..not.ifrest(0,jp)) call meshv (2)
150
151C If convection-diffusion of a passive scalar with a fixed velocity field,
152C make sure to fill up lagged arrays since this will not be done in
153C the time-stepping procedure (no flow calculation) (01/18/91 -EMR).
154
155 if (.not.ifflow.and.ifheat) then
156 ITEST=0
157 DO 400 IFIELD=2,NFIELD
158 IF (IFADVC(IFIELD)) ITEST=1
159 400 CONTINUE
160 IF (ITEST.EQ.1) THEN
161 NBDMAX = 3
162 NBDSAV = NBDINP
163 NBDINP = NBDMAX
164 DO 500 I=1,NBDMAX
165 CALL LAGVEL
166 500 CONTINUE
167 NBDINP = NBDSAV
168 ENDIF
169 ENDIF
170
171C Ensure that all processors have the same time as node 0.
172 if (nid.ne.0) time=0.0
173 time=glsum(time,1)
174
175 nxyz1=lx1*ly1*lz1
176 ntott=nelt*nxyz1
177 ntotv=nelv*nxyz1
178 nn = nxyz1
179 ntotg=nelgv*nn
180
181 if (.not.ifdg) then
182 ifield = 2
183 if (ifflow) ifield = 1
184 call rone(work,ntotv)
185 ifield = 1
186 call dssum (work,lx1,ly1,lz1)
187 call col2 (work,vmult,ntotv)
188 rdif = glsum(work,ntotv)
189 rtotg = ntotg
190 rdif = (rdif-rtotg)/rtotg
191 if (abs(rdif).gt.1e-14) then
192 if(nid.eq.0)write(*,*)'ERROR: dssum test has failed!',rdif
193 call exitt
194 endif
195 endif
196
197 call projfld_c0 ! ensure fields are contiguous
198
199C print min values
200 xxmax = glmin(xm1,ntott)
201 yymax = glmin(ym1,ntott)
202 zzmax = glmin(zm1,ntott)
203
204 vxmax = glmin(vx,ntotv)
205 vymax = glmin(vy,ntotv)
206 vzmax = glmin(vz,ntotv)
207 prmax = glmin(pr,ntot2)
208
209 ntot = nxyz1*nelfld(2)
210 ttmax = glmin(t ,ntott)
211
212 do i=1,ldimt-1
213 ntot = nxyz1*nelfld(i+2)
214 psmax(i) = glmin(T(1,1,1,1,i+1),ntot)
215 enddo
216
217 if (nio.eq.0) then
218 write(6,19) xxmax,yymax,zzmax
219 19 format(' xyz min ',5g13.5)
220 endif
221 if (nio.eq.0) then
222 write(6,20) vxmax,vymax,vzmax,prmax,ttmax
223 20 format(' uvwpt min',5g13.5)
224 endif
225 if (ldimt-1.gt.0) then
226 if (nio.eq.0) write(6,21) (psmax(i),i=1,LDIMT-1)
227 21 format(' PS min ',50g13.5)
228 endif
229
230c print max values
231 xxmax = glmax(xm1,ntott)
232 yymax = glmax(ym1,ntott)
233 zzmax = glmax(zm1,ntott)
234
235 vxmax = glmax(vx,ntotv)
236 vymax = glmax(vy,ntotv)
237 vzmax = glmax(vz,ntotv)
238 prmax = glmax(pr,ntot2)
239
240 ntot = nxyz1*nelfld(2)
241 ttmax = glmax(t ,ntott)
242
243 do i=1,ldimt-1
244 ntot = nxyz1*nelfld(i+2)
245 psmax(i) = glmax(T(1,1,1,1,i+1),ntot)
246 enddo
247
248 if (nio.eq.0) then
249 write(6,16) xxmax,yymax,zzmax
250 16 format(' xyz max ',5g13.5)
251 endif
252
253 if (nio.eq.0) then
254 write(6,17) vxmax,vymax,vzmax,prmax,ttmax
255 17 format(' uvwpt max',5g13.5)
256 endif
257
258 if (ldimt-1.gt.0) then
259 if (nio.eq.0) then
260 write(6,18) (psmax(i),i=1,ldimt-1)
261 18 format(' PS max ',50g13.5)
262 endif
263 endif
264
265 if (iflomach .and. ifdp0dt) then
266 if (p0th.le.0) call exitti('Invalid thermodynamic pressure!$',1)
267 if (gamma0.lt.0) call exitti('Invalid gamma0!$',1)
268 endif
269
270 if (ifrest(0,jp)) then ! mesh has been read in.
271 if (nio.eq.0) write(6,*) 'Restart: recompute geom. factors.'
272 call geom_reset(1) ! recompute geometric factors
273 endif
274
275 if(nio.eq.0) then
276 write(6,*) 'done :: set initial conditions'
277 write(6,*) ' '
278 endif
279
280 return
281 end
282c-----------------------------------------------------------------------
283 subroutine slogic (iffort,ifrest,ifprsl,nfiles)
284C---------------------------------------------------------------------
285C
286C Set up logicals for initial conditions.
287C
288C---------------------------------------------------------------------
289 INCLUDE 'SIZE'
290 INCLUDE 'INPUT'
291 INCLUDE 'RESTART'
292c
293 logical iffort( ldimt1,0:lpert)
294 $ , ifrest(0:ldimt1,0:lpert)
295 $ , ifprsl( ldimt1,0:lpert)
296c
297 character*132 line,fname,cdum
298 character*2 s2
299 character*1 line1(132)
300 equivalence (line1,line)
301C
302C Default is user specified fortran function (=0 if not specified)
303C
304 nfldt = nfield
305 if (ifmhd) nfldt = nfield+1
306
307 do jp=0,npert
308 ifrest(0,jp) = .false.
309 do ifld=1,nfldt
310 iffort(ifld,jp) = .true.
311 ifrest(ifld,jp) = .false.
312 ifprsl(ifld,jp) = .false.
313 enddo
314 enddo
315
316 jp = 0
317 nfiles=0
318
319c Check for Presolve options
320
321 DO 1000 ILINE=1,15
322 LINE=INITC(ILINE)
323 CALL CAPIT(LINE,132)
324 IF (INDX1(LINE,'PRESOLV',7).NE.0) THEN
325C found a presolve request
326 CALL BLANK(INITC(ILINE),132)
327 CALL LJUST(LINE)
328 CALL CSPLIT(CDUM,LINE,' ',1)
329C
330 IF (LTRUNC(LINE,132).EQ.0) THEN
331 IF (NIO.EQ.0) WRITE(6,700)
332 700 FORMAT(/,2X,'Presolve options: ALL')
333C default - all fields are presolved.
334 DO 800 IFIELD=1,nfldt
335 ifprsl(ifield,jp) = .true.
336 iffort(ifield,jp) = .false.
337 800 CONTINUE
338 ELSE
339C check line for arguments
340C
341 LL=LTRUNC(LINE,132)
342 IF (NIO.EQ.0) WRITE(6,810) (LINE1(L),L=1,LL)
343 810 FORMAT(/,2X,'Presolve options: ',132A1)
344C
345 IF (INDX_CUT(LINE,'U',1).NE.0) THEN
346 ifprsl(1,jp) = .true.
347 iffort(1,jp) = .false.
348 ENDIF
349C
350 IF (INDX_CUT(LINE,'T',1).NE.0) THEN
351 ifprsl(2,jp) = .true.
352 iffort(2,jp) = .false.
353 ENDIF
354C
355 DO 900 IFIELD=3,NPSCAL+2
356 IP=IFIELD-2
357 WRITE(S2,901) IP
358 IF (INDX_CUT(LINE,S2,2).NE.0) THEN
359 ifprsl(ifield,jp) = .true.
360 iffort(ifield,jp) = .false.
361 ENDIF
362 900 CONTINUE
363 901 FORMAT('P',I1)
364 ENDIF
365 ENDIF
366 1000 CONTINUE
367C
368C Check for restart options
369C
370 jp = 0
371 DO 2000 ILINE=1,15
372 if (ifpert) jp=iline-1
373 LINE=INITC(ILINE)
374 IF (LTRUNC(LINE,132).NE.0) THEN
375C found a filename
376 NFILES=NFILES+1
377 INITC(NFILES)=LINE
378C
379C IF (NIO.EQ.0.AND.NFILES.EQ.1) WRITE(6,1010) LINE
380 IF (NIO.EQ.0.) WRITE(6,1010) LINE
381 1010 FORMAT(1X,'Checking restart options: ',A132)
382c IF (NID.EQ.0) WRITE(6,'(A132)') LINE
383C
384C Parse restart options
385
386 call sioflag(ndumps,fname,line)
387
388 if (ifgetx) then
389 ifrest(0,jp) = .true.
390 endif
391 if (ifgetu) then
392 iffort(1,jp) = .false.
393 ifprsl(1,jp) = .false.
394 ifrest(1,jp) = .true.
395 endif
396 if (ifgett) then
397 iffort(2,jp) = .false.
398 ifprsl(2,jp) = .false.
399 ifrest(2,jp) = .true.
400 endif
401 do 1900 ifield=3,nfldt
402c write(6,*) 'ifgetps:',(ifgtps(k),k=1,ldimt-1)
403 if (ifgtps(ifield-2)) then
404 iffort(ifield,jp) = .false.
405 ifprsl(ifield,jp) = .false.
406 ifrest(ifield,jp) = .true.
407 endif
408 1900 continue
409 endif
410 2000 continue
411
412 return
413 end
414c-----------------------------------------------------------------------
415 subroutine restart(nfiles)
416C----------------------------------------------------------------------
417C
418C (1) Open restart file(s)
419C (2) Check previous spatial discretization
420C (3) Map (K1,N1) => (K2,N2) if necessary
421C
422C nfiles > 1 has several implications:
423C
424C i. For std. run, data is taken from last file in list, unless
425C explicitly specified in argument list of filename
426C
427C ii. For MHD and perturbation cases, 1st file is for U,P,T;
428C subsequent files are for B-field or perturbation fields
429C
430C
431C----------------------------------------------------------------------
432 INCLUDE 'SIZE'
433 INCLUDE 'TOTAL'
434 INCLUDE 'RESTART'
435
436 common /inelr/ nelrr
437
438 PARAMETER (LXR=LX1+6)
439 PARAMETER (LYR=LY1+6)
440 PARAMETER (LZR=LZ1+6)
441 PARAMETER (LXYZR=LXR*LYR*LZR)
442 PARAMETER (LXYZT=LX1*LY1*LZ1*LELT)
443 PARAMETER (LPSC9=LDIMT+9)
444
445 common /scrcg/ pm1(lx1*ly1*lz1,lelv)
446 COMMON /SCRNS/ SDUMP(LXYZT,7)
447 integer mesg(40)
448
449C note, this usage of CTMP1 will be less than elsewhere if NELT ~> 9.
450 COMMON /CTMP1/ TDUMP(LXYZR,LPSC9)
451 real*4 tdump
452c
453 REAL SDMP2(LXYZT,LDIMT)
454
455c cdump comes in via PARALLEL (->TOTAL)
456
457 character*30 excoder
458 character*1 excoder1(30)
459 equivalence (excoder,excoder1)
460
461
462 character*132 fname
463 character*1 fname1(132)
464 equivalence (fname1,fname)
465
466 integer hnami (30)
467 character*132 hname
468 character*1 hname1(132)
469 equivalence (hname,hname1)
470 equivalence (hname,hnami )
471
472 CHARACTER*132 header
473
474C Local logical flags to determine whether to copy data or not.
475 logical ifok,iffmat
476 integer iposx,iposz,iposu,iposw,iposp,ipost,ipsps(ldimt1)
477
478 logical ifbytsw, if_byte_swap_test
479 real*4 bytetest
480
481c
482 ifok=.false.
483 ifbytsw = .false.
484
485 if(nfiles.lt.1) return
486
487 if(nio.eq.0) write(6,*) 'Reading checkpoint data '
488
489c use new reader (only binary support)
490 p67 = abs(param(67))
491 if (p67.eq.6.0) then
492 do ifile=1,nfiles
493 call sioflag(ndumps,fname,initc(ifile))
494 call mfi(fname,ifile)
495 enddo
496 call bcast(time,wdsize)! Sync time across processors
497 return
498 endif
499
500c use old reader (for ASCII + old binary support)
501
502 if (param(67).lt.1.0) then ! zero only. should be abs.
503 iffmat=.true. ! ascii
504 else
505 iffmat=.false. ! binary
506 endif
507
508 do 6000 ifile=1,nfiles
509 call sioflag(ndumps,fname,initc(ifile))
510 ierr = 0
511 if (nid.eq.0) then
512
513 if (iffmat) then
514 open (unit=91,file=fname,status='old',err=500)
515 else
516 len= ltrunc(fname,79)
517 call izero (hnami,20)
518 call chcopy(hname1,fname,len)
519c test for presence of file
520 open (unit=91,file=hname
521 $ ,form='unformatted',status='old',err=500)
522 close(unit=91)
523 call byte_open(hname,ierr)
524 if(ierr.ne.0) goto 500
525 ENDIF
526 ifok = .true.
527 endif
528
529 500 continue
530 call lbcast(ifok)
531 if (.not.ifok) goto 5000
532
533 ndumps = 1
534C
535C Only NODE 0 reads from the disk.
536C
537 DO 1000 IDUMP=1,NDUMPS
538
539 IF (NID.EQ.0) THEN
540 ! read header
541 if (iffmat) then
542 ierr = 2
543 if(mod(param(67),1.0).eq.0) then ! old header format
544 if(nelgt.lt.10000) then
545 read(91,91,err=10,end=10)
546 $ neltr,nxr,nyr,nzr,rstime,istepr,(excoder1(i),i=1,30)
547 91 format(4i4,1x,g13.4,i5,1x,30a1)
548 ierr=0
549 else
550 read(91,92,err=10,end=10)
551 $ neltr,nxr,nyr,nzr,rstime,istepr,(excoder1(i),i=1,30)
552 92 format(i10,3i4,1P1e18.9,i9,1x,30a1)
553 ierr=0
554 endif
555 else ! new head format
556 read(91,'(A132)',err=10,end=10) header
557 read(header,*)
558 & neltr,nxr,nyr,nzr,rstime,istepr,excoder
559 ierr=0
560 endif
561 else
562 if(mod(param(67),1.0).eq.0) then ! old header format
563 call byte_read(hnami,20,ierr)
564 if(ierr.ne.0) goto 10
565 icase = 2
566 if (nelgt.lt.10000) icase = 1
567 ipass = 0
568 93 continue ! test each possible case UGLY (7/31/07)
569 if(ipass.lt.2) then
570 ipass = ipass+1
571 if(icase.eq.1) then
572 read(hname,'(4i4,1x,g13.4,i5,1x,30a1)',
573 $ err=94,end=94)
574 $ neltr,nxr,nyr,nzr,rstime,istepr,
575 $ (excoder1(i),i=1,30)
576 goto 95
577 else
578 read(hname,'(i10,3i4,1P1e18.9,i9,1x,30a1)',
579 $ err=94,end=94)
580 $ neltr,nxr,nyr,nzr,rstime,istepr,
581 $ (excoder1(i),i=1,30)
582 goto 95
583 endif
584 94 icase = 3-icase ! toggle: 2-->1 1-->2
585 goto 93
586 else
587 ierr=2
588 goto 10
589 endif
590 95 continue
591 else ! new head format
592 call byte_read(header,20,ierr)
593 if(ierr.ne.0) goto 10
594 read(header,*)
595 & neltr,nxr,nyr,nzr,rstime,istepr,excoder
596 endif
597 call byte_read(bytetest,1,ierr)
598c call byte_read2(bytetest,1,ierr)
599 if(ierr.ne.0) goto 10
600 ifbytsw = if_byte_swap_test(bytetest,ierr)
601 if(ierr.ne.0) goto 10
602 endif
603 mesg(1) = neltr
604 mesg(2) = nxr
605 mesg(3) = nyr
606 mesg(4) = nzr
607 write(6,*) 'Read mode: ', param(67)
608 write(6,333)'neltr,nxr,nyr,nzr: ', neltr,nxr,nyr,nzr
609 333 format(A,i9,3i4)
610 call chcopy(mesg(5),excoder1,20)
611 len = 14*isize
612 endif
613 10 call err_chk(ierr,'Error reading restart header. Abort.$')
614
615 IF (IDUMP.EQ.1) THEN
616 len = 14*isize
617 call bcast(mesg,len)
618 neltr = mesg(1)
619 nxr = mesg(2)
620 nyr = mesg(3)
621 nzr = mesg(4)
622 call chcopy(excoder1,mesg(5),20)
623
624 call lbcast(ifbytsw)
625
626C Bounds checking on mapped data.
627 IF (NXR.GT.LXR) THEN
628 WRITE(6,20) NXR,lx1
629 20 FORMAT(//,2X,
630 $ 'ABORT: Attempt to map from',I3,
631 $ ' to N=',I3,'.',/,2X,
632 $ 'NEK5000 currently supports mapping from N+6 or less.'
633 $ ,/,2X,'Increase N or LXR in IC.FOR.')
634 CALL EXITT
635 ENDIF
636C
637C Figure out position of data in file "IFILE"
638C
639 NOUTS=0
640 IPOSX=0
641 IPOSY=0
642 IPOSZ=0
643 IPOSU=0
644 IPOSV=0
645 IPOSW=0
646 IPOSP=0
647 IPOST=0
648 DO 40 I=1,NPSCAL
649 IPSPS(I)=0
650 40 CONTINUE
651
652 IPS = 0
653 NPS = 0
654 DO 50 I=1, 30
655 IF (excoder1(i).EQ.'X') THEN
656 NOUTS=NOUTS + 1
657 IPOSX=NOUTS
658 NOUTS=NOUTS+1
659 IPOSY=NOUTS
660 IF (IF3D) THEN
661 NOUTS=NOUTS + 1
662 IPOSZ=NOUTS
663 ENDIF
664 ENDIF
665 IF (excoder1(i).EQ.'U') THEN
666 NOUTS=NOUTS + 1
667 IPOSU=NOUTS
668 NOUTS=NOUTS+1
669 IPOSV=NOUTS
670 IF (IF3D) THEN
671 NOUTS=NOUTS + 1
672 IPOSW=NOUTS
673 ENDIF
674 ENDIF
675 IF (excoder1(i).EQ.'P') THEN
676 NOUTS=NOUTS + 1
677 IPOSP=NOUTS
678 ENDIF
679 IF (excoder1(i).EQ.'T') THEN
680 NOUTS=NOUTS + 1
681 IPOST=NOUTS
682 ENDIF
683 IF(mod(param(67),1.0).eq.0.0) THEN
684 i1 = i1_from_char(excoder1(i))
685 if (0.lt.i1.and.i1.lt.10) then
686 if (i1.le.ldimt1) then
687 nouts=nouts + 1
688 ipsps(i1)=nouts
689 else
690 if (nid.eq.0) write(6,2) i1,i,excoder1(i)
691 2 format(2i4,a1,' PROBLEM W/ RESTART DATA')
692 endif
693 endif
694 ELSE
695 IF(excoder1(i).EQ.'S') THEN
696 READ(excoder1(i+1),'(I1)') NPS1
697 READ(excoder1(i+2),'(I1)') NPS0
698 NPS=10*NPS1 + NPS0
699 DO IS = 1, NPS
700 NOUTS=NOUTS + 1
701 IPSPS(IS)=NOUTS
702 ENDDO
703 GOTO 50
704 ENDIF
705 ENDIF
706 50 CONTINUE
707
708 IF (NPS.GT.(LDIMT-1)) THEN
709 IF (NID.EQ.0) THEN
710 WRITE(*,'(A)')
711 & 'ERROR: restart file has a NSPCAL > LDIMT'
712 WRITE(*,'(A,I2)')
713 & 'Change LDIMT in SIZE'
714 ENDIF
715 CALL EXITT
716 ENDIF
717
718 lname=ltrunc(fname,132)
719 if (nio.eq.0) write(6,61) (fname1(i),i=1,lname)
720 if (nio.eq.0) write(6,62)
721 $ iposu,iposv,iposw,iposp,ipost,nps,nouts
722 61 FORMAT(/,2X,'Restarting from file ',132A1)
723 62 FORMAT(2X,'Columns for restart data U,V,W,P,T,S,N: ',7I4)
724
725C Make sure the requested data is present in this file....
726 if (iposx.eq.0) ifgetx=.false.
727 if (iposy.eq.0) ifgetx=.false.
728 if (iposz.eq.0) ifgetz=.false.
729 if (iposu.eq.0) ifgetu=.false.
730 if (iposv.eq.0) ifgetu=.false.
731 if (iposw.eq.0) ifgetw=.false.
732 if (iposp.eq.0) ifgetp=.false.
733 if (ipost.eq.0) ifgett=.false.
734 do 65 i=1,npscal
735 if (ipsps(i).eq.0) ifgtps(i)=.false.
736 65 continue
737
738C End of restart file header evaluation.
739
740 ENDIF
741C
742C Read the error estimators
743C not supported at the moment => just do dummy reading
744C
745 ifok = .false.
746 IF(NID.EQ.0)THEN
747 if (iffmat)
748 & READ(91,'(6G11.4)',END=15)(CDUMP,I=1,NELTR)
749 ifok = .true.
750 ENDIF
751 15 continue
752 call lbcast(ifok)
753 if(.not.ifok) goto 1600
754C
755C Read the current dump, double buffer so that we can
756C fit the data on a distributed memory machine,
757C and so we won't have to read the restart file twice
758C in case of an incomplete data file.
759C
760 NXYZR = NXR*NYR*NZR
761C
762C Read the data
763C
764 nelrr = min(nelgt,neltr) ! # of elements to _really_read
765 ! why not just neltr?
766 do 200 ieg=1,nelrr
767 ifok = .false.
768 IF (NID.EQ.0) THEN
769 IF (MOD(IEG,10000).EQ.1) WRITE(6,*) 'Reading',IEG
770 IF (iffmat) THEN
771 READ(91,*,ERR=70,END=70)
772 $ ((tdump(IXYZ,II),II=1,NOUTS),IXYZ=1,NXYZR)
773 ELSE
774 do ii=1,nouts
775 call byte_read(tdump(1,II),nxyzr,ierr)
776 if(ierr.ne.0) then
777 write(6,*) "Error reading xyz restart data"
778 goto 70
779 endif
780 enddo
781 ENDIF
782 IFOK=.TRUE.
783 ENDIF
784C
785C Notify other processors that we've read the data OK.
786C
787 70 continue
788 call lbcast(ifok)
789 IF (.NOT.IFOK) GOTO 1600
790C
791C MAPDMP maps data from NXR to lx1
792C (and sends data to the appropriate processor.)
793C
794C The buffer SDUMP is used so that if an incomplete dump
795C file is found (e.g. due to UNIX io buffering!), then
796C the previous read data stored in VX,VY,.., is not corrupted.
797C
798 IF (IFGETX) CALL MAPDMP
799 $ (SDUMP(1,1),TDUMP(1,IPOSX),IEG,NXR,NYR,NZR,ifbytsw)
800 IF (IFGETX) CALL MAPDMP
801 $ (SDUMP(1,2),TDUMP(1,IPOSY),IEG,NXR,NYR,NZR,ifbytsw)
802 IF (IFGETZ) CALL MAPDMP
803 $ (SDUMP(1,3),TDUMP(1,IPOSZ),IEG,NXR,NYR,NZR,ifbytsw)
804 IF (IFGETU) CALL MAPDMP
805 $ (SDUMP(1,4),TDUMP(1,IPOSU),IEG,NXR,NYR,NZR,ifbytsw)
806 IF (IFGETU) CALL MAPDMP
807 $ (SDUMP(1,5),TDUMP(1,IPOSV),IEG,NXR,NYR,NZR,ifbytsw)
808 IF (IFGETW) CALL MAPDMP
809 $ (SDUMP(1,6),TDUMP(1,IPOSW),IEG,NXR,NYR,NZR,ifbytsw)
810 IF (IFGETP) CALL MAPDMP
811 $ (SDUMP(1,7),TDUMP(1,IPOSP),IEG,NXR,NYR,NZR,ifbytsw)
812 if (ifgett) call mapdmp
813 $ (SDMP2(1,1),TDUMP(1,IPOST),IEG,NXR,NYR,NZR,ifbytsw)
814
815C passive scalars
816 do 100 ips=1,npscal
817 if (ifgtps(ips)) call mapdmp(sdmp2(1,ips+1)
818 $ ,tdump(1,ipsps(ips)),ieg,nxr,nyr,nzr,IFBYTSW)
819 100 continue
820
821 200 continue
822C
823C Successfully read a complete field, store it.
824C
825 nerr = 0 ! Count number of elements rec'd by nid
826 do ieg=1,nelrr
827 mid = gllnid(ieg)
828 if (mid.eq.nid) nerr = nerr+1
829 enddo
830
831 nxyz2=lx2*ly2*lz2
832 nxyz1=lx1*ly1*lz1
833 ntott=nerr*nxyz1
834 ntotv=nerr*nxyz1 ! Problem for differing Vel. and Temp. counts!
835 ! for now we read nelt dataset
836
837 if (ifmhd.and.ifile.eq.2) then
838 if (ifgetu) call copy(bx,sdump(1,4),ntott)
839 if (ifgetu) call copy(by,sdump(1,5),ntott)
840 if (ifgetw) call copy(bz,sdump(1,6),ntott)
841 if (ifgetp) then
842 if (nio.eq.0) write(6,*) 'getting restart pressure'
843 if (ifsplit) then
844 call copy( pm,sdump(1,7),ntotv)
845 else
846 do iel=1,nelv
847 iiel = (iel-1)*nxyz1+1
848 call map12 (pm(1,1,1,iel),sdump(iiel,7),iel)
849 enddo
850 endif
851 endif
852 if (ifaxis.and.ifgett)
853 $ call copy(t(1,1,1,1,2),sdmp2(1,1),ntott)
854 elseif (ifpert.and.ifile.ge.2) then
855 j=ifile-1 ! pointer to perturbation field
856 if (ifgetu) call copy(vxp(1,j),sdump(1,4),ntotv)
857 if (ifgetu) call copy(vyp(1,j),sdump(1,5),ntotv)
858 if (ifgetw) call copy(vzp(1,j),sdump(1,6),ntotv)
859 if (ifgetp) then
860 if (nio.eq.0) write(6,*) 'getting restart pressure'
861 if (ifsplit) then
862 call copy(prp(1,j),sdump(1,7),ntotv)
863 else
864 do ie=1,nelv
865 ie1 = (ie-1)*nxyz1+1
866 ie2 = (ie-1)*nxyz2+1
867 call map12 (prp(ie2,j),sdump(ie1,7),ie)
868 enddo
869 endif
870 endif
871 if (ifgett) call copy(tp(1,1,j),sdmp2(1,1),ntott)
872C passive scalars
873 do ips=1,NPSCAL
874 if (ifgtps(ips))
875 $ call copy(tp(1,ips+1,j),sdmp2(1,ips+1),ntott)
876 enddo
877
878 else ! Std. Case
879 if (ifgetx) call copy(xm1,sdump(1,1),ntott)
880 if (ifgetx) call copy(ym1,sdump(1,2),ntott)
881 if (ifgetz) call copy(zm1,sdump(1,3),ntott)
882 if (ifgetu) call copy(vx ,sdump(1,4),ntotv)
883 if (ifgetu) call copy(vy ,sdump(1,5),ntotv)
884 if (ifgetw) call copy(vz ,sdump(1,6),ntotv)
885 if (ifgetp) call copy(pm1,sdump(1,7),ntotv)
886 if (ifgett) call copy(t,sdmp2(1,1),ntott)
887C passive scalars
888 do i=1,NPSCAL
889 if (ifgtps(i))
890 $ call copy(t(1,1,1,1,i+1),sdmp2(1,i+1),ntott)
891 enddo
892
893 if (ifaxis) call axis_interp_ic(pm1) ! Interpolate to axi mesh
894 if (ifgetp) call map_pm1_to_pr(pm1,ifile) ! Interpolate pressure
895
896 if (ifgtim) time=rstime
897 endif
898
899 1000 CONTINUE
900 GOTO 1600
901
902 1600 CONTINUE
903C
904 IF (IDUMP.EQ.1.AND.NID.EQ.0) THEN
905 write(6,1700) fname
906 write(6,1701) ieg,ixyz
907 write(6,1702)
908 $ ((tdump(jxyz,ii),ii=1,nouts),jxyz=ixyz-1,ixyz)
909 1700 FORMAT(5X,'WARNING: No data read in for file ',A132)
910 1701 FORMAT(5X,'Failed on element',I4,', point',I5,'.')
911 1702 FORMAT(5X,'Last read dump:',/,5G15.7)
912 write(6,*) nid,'call exitt 1702a',idump
913 call exitt
914 ELSEIF (IDUMP.EQ.1) THEN
915 write(6,*) nid,'call exitt 1702b',idump
916 call exitt
917 ELSE
918 IDUMP=IDUMP-1
919 IF (NIO.EQ.0) WRITE(6,1800) IDUMP
920 1800 FORMAT(2X,'Successfully read data from dump number',I3,'.')
921 ENDIF
922 if (iffmat) then
923 if (nid.eq.0) close(unit=91)
924 else
925 ierr = 0
926 if (nid.eq.0) call byte_close(ierr)
927 call err_chk(ierr,'Error closing restart file in restart$')
928 endif
929 GOTO 6000
930
931C Can't open file...
932 5000 CONTINUE
933 if (nid.eq.0) write(6,5001) fname
934 5001 FORMAT(2X,' ******* ERROR ******* '
935 $ ,/,2X,' ******* ERROR ******* '
936 $ ,/,2X,' Could not open restart file:'
937 $ ,/,A132
938 $ ,//,2X,'Quitting in routine RESTART.')
939 CLOSE(UNIT=91)
940 call exitt
941 5002 CONTINUE
942 IF (NIO.EQ.0) WRITE(6,5001) HNAME
943 call exitt
944C
945C
946C End of IFILE loop
947 6000 CONTINUE
948C
949 return
950 end
951C
952c-----------------------------------------------------------------------
953 subroutine sioflag(ndumps,fname,rsopts)
954C
955C Set IO flags according to Restart Options File, RSOPTS
956C
957 INCLUDE 'SIZE'
958 INCLUDE 'INPUT'
959 INCLUDE 'RESTART'
960 INCLUDE 'TSTEP'
961
962 character*132 rsopts,fname
963 character*2 s2
964 logical ifgtrl
965
966C Scratch variables..
967 logical ifdeft,ifanyc
968 CHARACTER*132 RSOPT ,LINE
969 CHARACTER*1 RSOPT1(132),LINE1(132)
970 EQUIVALENCE (RSOPT1,RSOPT)
971 EQUIVALENCE (LINE1,LINE)
972C
973C Parse filename
974C
975C CSPLIT splits S1 into two parts, delimited by S2.
976C S1 returns with 2nd part of S1. CSPLIT returns 1st part.
977C
978 rsopt=rsopts
979 call ljust(rsopt)
980 call csplit(fname,rsopt,' ',1)
981C check fname for user supplied extension.
982 if (indx1(fname,'.',1).eq.0) then
983 len=ltrunc(fname,132)
984 len1=len+1
985 len4=len+4
986 fname(len1:len4)='.fld'
987 endif
988C
989C Parse restart options
990C
991C set default flags
992C
993 ifgetx=.false.
994 ifgetz=.false.
995 ifgetu=.false.
996 ifgetw=.false.
997 ifgetp=.false.
998 ifgett=.false.
999 do 100 i=1,ldimt-1
1000 ifgtps(i)=.false.
1001 100 continue
1002 ifgtim=.true.
1003 ndumps=0
1004C
1005C Check for default case - just a filename given, no i/o options specified
1006C
1007 ifdeft=.true.
1008C
1009C Parse file for i/o options and/or dump number
1010C
1011 CALL CAPIT(RSOPT,132)
1012
1013 IF (LTRUNC(RSOPT,132).NE.0) THEN
1014C
1015C Check for explicit specification of restart TIME.
1016C
1017 ITO=INDX_CUT(RSOPT,'TIME',4)
1018 IFGTIM=.TRUE.
1019 IF (ITO.NE.0) THEN
1020C user has specified the time explicitly.
1021 IT1=INDX_CUT(RSOPT,'=',1)
1022 IT8=132-IT1
1023 CALL BLANK(LINE,132)
1024 CALL CHCOPY(LINE,RSOPT1(IT1),IT8)
1025 IF (IFGTRL(TTIME,LINE)) THEN
1026 IFGTIM=.FALSE.
1027 TIME=TTIME
1028 ENDIF
1029C remove the user specified time from the RS options line.
1030 ITA=132-ITO+1
1031 CALL BLANK(RSOPT1(ITO),ITA)
1032 CALL LJUST(LINE)
1033 IT1=INDX1(LINE,' ',1)
1034 ITB=132-IT1+1
1035 CALL CHCOPY(RSOPT1(ITO),LINE1(IT1),ITB)
1036 ENDIF
1037
1038C Parse field specifications.
1039
1040 IXO=INDX_CUT(RSOPT,'X',1)
1041 IF (IXO.NE.0) THEN
1042 ifdeft=.false.
1043 IFGETX=.TRUE.
1044 IF (IF3D) IFGETZ=.TRUE.
1045 ENDIF
1046
1047 IVO=INDX_CUT(RSOPT,'U',1)
1048 IF (IVO.NE.0) THEN
1049 ifdeft=.false.
1050 IFGETU=.TRUE.
1051 IF (IF3D) IFGETW=.TRUE.
1052 ENDIF
1053
1054 IPO=INDX_CUT(RSOPT,'P',1)
1055 IF (IPO.NE.0) THEN
1056 ifdeft=.false.
1057 IFGETP=.TRUE.
1058 ENDIF
1059
1060 ITO=INDX_CUT(RSOPT,'T',1)
1061 IF (ITO.NE.0) THEN
1062 ifdeft=.false.
1063 ifgett=.true.
1064 ENDIF
1065
1066 do 300 i=1,ldimt-1
1067 write (s2,301) i
1068 ipo=indx_cut(rsopt,s2,2)
1069 if (ipo.ne.0) then
1070 ifdeft=.false.
1071 ifgtps(i)=.true.
1072 endif
1073 300 continue
1074 301 format('S',i1)
1075
1076C Get number of dumps from remainder of user supplied line.
1077 if (ifgtrl(tdumps,rsopt)) ndumps=int(tdumps)
1078 endif
1079
1080C If no fields were explicitly specified, assume getting all fields.
1081 if (ifdeft) then
1082 IFGETX=.TRUE.
1083 IF (IF3D) IFGETZ=.TRUE.
1084 IFANYC=.FALSE.
1085 DO 400 I=1,NFIELD
1086 IF (IFADVC(I)) IFANYC=.TRUE.
1087 400 CONTINUE
1088 IF (IFFLOW.OR.IFANYC) THEN
1089 IFGETU=.TRUE.
1090 IF (IF3D) IFGETW=.TRUE.
1091 ENDIF
1092 if (ifflow) ifgetp=.true.
1093 if (ifheat) ifgett=.true.
1094#ifdef CMTNEK
1095 ifgett=.true. ! CMT-nek still not compatible with IFHEAT
1096#endif
1097 do 410 i=1,ldimt-1
1098 ifgtps(i)=.TRUE.
1099 410 continue
1100 endif
1101
1102
1103
1104
1105 return
1106 END
1107c-----------------------------------------------------------------------
1108 subroutine mapdmp(sdump,tdump,ieg,nxr,nyr,nzr,if_byte_sw)
1109C----------------------------------------------------------------------
1110C
1111C----------------------------------------------------------------------
1112 INCLUDE 'SIZE'
1113 INCLUDE 'PARALLEL'
1114C
1115 PARAMETER (LXYZ1=LX1*LY1*LZ1)
1116 PARAMETER (LXR=LX1+6)
1117 PARAMETER (LYR=LY1+6)
1118 PARAMETER (LZR=LZ1+6)
1119 PARAMETER (LXYZR=LXR*LYR*LZR)
1120C
1121 REAL SDUMP(LXYZ1,LELT)
1122 REAL*4 TDUMP(LXYZR)
1123c
1124 logical if_byte_sw
1125c
1126 NXYZ=lx1*ly1*lz1
1127 NXYR=NXR*NYR*NZR
1128 ierr=0
1129C
1130C Serial processor code:
1131C
1132 IF (NP.EQ.1) THEN
1133
1134 IF (if_byte_sw) call byte_reverse(TDUMP,NXYR,ierr)
1135 if(ierr.ne.0) call exitti("Error in mapdmp")
1136 IF (NXR.EQ.lx1.AND.NYR.EQ.ly1.AND.NZR.EQ.lz1) THEN
1137 CALL COPY4r(SDUMP(1,IEG),TDUMP,NXYZ)
1138 ELSE
1139C do the map (assumes that NX=NY=NZ, or NX=NY, NZ=1)
1140 call mapab4r(sdump(1,ieg),tdump,nxr,1)
1141 ENDIF
1142
1143 ELSE
1144C
1145C Parallel code - send data to appropriate processor and map.
1146C
1147 JNID=GLLNID(IEG)
1148 MTYPE=3333+GLLEL(IEG)
1149 LEN=4*NXYR
1150 LE1=4
1151 IF (NID.EQ.0.AND.JNID.NE.0) THEN
1152c hand-shake
1153 CALL CSEND(MTYPE,TDUMP,LE1,JNID,NULLPID)
1154 CALL CRECV2(MTYPE,dummy,LE1,JNID)
1155 CALL CSEND(MTYPE,TDUMP,LEN,JNID,NULLPID)
1156 ELSEIF (NID.NE.0.AND.JNID.EQ.NID) THEN
1157C Receive data from node 0
1158 CALL CRECV2(MTYPE,dummy,LE1,0)
1159 CALL CSEND(MTYPE,TDUMP,LE1,0,NULLPID)
1160 CALL CRECV2(MTYPE,TDUMP,LEN,0)
1161 ENDIF
1162C
1163C If the data is targeted for this processor, then map
1164C to appropriate element.
1165C
1166 IF (JNID.EQ.NID) THEN
1167 IE=GLLEL(IEG)
1168 IF (if_byte_sw) call byte_reverse(TDUMP,NXYR,ierr)
1169 IF (NXR.EQ.lx1.AND.NYR.EQ.ly1.AND.NZR.EQ.lz1) THEN
1170 CALL COPY4r(SDUMP(1,IE),TDUMP,NXYZ)
1171 ELSE
1172 call mapab4r(sdump(1,ie),tdump,nxr,1)
1173 ENDIF
1174 ENDIF
1175 call err_chk(ierr,'Error using byte_reverse in mapdmp.$')
1176C
1177C End of parallel distribution/map routine.
1178C
1179 ENDIF
1180 return
1181 END
1182c-----------------------------------------------------------------------
1183 subroutine mapab(x,y,nxr,nel)
1184C---------------------------------------------------------------
1185C
1186C Interpolate Y(NXR,NYR,NZR,NEL) to X(lx1,ly1,lz1,NEL)
1187C (assumes that NXR=NYR=NZR, or NXR=NYR, NZR=1)
1188C---------------------------------------------------------------
1189C
1190 INCLUDE 'SIZE'
1191 INCLUDE 'IXYZ'
1192 INCLUDE 'WZ'
1193C
1194 PARAMETER (LXR=LX1+6)
1195 PARAMETER (LYR=LY1+6)
1196 PARAMETER (LZR=LZ1+6)
1197 PARAMETER (LXYZR=LXR*LYR*LZR)
1198 PARAMETER (LXYZ1=LX1*LY1*LZ1)
1199 DIMENSION X(lx1,ly1,lz1,NEL)
1200 DIMENSION Y(NXR,NXR,NXR,NEL)
1201
1202 common /ctmp0/ xa(lxyzr) ,xb(lx1,ly1,lzr) ,xc(lxyzr)
1203 $ , zgmr(lxr) ,wgtr(lxr)
1204 common /ctmpab/ ires(lxr,lxr) ,itres(lxr,lxr)
1205 real ires,itres
1206
1207 INTEGER NOLD
1208 SAVE NOLD
1209 DATA NOLD /0/
1210
1211C Bounds checking on mapped data.
1212 if (nxr.gt.lxr) then
1213 if (nid.eq.0) write(6,20) nxr,lx1
1214 20 FORMAT(//,2X,
1215 $ 'ABORT: Attempt to map from',I3,
1216 $ ' to N=',I3,'.',/,2X,
1217 $ 'NEK5000 currently supports mapping from N+6 or less.'
1218 $ ,/,2X,'Increase N')
1219 call exitt
1220 endif
1221
1222 NZR = NXR
1223 IF(lz1.EQ.1) NZR=1
1224 NYZR = NXR*NZR
1225 NXY1 = lx1*ly1
1226
1227 IF (NXR.NE.NOLD) THEN
1228 NOLD=NXR
1229 CALL ZWGLL (ZGMR,WGTR,NXR)
1230 CALL IGLLM (IRES,ITRES,ZGMR,ZGM1,NXR,lx1,NXR,lx1)
1231 IF (NIO.EQ.0) WRITE(6,10) NXR,lx1
1232 10 FORMAT(2X,'Mapping restart data from Nold=',I2
1233 $ ,' to Nnew=',I2,'.')
1234 ENDIF
1235C
1236 DO 1000 IE=1,NEL
1237 CALL MXM (IRES,lx1,Y(1,1,1,IE),NXR,XA,NYZR)
1238 DO 100 IZ=1,NZR
1239 IZOFF = 1 + (IZ-1)*lx1*NXR
1240 CALL MXM (XA(IZOFF),lx1,ITRES,NXR,XB(1,1,IZ),ly1)
1241 100 CONTINUE
1242 IF (ldim.EQ.3) THEN
1243 CALL MXM (XB,NXY1,ITRES,NZR,X(1,1,1,IE),lz1)
1244 ELSE
1245 CALL COPY(X(1,1,1,IE),XB,NXY1)
1246 ENDIF
1247 1000 CONTINUE
1248C
1249 return
1250 END
1251c-----------------------------------------------------------------------
1252 subroutine mapab4R(x,y,nxr,nel)
1253C---------------------------------------------------------------
1254C
1255C Interpolate Y(NXR,NYR,NZR,NEL) to X(lx1,ly1,lz1,NEL)
1256C (assumes that NXR=NYR=NZR, or NXR=NYR, NZR=1)
1257c
1258c Input: real*4, Output: default precision
1259c
1260C---------------------------------------------------------------
1261C
1262 INCLUDE 'SIZE'
1263 INCLUDE 'IXYZ'
1264 INCLUDE 'WZ'
1265C
1266 PARAMETER (LXR=LX1+6)
1267 PARAMETER (LYR=LY1+6)
1268 PARAMETER (LZR=LZ1+6)
1269 PARAMETER (LXYZR=LXR*LYR*LZR)
1270 PARAMETER (LXYZ1=LX1*LY1*LZ1)
1271 REAL*4 X(lx1,ly1,lz1,NEL)
1272 REAL Y(NXR,NXR,NXR,NEL)
1273
1274 common /ctmp0/ xa(lxyzr) ,xb(lx1,ly1,lzr) ,xc(lxyzr)
1275 $ , zgmr(lxr) ,wgtr(lxr)
1276 common /ctmpa4/ ires(lxr,lxr) ,itres(lxr,lxr)
1277 real ires,itres
1278
1279 INTEGER NOLD
1280 SAVE NOLD
1281 DATA NOLD /0/
1282
1283C Bounds checking on mapped data.
1284 if (nxr.gt.lxr) then
1285 if (nid.eq.0) write(6,20) nxr,lx1
1286 20 FORMAT(//,2X,
1287 $ 'ABORT: Attempt to map from',I3,
1288 $ ' to N=',I3,'.',/,2X,
1289 $ 'NEK5000 currently supports mapping from N+6 or less.'
1290 $ ,/,2X,'Increase N')
1291 call exitt
1292 endif
1293
1294 NZR = NXR
1295 IF(lz1.EQ.1) NZR=1
1296 NYZR = NXR*NZR
1297 NXY1 = lx1*ly1
1298 nxyzr = nxr*nxr*nzr
1299
1300 IF (NXR.NE.NOLD) THEN
1301 NOLD=NXR
1302 CALL ZWGLL (ZGMR,WGTR,NXR)
1303 CALL IGLLM (IRES,ITRES,ZGMR,ZGM1,NXR,lx1,NXR,lx1)
1304 IF (NIO.EQ.0) WRITE(6,10) NXR,lx1
1305 10 FORMAT(2X,'Mapping restart data from Nold=',I2
1306 $ ,' to Nnew=',I2,'.')
1307 ENDIF
1308C
1309 DO 1000 IE=1,NEL
1310 call copy4r(xc,y(1,1,1,ie),nxyzr)
1311 CALL MXM (IRES,lx1,xc,NXR,XA,NYZR)
1312 DO 100 IZ=1,NZR
1313 IZOFF = 1 + (IZ-1)*lx1*NXR
1314 CALL MXM (XA(IZOFF),lx1,ITRES,NXR,XB(1,1,IZ),ly1)
1315 100 CONTINUE
1316 IF (ldim.EQ.3) THEN
1317 CALL MXM (XB,NXY1,ITRES,NZR,X(1,1,1,IE),lz1)
1318 ELSE
1319 CALL COPY(X(1,1,1,IE),XB,NXY1)
1320 ENDIF
1321 1000 CONTINUE
1322C
1323 return
1324 END
1325c-----------------------------------------------------------------------
1326 function i1_from_char(s1)
1327 character*1 s1
1328
1329 character*10 n10
1330 save n10
1331 data n10 / '0123456789' /
1332
1333 i1_from_char = indx2(n10,10,s1,1)-1
1334
1335 return
1336 end
1337c-----------------------------------------------------------------------
1338 integer function indx2(s1,l1,s2,l2)
1339 character*132 s1,s2
1340
1341 n1=l1-l2+1
1342 indx2=0
1343 if (n1.lt.1) return
1344
1345 do i=1,n1
1346 i2=i+l2-1
1347 if (s1(i:i2).eq.s2(1:l2)) then
1348 indx2=i
1349 return
1350 endif
1351 enddo
1352
1353 return
1354 end
1355c-----------------------------------------------------------------------
1356 INTEGER FUNCTION INDX1(S1,S2,L2)
1357 CHARACTER*132 S1,S2
1358C
1359 N1=132-L2+1
1360 INDX1=0
1361 IF (N1.LT.1) return
1362C
1363 DO 100 I=1,N1
1364 I2=I+L2-1
1365 IF (S1(I:I2).EQ.S2(1:L2)) THEN
1366 INDX1=I
1367 return
1368 ENDIF
1369 100 CONTINUE
1370C
1371 return
1372 END
1373c-----------------------------------------------------------------------
1374 INTEGER FUNCTION INDX_CUT(S1,S2,L2)
1375C
1376C INDX_CUT is returned with the location of S2 in S1 (0 if not found)
1377C S1 is returned with 1st occurance of S2 removed.
1378C
1379 CHARACTER*1 S1(132),S2(132)
1380C
1381 I1=INDX1(S1,S2,L2)
1382C
1383 IF (I1.NE.0) THEN
1384C
1385 N1=132-L2
1386 DO 100 I=I1,N1
1387 I2=I+L2
1388C remove the 1st occurance of S2 from S1.
1389 S1(I)=S1(I2)
1390 100 CONTINUE
1391 N2=N1+1
1392 DO 200 I=N2,132
1393 S1(I)=' '
1394 200 CONTINUE
1395 ENDIF
1396C
1397 INDX_CUT=I1
1398 return
1399 END
1400c-----------------------------------------------------------------------
1401 subroutine csplit(s0,s1,s2,l0)
1402 CHARACTER*132 S0,S1,S2
1403C split string S1 into two parts, delimited by S2.
1404C
1405 I2=INDX_CUT(S1,S2,L0)
1406 IF (I2.EQ.0) return
1407C
1408 I1=I2-1
1409 CALL BLANK(S0,132)
1410 S0(1:I1)=S1(1:I1)
1411 CALL LSHFT(S1,I2)
1412C
1413 return
1414 END
1415c-----------------------------------------------------------------------
1416 subroutine lshft(string,ipt)
1417C shift string from IPT to the left
1418C INPUT : "abcde...... test "
1419C OUTPUT: "e...... test " if ipt.eq.5
1420 CHARACTER*1 STRING(132)
1421C
1422 DO 20 J=1,133-IPT
1423 IJ=IPT+J-1
1424 STRING(J)=STRING(IJ)
1425 20 CONTINUE
1426 DO 30 J=134-IPT,132
1427 STRING(J)=' '
1428 30 CONTINUE
1429 return
1430 END
1431c-----------------------------------------------------------------------
1432 subroutine ljust(string)
1433C left justify string
1434 CHARACTER*1 STRING(132)
1435C
1436 IF (STRING(1).NE.' ') return
1437C
1438 DO 100 I=2,132
1439C
1440 IF (STRING(I).NE.' ') THEN
1441 DO 20 J=1,133-I
1442 IJ=I+J-1
1443 STRING(J)=STRING(IJ)
1444 20 CONTINUE
1445 DO 30 J=134-I,132
1446 STRING(J)=' '
1447 30 CONTINUE
1448 return
1449 ENDIF
1450C
1451 100 CONTINUE
1452 return
1453 END
1454c-----------------------------------------------------------------------
1455 subroutine chknorm (ifzero)
1456C----------------------------------------------------------------------
1457C
1458C Check if trivial user specified initial conditions
1459C
1460C----------------------------------------------------------------------
1461 INCLUDE 'SIZE'
1462 INCLUDE 'INPUT'
1463 INCLUDE 'TSTEP'
1464 LOGICAL IFZERO
1465C
1466 IFZERO = .TRUE.
1467C
1468 IF (IFFLOW) THEN
1469 IFIELD = 1
1470 IMESH = 1
1471 CALL UNORM
1472 IF (VNRML8.GT.0.) IFZERO = .FALSE.
1473 ENDIF
1474 IF (IFHEAT) THEN
1475 DO 100 IFIELD=2,NFIELD
1476 IMESH = 1
1477 IF (IFTMSH(IFIELD)) IMESH = 2
1478 CALL UNORM
1479 IF (TNRML8(IFIELD).GT.0.) IFZERO = .FALSE.
1480 100 CONTINUE
1481 ENDIF
1482c
1483 return
1484 END
1485C
1486c-----------------------------------------------------------------------
1487 subroutine prsolvt
1488C----------------------------------------------------------------------
1489C
1490C Use steady state solution as initial condition
1491C for temperatur/passive scalar
1492C
1493C----------------------------------------------------------------------
1494 INCLUDE 'SIZE'
1495 INCLUDE 'INPUT'
1496 INCLUDE 'TSTEP'
1497 LOGICAL IFSAV1,IFSAV2(LDIMT1)
1498C
1499 IF (NIO.EQ.0) WRITE(6,*) ' '
1500 IF (NIO.EQ.0) WRITE(6,*) 'Conduction pre-solver activated'
1501C
1502C Set logical IFTRAN to false (steady state)
1503C Save logicals for convection
1504C Turn convection off
1505C
1506 IFSAV1 = IFTRAN
1507 IFTRAN = .FALSE.
1508 DO 100 IFIELD=2,NFIELD
1509 IFSAV2(IFIELD) = IFADVC(IFIELD)
1510 IFADVC(IFIELD) = .FALSE.
1511 100 CONTINUE
1512C
1513 CALL SETPROP
1514 CALL SETSOLV
1515C
1516 IF(NIO.EQ.0)WRITE(6,*)'Steady conduction/passive scalar problem'
1517C
1518 DO 200 IGEOM=1,2
1519 CALL HEAT (IGEOM)
1520 200 CONTINUE
1521C
1522C Set IFTRAN to true again
1523C Turn convection on again
1524C
1525 IFTRAN = IFSAV1
1526 DO 300 IFIELD=2,NFIELD
1527 IFADVC(IFIELD) = IFSAV2(IFIELD)
1528 300 CONTINUE
1529C
1530 return
1531 END
1532C
1533c-----------------------------------------------------------------------
1534 subroutine prsolvv
1535C----------------------------------------------------------------------
1536C
1537C Use steady Stokes solution as initial condition
1538C for flow problem
1539C
1540C----------------------------------------------------------------------
1541 INCLUDE 'SIZE'
1542 INCLUDE 'INPUT'
1543 INCLUDE 'SOLN'
1544 INCLUDE 'TSTEP'
1545 LOGICAL IFSAV1,IFSAV2
1546C
1547 IF (NIO.EQ.0) WRITE(6,*) ' '
1548 IF (NIO.EQ.0) WRITE(6,*) 'Velocity pre-solver activated'
1549C
1550C Initialize velocity to some non-trivial RHS to avoid FP trap in i860.
1551C
1552 IF (PARAM(60).NE.0.0) THEN
1553 SMALL=10.0E-10
1554 CALL CFILL(VX,SMALL,NTOTV)
1555 CALL CFILL(VY,SMALL,NTOTV)
1556 CALL CFILL(VZ,SMALL,NTOTV)
1557 ENDIF
1558C
1559C Set logical IFTRAN to false (steady state)
1560C Save logicals for convection
1561C Turn convection off
1562C
1563 IF (IFSPLIT) THEN
1564 WRITE(6,10)
1565 10 FORMAT(
1566 $ /,2X,'ERROR: Steady Stokes Flow initial condition cannot'
1567 $,/,2X,' be computed using the splitting formulation.'
1568 $,/,2X,' Either compute using UZAWA, or remove PRESOLVE'
1569 $,/,2X,' request for velocity.'
1570 $,/,2X
1571 $,/,2X,' ABORTING IN PRSOLVV.')
1572 CALL EXITT
1573 ENDIF
1574C
1575 IFSAV1 = IFTRAN
1576 IFSAV2 = IFADVC(IFIELD)
1577 IFTRAN = .FALSE.
1578 IFADVC(IFIELD) = .FALSE.
1579C
1580 CALL SETPROP
1581 CALL SETSOLV
1582C
1583 IF (NIO.EQ.0) WRITE (6,*) 'Steady Stokes problem'
1584 DO 100 IGEOM=1,2
1585 IF (.NOT.IFSPLIT) CALL FLUID (IGEOM)
1586 100 CONTINUE
1587C
1588C Set IFTRAN to true again
1589C Turn convection on again
1590C
1591 IFTRAN = IFSAV1
1592 IFADVC(IFIELD) = IFSAV2
1593C
1594 return
1595 END
1596C
1597c-----------------------------------------------------------------------
1598 subroutine nekuic ! user specified fortran function (=0 if not specified)
1599
1600 include 'SIZE'
1601 include 'INPUT'
1602 include 'SOLN'
1603 include 'TSTEP'
1604 include 'PARALLEL'
1605 include 'NEKUSE'
1606
1607 integer e,eg
1608
1609 nel = nelfld(ifield)
1610
1611 do e=1,nel
1612 eg = lglel(e)
1613 do 300 k=1,lz1
1614 do 300 j=1,ly1
1615 do 300 i=1,lx1
1616 call nekasgn (i,j,k,e)
1617 call useric (i,j,k,eg)
1618 if (jp.eq.0) then
1619 if (ifield.eq.1) then
1620 vx(i,j,k,e) = ux
1621 vy(i,j,k,e) = uy
1622 vz(i,j,k,e) = uz
1623 elseif (ifield.eq.ifldmhd .and. ifmhd) then
1624 bx(i,j,k,e) = ux
1625 by(i,j,k,e) = uy
1626 bz(i,j,k,e) = uz
1627 else
1628 t(i,j,k,e,ifield-1) = temp
1629 endif
1630 else
1631 ijke = i+lx1*((j-1)+ly1*((k-1) + lz1*(e-1)))
1632 if (ifield.eq.1) then
1633 vxp(ijke,JP) = ux
1634 vyp(ijke,JP) = uy
1635 vzp(ijke,JP) = uz
1636 else
1637 tp(ijke,ifield-1,JP) = temp
1638 endif
1639 endif
1640
1641 300 continue
1642 enddo
1643
1644 return
1645 END
1646c-----------------------------------------------------------------------
1647 subroutine capit(lettrs,n)
1648C Capitalizes string of length n
1649 CHARACTER LETTRS(N)
1650C
1651 DO 5 I=1,N
1652 INT=ICHAR(LETTRS(I))
1653 IF(INT.GE.97 .AND. INT.LE.122) THEN
1654 INT=INT-32
1655 LETTRS(I)=CHAR(INT)
1656 ENDIF
16575 CONTINUE
1658 return
1659 END
1660c-----------------------------------------------------------------------
1661 LOGICAL FUNCTION IFGTRL(VALUE,LINE)
1662C
1663C Read VALUE from LINE and set IFGTRL to .TRUE. if successful,
1664C IFGTRL to .FALSE. otherwise.
1665C
1666C This complicated function is necessary thanks to the Ardent,
1667C which won't allow free formatted reads (*) from internal strings!
1668C
1669 CHARACTER*132 LINE
1670 CHARACTER*132 WORK
1671 CHARACTER*8 FMAT
1672C
1673C Note that the format Fn.0 is appropriate for fields of type:
1674C 34 34.0 34.0e+00
1675C The only difficulty would be with '34' but since we identify
1676C the field width exactly, there is no problem.
1677C
1678 IFGTRL=.FALSE.
1679 VALUE=0.0
1680C
1681 WORK=LINE
1682 CALL LJUST(WORK)
1683 IFLDW=INDX1(WORK,' ',1)-1
1684C
1685 IF (IFLDW.GT.0) THEN
1686 WRITE(FMAT,10) IFLDW
1687 10 FORMAT('(F',I3.3,'.0)')
1688 READ(WORK,FMAT,ERR=100,END=100) TVAL
1689 VALUE=TVAL
1690 IFGTRL=.TRUE.
1691 return
1692 ENDIF
1693C
1694 100 CONTINUE
1695 return
1696 END
1697c-----------------------------------------------------------------------
1698 LOGICAL FUNCTION IFGTIL(IVALUE,LINE)
1699C
1700C Read IVALUE from LINE and set IFGTIL to .TRUE. if successful,
1701C IFGTIL to .FALSE. otherwise.
1702C
1703C This complicated function is necessary thanks to the Ardent,
1704C which won't allow free formatted reads (*) from internal strings!
1705C
1706 CHARACTER*132 LINE
1707 CHARACTER*132 WORK
1708 CHARACTER*8 FMAT
1709C
1710 IFGTIL=.FALSE.
1711 IVALUE=0
1712C
1713 WORK=LINE
1714 CALL LJUST(WORK)
1715 IFLDW=INDX1(WORK,' ',1)-1
1716C
1717 IF (IFLDW.GT.0) THEN
1718 WRITE(FMAT,10) IFLDW
1719 10 FORMAT('(F',I3.3,'.0)')
1720 READ(WORK,FMAT,ERR=100,END=100) TVAL
1721 IVALUE=INT(TVAL)
1722 IFGTIL=.TRUE.
1723 return
1724 ENDIF
1725C
1726 100 CONTINUE
1727 return
1728 END
1729c-----------------------------------------------------------------------
1730 subroutine perturb(tt,ifld,eps)
1731 include 'SIZE'
1732 include 'TOTAL'
1733c
1734 real tt(1)
1735 integer ifld
1736
1737 ifield = ifld
1738
1739 n = lx1*ly1*lz1*nelfld(ifield)
1740 call vcospf(tt,bm1,n)
1741 call cmult(tt,eps,n)
1742 call dssum(tt,lx1,ly1,lz1)
1743
1744 return
1745 end
1746c-----------------------------------------------------------------------
1747 subroutine vcospf(x,y,n)
1748 real x(1),y(1)
1749 do i=1,n
1750 x(i) = cos(1000.*y(i))
1751 enddo
1752 return
1753 end
1754c-----------------------------------------------------------------------
1755 subroutine vbyte_swap(x,n)
1756 character*1 x(0:3,1),tmp0,tmp1
1757 character*1 in (0:3), out(0:3)
1758 real*4 in4 , out4
1759 equivalence (in ,in4 )
1760 equivalence (out,out4)
1761c
1762 do i=1,n
1763 do j=0,3
1764 in (j) = x(j,i)
1765 enddo
1766 tmp0 = x(0,i)
1767 tmp1 = x(1,i)
1768 x(0,i) = x(3,i)
1769 x(1,i) = x(2,i)
1770 x(2,i) = tmp1
1771 x(3,i) = tmp0
1772 do j=0,3
1773 out(j) = x(j,i)
1774 enddo
1775 write(6,*) 'swap:',i,in4,out4
1776 enddo
1777c
1778 return
1779 end
1780c-----------------------------------------------------------------------
1781 logical function if_byte_swap_test(bytetest,ierr)
1782 include 'SIZE'
1783
1784 real*4 bytetest,test2
1785 real*4 test_pattern
1786 save test_pattern
1787
1788 test_pattern = 6.54321
1789 eps = 0.00020
1790 etest = abs(test_pattern-bytetest)
1791 if_byte_swap_test = .true.
1792 if (etest.le.eps) if_byte_swap_test = .false.
1793
1794 test2 = bytetest
1795 call byte_reverse(test2,1,ierr)
1796 if (nid.eq.0 .and. loglevel.gt.2)
1797 $ write(6,*) 'byte swap:',if_byte_swap_test,bytetest,test2
1798 return
1799 end
1800c-----------------------------------------------------------------------
1801 subroutine geom_reset(icall)
1802C
1803C Generate geometry data
1804C
1805 INCLUDE 'SIZE'
1806 INCLUDE 'INPUT'
1807 INCLUDE 'GEOM'
1808 INCLUDE 'SOLN'
1809 INCLUDE 'TSTEP'
1810 include 'WZ'
1811c
1812 COMMON /scruz/ XM3 (LX1,LY1,LZ1,LELT)
1813 $ , YM3 (LX1,LY1,LZ1,LELT)
1814 $ , ZM3 (LX1,LY1,LZ1,LELT)
1815C
1816c
1817 if(nio.eq.0) write(6,*) 'regenerate geometry data',icall
1818
1819 ntot = lx1*ly1*lz1*nelt
1820c
1821 if (lx3.eq.lx1) then
1822 call copy(xm3,xm1,ntot)
1823 call copy(ym3,ym1,ntot)
1824 call copy(zm3,zm1,ntot)
1825 else
1826 call map13_all(xm3,xm1)
1827 call map13_all(ym3,ym1)
1828 if (if3d) call map13_all(zm3,zm1)
1829 endif
1830c
1831 CALL GEOM1 (XM3,YM3,ZM3)
1832 CALL GEOM2
1833 CALL UPDMSYS (1)
1834 CALL VOLUME
1835 CALL SETINVM
1836 CALL SETDEF
1837 CALL SFASTAX
1838c
1839 do ie = 1,nelt
1840 ! x
1841 xc(1,ie) = XM1(1 ,1 ,1 ,ie)
1842 xc(2,ie) = XM1(lx1,1 ,1 ,ie)
1843 xc(3,ie) = XM1(lx1,ly1,1 ,ie)
1844 xc(4,ie) = XM1(1 ,ly1,1 ,ie)
1845 xc(5,ie) = XM1(1 ,1 ,lz1,ie)
1846 xc(6,ie) = XM1(lx1,1 ,lz1,ie)
1847 xc(7,ie) = XM1(lx1,ly1,lz1,ie)
1848 xc(8,ie) = XM1(1 ,ly1,lz1,ie)
1849 ! y
1850 yc(1,ie) = YM1(1 ,1 ,1 ,ie)
1851 yc(2,ie) = YM1(lx1,1 ,1 ,ie)
1852 yc(3,ie) = YM1(lx1,ly1,1 ,ie)
1853 yc(4,ie) = YM1(1 ,ly1,1 ,ie)
1854 yc(5,ie) = YM1(1 ,1 ,lz1,ie)
1855 yc(6,ie) = YM1(lx1,1 ,lz1,ie)
1856 yc(7,ie) = YM1(lx1,ly1,lz1,ie)
1857 yc(8,ie) = YM1(1 ,ly1,lz1,ie)
1858 ! z
1859 zc(1,ie) = ZM1(1 ,1 ,1 ,ie)
1860 zc(2,ie) = ZM1(lx1,1 ,1 ,ie)
1861 zc(3,ie) = ZM1(lx1,ly1,1 ,ie)
1862 zc(4,ie) = ZM1(1 ,ly1,1 ,ie)
1863 zc(5,ie) = ZM1(1 ,1 ,lz1,ie)
1864 zc(6,ie) = ZM1(lx1,1 ,lz1,ie)
1865 zc(7,ie) = ZM1(lx1,ly1,lz1,ie)
1866 zc(8,ie) = ZM1(1 ,ly1,lz1,ie)
1867 enddo
1868
1869 if(nio.eq.0) then
1870 write(6,*) 'done :: regenerate geometry data',icall
1871 write(6,*) ' '
1872 endif
1873
1874 return
1875 end
1876c-----------------------------------------------------------------------
1877 subroutine dsavg(u)
1878c
1879c
1880 include 'SIZE'
1881 include 'TOTAL'
1882 real u(lx1,ly1,lz1,lelt)
1883c
1884c
1885c Take direct stiffness avg of u
1886c
1887c
1888 ifieldo = ifield
1889 if (ifflow) then
1890 ifield = 1
1891 ntot = lx1*ly1*lz1*nelv
1892 call dssum(u,lx1,ly1,lz1)
1893 call col2 (u,vmult,ntot)
1894 else
1895 ifield = 2
1896 ntot = lx1*ly1*lz1*nelt
1897 call dssum(u,lx1,ly1,lz1)
1898 call col2 (u,tmult,ntot)
1899 endif
1900 ifield = ifieldo
1901c
1902 return
1903 end
1904c-----------------------------------------------------------------------
1905 subroutine map13_all(x3,x1)
1906c
1907 include 'SIZE'
1908 include 'TOTAL'
1909c
1910 real x3(lx3,ly3,lz3,lelt)
1911 real x1(lx1,ly1,lz1,lelt)
1912c
1913 integer e
1914c
1915 do e=1,nelt
1916 call map13 (x3(1,1,1,e),x1(1,1,1,e),e)
1917 enddo
1918c
1919 return
1920 end
1921c-----------------------------------------------------------------------
1922 subroutine mfi_gets(u,wk,lwk,iskip)
1923
1924 include 'SIZE'
1925 include 'INPUT'
1926 include 'PARALLEL'
1927 include 'RESTART'
1928
1929
1930 real u(lx1*ly1*lz1,1)
1931
1932 real*4 wk(2*lwk) ! message buffer
1933 parameter(lrbs=20*lx1*ly1*lz1*lelt)
1934 common /vrthov/ w2(lrbs) ! read buffer
1935 real*4 w2
1936
1937 integer e,ei,eg,msg_id(lelt)
1938 logical iskip
1939 integer*8 i8tmp
1940
1941 call nekgsync() ! clear outstanding message queues.
1942
1943 nxyzr = nxr*nyr*nzr
1944 dnxyzr = nxyzr
1945 len = nxyzr*wdsizr ! message length
1946 if (wdsizr.eq.8) nxyzr = 2*nxyzr
1947
1948 ! check message buffer wk
1949 num_recv = nxyzr*nelt
1950 num_avail = size(wk)
1951 call lim_chk(num_recv,num_avail,' ',' ','mfi_gets a')
1952
1953 ! setup read buffer
1954 if (nid.eq.pid0r) then
1955 i8tmp = int(nxyzr,8)*int(nelr,8)
1956 nread = i8tmp/int(lrbs,8)
1957 if (mod(i8tmp,int(lrbs,8)).ne.0) nread = nread + 1
1958 if(ifmpiio) nread = iglmax(nread,1) ! needed because of collective read
1959 nelrr = nelr/nread
1960 endif
1961 call bcast(nelrr,4)
1962 call lim_chk(nxyzr*nelrr,lrbs,' ',' ','mfi_gets b')
1963
1964 ! pre-post recieves
1965 if (np.gt.1) then
1966 l = 1
1967 do e=1,nelt
1968 msg_id(e) = irecv(e,wk(l),len)
1969 l = l+nxyzr
1970 enddo
1971 endif
1972
1973 ierr = 0
1974 if (nid.eq.pid0r.and.np.gt.1) then ! only i/o nodes will read
1975 ! read blocks of size nelrr
1976 k = 0
1977 do i = 1,nread
1978 if(i.eq.nread) then ! clean-up
1979 nelrr = nelr - (nread-1)*nelrr
1980 if(nelrr.lt.0) nelrr = 0
1981 endif
1982
1983 if(ierr.eq.0) then
1984 if(ifmpiio) then
1985 call byte_read_mpi(w2,nxyzr*nelrr,-1,ifh_mbyte,ierr)
1986 else
1987 call byte_read (w2,nxyzr*nelrr,ierr)
1988 endif
1989 endif
1990
1991 ! distribute data across target processors
1992 l = 1
1993 do e = k+1,k+nelrr
1994 jnid = gllnid(er(e)) ! where is er(e) now?
1995 jeln = gllel(er(e))
1996 if(ierr.ne.0) call rzero(w2(l),len)
1997 call csend(jeln,w2(l),len,jnid,0) ! blocking send
1998 l = l+nxyzr
1999 enddo
2000 k = k + nelrr
2001 enddo
2002 elseif (np.eq.1) then
2003 if(ifmpiio) then
2004 call byte_read_mpi(wk,nxyzr*nelr,-1,ifh_mbyte,ierr)
2005 else
2006 call byte_read(wk,nxyzr*nelr,ierr)
2007 endif
2008 endif
2009
2010 if (iskip) then
2011 call nekgsync() ! clear outstanding message queues.
2012 goto 100 ! don't use the data
2013 endif
2014
2015 nxyzr = nxr*nyr*nzr
2016 nxyzv = nxr*nyr*nzr
2017 nxyzw = nxr*nyr*nzr
2018 if (wdsizr.eq.8) nxyzw = 2*nxyzw
2019
2020 l = 1
2021 do e=1,nelt
2022 if (np.gt.1) then
2023 call msgwait(msg_id(e))
2024 ei = e
2025 elseif(np.eq.1) then
2026 ei = er(e)
2027 endif
2028 if (if_byte_sw) then
2029 if(wdsizr.eq.8) then
2030 call byte_reverse8(wk(l),nxyzv*2,ierr)
2031 else
2032 call byte_reverse(wk(l),nxyzv,ierr)
2033 endif
2034 endif
2035 if (nxr.eq.lx1.and.nyr.eq.ly1.and.nzr.eq.lz1) then
2036 if (wdsizr.eq.4) then ! COPY
2037 call copy4r(u(1,ei),wk(l ),nxyzr)
2038 else
2039 call copy (u(1,ei),wk(l ),nxyzr)
2040 endif
2041 else ! INTERPOLATE
2042 if (wdsizr.eq.4) then
2043 call mapab4r(u(1,ei),wk(l ),nxr,1)
2044 else
2045 call mapab (u(1,ei),wk(l ),nxr,1)
2046 endif
2047 endif
2048 l = l+nxyzw
2049 enddo
2050
2051
2052 100 call err_chk(ierr,'Error reading restart data,in gets.$')
2053 return
2054 end
2055c-----------------------------------------------------------------------
2056 subroutine mfi_getv(u,v,w,wk,lwk,iskip)
2057
2058 include 'SIZE'
2059 include 'INPUT'
2060 include 'PARALLEL'
2061 include 'RESTART'
2062
2063 real u(lx1*ly1*lz1,1),v(lx1*ly1*lz1,1),w(lx1*ly1*lz1,1)
2064 logical iskip
2065
2066 real*4 wk(2*lwk) ! message buffer
2067 parameter(lrbs=20*lx1*ly1*lz1*lelt)
2068 common /vrthov/ w2(lrbs) ! read buffer
2069 real*4 w2
2070
2071 integer e,ei,eg,msg_id(lelt)
2072 integer*8 i8tmp
2073
2074 call nekgsync() ! clear outstanding message queues.
2075
2076 nxyzr = ldim*nxr*nyr*nzr
2077 len = nxyzr*wdsizr ! message length in bytes
2078 if (wdsizr.eq.8) nxyzr = 2*nxyzr
2079
2080 ! check message buffer wk
2081 num_recv = nxyzr*nelt
2082 num_avail = size(wk)
2083 call lim_chk(num_recv,num_avail,' ',' ','mfi_getv a')
2084
2085 ! setup read buffer
2086 if(nid.eq.pid0r) then
2087 i8tmp = int(nxyzr,8)*int(nelr,8)
2088 nread = i8tmp/int(lrbs,8)
2089 if (mod(i8tmp,int(lrbs,8)).ne.0) nread = nread + 1
2090 if(ifmpiio) nread = iglmax(nread,1) ! needed because of collective read
2091 nelrr = nelr/nread
2092 endif
2093 call bcast(nelrr,4)
2094 call lim_chk(nxyzr*nelrr,lrbs,' ',' ','mfi_getv b')
2095
2096 ! pre-post recieves (one mesg per element)
2097 ! this assumes we never pre post more messages than supported
2098 if (np.gt.1) then
2099 l = 1
2100 do e=1,nelt
2101 msg_id(e) = irecv(e,wk(l),len)
2102 l = l+nxyzr
2103 enddo
2104 endif
2105
2106 ierr = 0
2107 if (nid.eq.pid0r .and. np.gt.1) then ! only i/o nodes
2108 k = 0
2109 do i = 1,nread
2110 if(i.eq.nread) then ! clean-up
2111 nelrr = nelr - (nread-1)*nelrr
2112 if(nelrr.lt.0) nelrr = 0
2113 endif
2114
2115 if(ierr.eq.0) then
2116 if(ifmpiio) then
2117 call byte_read_mpi(w2,nxyzr*nelrr,-1,ifh_mbyte,ierr)
2118 else
2119 call byte_read (w2,nxyzr*nelrr,ierr)
2120 endif
2121 endif
2122
2123 ! redistribute data based on the current el-proc map
2124 l = 1
2125 do e = k+1,k+nelrr
2126 jnid = gllnid(er(e)) ! where is er(e) now?
2127 jeln = gllel(er(e))
2128 if(ierr.ne.0) call rzero(w2(l),len)
2129 call csend(jeln,w2(l),len,jnid,0) ! blocking send
2130 l = l+nxyzr
2131 enddo
2132 k = k + nelrr
2133 enddo
2134 elseif (np.eq.1) then
2135 if(ifmpiio) then
2136 call byte_read_mpi(wk,nxyzr*nelr,-1,ifh_mbyte,ierr)
2137 else
2138 call byte_read(wk,nxyzr*nelr,ierr)
2139 endif
2140 endif
2141
2142 if (iskip) then
2143 call nekgsync() ! clear outstanding message queues.
2144 goto 100 ! don't assign the data we just read
2145 endif
2146
2147 nxyzr = nxr*nyr*nzr
2148 nxyzv = ldim*nxr*nyr*nzr
2149 nxyzw = nxr*nyr*nzr
2150 if (wdsizr.eq.8) nxyzw = 2*nxyzw
2151
2152 l = 1
2153 do e=1,nelt
2154 if (np.gt.1) then
2155 call msgwait(msg_id(e))
2156 ei = e
2157 else if(np.eq.1) then
2158 ei = er(e)
2159 endif
2160 if (if_byte_sw) then
2161 if(wdsizr.eq.8) then
2162 call byte_reverse8(wk(l),nxyzv*2,ierr)
2163 else
2164 call byte_reverse(wk(l),nxyzv,ierr)
2165 endif
2166 endif
2167 if (nxr.eq.lx1.and.nyr.eq.ly1.and.nzr.eq.lz1) then
2168 if (wdsizr.eq.4) then ! COPY
2169 call copy4r(u(1,ei),wk(l ),nxyzr)
2170 call copy4r(v(1,ei),wk(l+ nxyzw),nxyzr)
2171 if (if3d)
2172 $ call copy4r(w(1,ei),wk(l+2*nxyzw),nxyzr)
2173 else
2174 call copy (u(1,ei),wk(l ),nxyzr)
2175 call copy (v(1,ei),wk(l+ nxyzw),nxyzr)
2176 if (if3d)
2177 $ call copy (w(1,ei),wk(l+2*nxyzw),nxyzr)
2178 endif
2179 else ! INTERPOLATE
2180 if (wdsizr.eq.4) then
2181 call mapab4r(u(1,ei),wk(l ),nxr,1)
2182 call mapab4r(v(1,ei),wk(l+ nxyzw),nxr,1)
2183 if (if3d)
2184 $ call mapab4r(w(1,ei),wk(l+2*nxyzw),nxr,1)
2185 else
2186 call mapab (u(1,ei),wk(l ),nxr,1)
2187 call mapab (v(1,ei),wk(l+ nxyzw),nxr,1)
2188 if (if3d)
2189 $ call mapab (w(1,ei),wk(l+2*nxyzw),nxr,1)
2190 endif
2191 endif
2192 l = l+ldim*nxyzw
2193 enddo
2194
2195 100 call err_chk(ierr,'Error reading restart data, in getv.$')
2196 return
2197 end
2198c-----------------------------------------------------------------------
2199 subroutine mfi_parse_hdr(hdr,ierr)
2200 include 'SIZE'
2201
2202 character*132 hdr
2203
2204 if (indx2(hdr,132,'#std',4).eq.1) then
2205 call parse_std_hdr(hdr)
2206 else
2207 if (nio.eq.0) write(6,80) hdr
2208 if (nio.eq.0) write(6,80) 'NONSTD HDR, parse_hdr, abort.'
2209 80 format(a132)
2210 ierr = 1
2211 endif
2212
2213 return
2214 end
2215c-----------------------------------------------------------------------
2216 subroutine parse_std_hdr(hdr)
2217 include 'SIZE'
2218 include 'INPUT'
2219 include 'SOLN'
2220 include 'PARALLEL'
2221 include 'RESTART'
2222 include 'TSTEP'
2223
2224 character*132 hdr
2225 character*4 dummy
2226 logical if_press_mesh
2227
2228 p0thr = -1
2229 if_press_mesh = .false.
2230
2231 read(hdr,*,iostat=ierr) dummy
2232 $ , wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
2233 $ , ifiler,nfiler
2234 $ , rdcode ! 74+20=94
2235 $ , p0thr, if_press_mesh
2236
2237 if (ierr.gt.0) then ! try again without pressure format flag
2238 read(hdr,*,iostat=ierr) dummy
2239 $ , wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
2240 $ , ifiler,nfiler
2241 $ , rdcode ! 74+20=94
2242 $ , p0thr
2243 endif
2244
2245 if (ierr.gt.0) then ! try again without mean pressure
2246 read(hdr,*,err=99) dummy
2247 $ , wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
2248 $ , ifiler,nfiler
2249 $ , rdcode ! 74+20=94
2250 endif
2251
2252c set if_full_pres flag
2253 if_full_pres = .false.
2254 if (.not.ifsplit) if_full_pres = if_press_mesh
2255
2256c ifgtim = .true. ! always get time
2257 ifgetxr = .false.
2258 ifgetur = .false.
2259 ifgetpr = .false.
2260 ifgettr = .false.
2261 do k=1,ldimt-1
2262 ifgtpsr(k) = .false.
2263 enddo
2264
2265 NPSR = 0
2266 NPS = 0
2267 do i=1,10
2268 if (rdcode1(i).eq.'X') ifgetxr = .true.
2269 if (rdcode1(i).eq.'U') ifgetur = .true.
2270 if (rdcode1(i).eq.'P') ifgetpr = .true.
2271 if (rdcode1(i).eq.'T') ifgettr = .true.
2272 if (rdcode1(i).eq.'S') then
2273 read(rdcode1(i+1),'(I1)') NPS1
2274 read(rdcode1(i+2),'(I1)') NPS0
2275 NPSR = 10*NPS1+NPS0
2276 NPS = NPSR
2277 if(NPSR.gt.ldimt-1) NPS=ldimt-1
2278 do k=1,NPS
2279 ifgtpsr(k) = .true.
2280 enddo
2281 ! nothing will follow
2282 GOTO 50
2283 endif
2284 enddo
2285
2286 50 if (NPS.lt.NPSR) then
2287 if (nid.eq.0) then
2288 write(*,'(A,/,A)')
2289 & 'WARNING: restart file has a NSPCAL > LDIMT',
2290 & 'read only part of the fld-data!'
2291 endif
2292 endif
2293
2294 if (NPS.lt.NPSCAL) then
2295 if (nid.eq.0) then
2296 write(*,'(A,/,A)')
2297 & 'WARNING: NPSCAL read from restart file differs from ',
2298 & 'currently used NPSCAL!'
2299 endif
2300 endif
2301
2302 p0th = 1
2303 if (p0thr.gt.0) p0th = p0thr
2304
2305 return
2306
2307 99 continue ! If we got here, then the May 2008 variant of std hdr
2308 ! failed and we may have an older input file.
2309
2310 call parse_std_hdr_2006(hdr,rdcode) ! try the original header format
2311
2312 return
2313 end
2314c-----------------------------------------------------------------------
2315 subroutine parse_std_hdr_2006(hdr,rlcode)
2316 include 'SIZE'
2317 include 'INPUT'
2318 include 'RESTART'
2319
2320 character*132 hdr
2321 character*1 rlcode(20)
2322
2323c 4 7 10 13 23 33 53 62 68 74
2324 read(hdr,1) wdsizr,nxr,nyr,nzr,nelr,nelgr,timer,istpr
2325 $ , ifiler,nfiler
2326 $ , (rlcode(k),k=1,20) ! 74+20=94
2327 1 format(4x,i2,3i3,2i10,e20.13,i9,2i6,20a1)
2328
2329 if (nid.eq.0) write(6,*) 'WARNING: reading depreacted header!'
2330
2331c Assign read conditions, according to rdcode
2332c NOTE: In the old hdr format: what you see in file is what you get.
2333c ifgtim = .true. ! always get time
2334 ifgetxr = .false.
2335 ifgetur = .false.
2336 ifgetpr = .false.
2337 ifgettr = .false.
2338 do k=1,npscal
2339 ifgtpsr(k) = .false.
2340 enddo
2341
2342 if (rlcode(1).eq.'X') ifgetxr = .true.
2343 if (rlcode(2).eq.'U') ifgetur = .true.
2344 if (rlcode(3).eq.'P') ifgetpr = .true.
2345 if (rlcode(4).eq.'T') ifgettr = .true.
2346 do k=1,npscal
2347 if (rlcode(4+k).ne.' ') ifgtpsr(k) = .true.
2348 enddo
2349
2350
2351 return
2352 end
2353c-----------------------------------------------------------------------
2354 subroutine mfi(fname_in,ifile)
2355c
2356c (1) Open restart file(s)
2357c (2) Check previous spatial discretization
2358c (3) Map (K1,N1) => (K2,N2) if necessary
2359c
2360c nfiles > 1 has several implications:
2361c
2362c i. For std. run, data is taken from last file in list, unless
2363c explicitly specified in argument list of filename
2364c
2365c ii. For MHD and perturbation cases, 1st file is for U,P,T;
2366c subsequent files are for B-field or perturbation fields
2367c
2368c
2369 include 'SIZE'
2370 include 'TOTAL'
2371 include 'RESTART'
2372 character*132 hdr
2373 character*132 fname_in
2374
2375 character*132 fname
2376 character*1 fnam1(132)
2377 equivalence (fnam1,fname)
2378
2379 character*1 frontc
2380
2381 parameter (lwk = 7*lx1*ly1*lz1*lelt)
2382 common /scrns/ wk(lwk)
2383 common /scrcg/ pm1(lx1*ly1*lz1,lelv)
2384 integer e
2385
2386 integer*8 offs0,offs,nbyte,stride,strideB,nxyzr8
2387
2388 tiostart=dnekclock()
2389
2390 ! add full path if required
2391 call blank(fname,132)
2392 call chcopy(frontc, fname_in, 1)
2393 if (frontc .ne. '/') then
2394 lenp = ltrunc(path,132)
2395 lenf = ltrunc(fname_in,132)
2396 call chcopy(fnam1(1),path,lenp)
2397 call chcopy(fnam1(lenp+1),fname_in,lenf)
2398 else
2399 lenf = ltrunc(fname_in,132)
2400 call chcopy(fnam1(1),fname_in,lenf)
2401 endif
2402
2403 call mfi_prepare(fname) ! determine reader nodes +
2404 ! read hdr + element mapping
2405
2406 offs0 = iHeadersize + 4 + isize*nelgr
2407 nxyzr8 = nxr*nyr*nzr
2408 strideB = nelBr* nxyzr8*wdsizr
2409 stride = nelgr* nxyzr8*wdsizr
2410
2411 iofldsr = 0
2412 if (ifgetxr) then ! if available
2413 offs = offs0 + ldim*strideB
2414 call byte_set_view(offs,ifh_mbyte)
2415 if (ifgetx) then
2416c if(nid.eq.0) write(6,*) 'Reading mesh'
2417 call mfi_getv(xm1,ym1,zm1,wk,lwk,.false.)
2418 else ! skip the data
2419 call mfi_getv(xm1,ym1,zm1,wk,lwk,.true.)
2420 endif
2421 iofldsr = iofldsr + ldim
2422 endif
2423
2424 if (ifgetur) then
2425 offs = offs0 + iofldsr*stride + ldim*strideB
2426 call byte_set_view(offs,ifh_mbyte)
2427 if (ifgetu) then
2428 if (ifmhd.and.ifile.eq.2) then
2429c if(nid.eq.0) write(6,*) 'Reading B field'
2430 call mfi_getv(bx,by,bz,wk,lwk,.false.)
2431 else
2432c if(nid.eq.0) write(6,*) 'Reading velocity field'
2433 call mfi_getv(vx,vy,vz,wk,lwk,.false.)
2434 endif
2435 else
2436 call mfi_getv(vx,vy,vz,wk,lwk,.true.)
2437 endif
2438 iofldsr = iofldsr + ldim
2439 endif
2440
2441 if (ifgetpr) then
2442 offs = offs0 + iofldsr*stride + strideB
2443 call byte_set_view(offs,ifh_mbyte)
2444 if (ifgetp) then
2445c if(nid.eq.0) write(6,*) 'Reading pressure field'
2446 call mfi_gets(pm1,wk,lwk,.false.)
2447 else
2448 call mfi_gets(pm1,wk,lwk,.true.)
2449 endif
2450 iofldsr = iofldsr + 1
2451 endif
2452
2453 if (ifgettr) then
2454 offs = offs0 + iofldsr*stride + strideB
2455 call byte_set_view(offs,ifh_mbyte)
2456 if (ifgett) then
2457c if(nid.eq.0) write(6,*) 'Reading temperature field'
2458 call mfi_gets(t,wk,lwk,.false.)
2459 else
2460 call mfi_gets(t,wk,lwk,.true.)
2461 endif
2462 iofldsr = iofldsr + 1
2463 endif
2464
2465 ierr = 0
2466 do k=1,ldimt-1
2467 if (ifgtpsr(k)) then
2468 offs = offs0 + iofldsr*stride + strideB
2469 call byte_set_view(offs,ifh_mbyte)
2470 if (ifgtps(k)) then
2471c if(nid.eq.0) write(6,'(A,I2,A)') ' Reading ps',k,' field'
2472 call mfi_gets(t(1,1,1,1,k+1),wk,lwk,.false.)
2473 else
2474 call mfi_gets(t(1,1,1,1,k+1),wk,lwk,.true.)
2475 endif
2476 iofldsr = iofldsr + 1
2477 endif
2478 enddo
2479 nbyte = 0
2480 if(nid.eq.pid0r) nbyte = iofldsr*nelr*wdsizr*nxr*nyr*nzr
2481
2482 if (ifgtim) time = timer
2483
2484 if(ifmpiio) then
2485 if(nid.eq.pid0r) call byte_close_mpi(ifh_mbyte,ierr)
2486 else
2487 if(nid.eq.pid0r) call byte_close(ierr)
2488 endif
2489 call err_chk(ierr,'Error closing restart file, in mfi.$')
2490 tio = dnekclock()-tiostart
2491
2492 dnbyte = nbyte
2493 nbyte = glsum(dnbyte,1)
2494 nbyte = nbyte + iHeaderSize + 4 + isize*nelgr
2495
2496 if (tio.eq.0) tio=1
2497 if (nio.eq.0) write(6,7) istep,time,
2498 & nbyte/tio/1024/1024/10,
2499 & nfiler
2500 7 format(/,i9,1pe12.4,' done :: Read checkpoint data',/,
2501 & 30X,'avg data-throughput = ',f7.1,'MBps',/,
2502 & 30X,'io-nodes = ',i5,/)
2503
2504
2505 if (ifaxis) call axis_interp_ic(pm1) ! Interpolate to axi mesh
2506 if (ifgetp) call map_pm1_to_pr(pm1,ifile) ! Interpolate pressure
2507
2508 return
2509 end
2510c-----------------------------------------------------------------------
2511 subroutine addfid(fname,fid)
2512 include 'SIZE'
2513 include 'TSTEP'
2514 include 'INPUT'
2515 include 'RESTART'
2516
2517
2518 character*1 fname(132)
2519 integer fid
2520
2521 character*8 eight,fmt,s8
2522 save eight
2523 data eight / "????????" /
2524
2525 do ipass=1,2 ! 2nd pass, in case 1 file/directory
2526 do k=8,1,-1
2527 i1 = indx1(fname,eight,k)
2528 if (i1.ne.0) then ! found k??? string
2529 write(fmt,1) k,k
2530 1 format('(i',i1,'.',i1,')')
2531 write(s8,fmt) fid
2532 call chcopy(fname(i1),s8,k)
2533 goto 10
2534 endif
2535 enddo
2536 10 continue
2537 enddo
2538
2539 return
2540 end
2541c-----------------------------------------------------------------------
2542 subroutine mfi_prepare(hname) ! determine which nodes are readers
2543 character*132 hname
2544
2545 include 'SIZE'
2546 include 'PARALLEL'
2547 include 'RESTART'
2548 include 'INPUT'
2549
2550 integer stride
2551 character*132 hdr, hname_
2552 logical if_byte_swap_test
2553 real*4 bytetest
2554
2555 integer*8 offs0,offs
2556
2557 ierr = 0
2558 ! rank0 (i/o master) will do a pre-read to get some infos
2559 ! we need to have in advance
2560 if (nid.eq.0) then
2561 call chcopy(hname_,hname,132)
2562 call addfid(hname_,0)
2563 call byte_open(hname_,ierr)
2564
2565 if(ierr.ne.0) goto 101
2566 call blank (hdr,iHeaderSize)
2567 call byte_read (hdr,iHeaderSize/4,ierr)
2568 if(ierr.ne.0) goto 101
2569 call byte_read (bytetest,1,ierr)
2570 if(ierr.ne.0) goto 101
2571 if_byte_sw = if_byte_swap_test(bytetest,ierr) ! determine endianess
2572 if(ierr.ne.0) goto 101
2573 call byte_close(ierr)
2574 endif
2575
2576 101 continue
2577 call err_chk(ierr,'Error reading restart header in mfi_prepare$')
2578
2579 call bcast(if_byte_sw,lsize)
2580 call bcast(hdr,iHeaderSize)
2581 call mfi_parse_hdr(hdr,ierr)
2582
2583 ifmpiio = .false.
2584 if (nfiler.eq.1 .and. abs(param(67)).eq.6) ifmpiio = .true.
2585#ifdef NOMPIIO
2586 ifmpiio = .false.
2587#endif
2588
2589 if(.not.ifmpiio) then
2590
2591 stride = np / nfiler
2592 if (stride.lt.1) then
2593 write(6,*) nfiler,np,' TOO MANY FILES, mfi_prepare'
2594 call exitt
2595 endif
2596
2597 if (mod(nid,stride).eq.0) then ! i/o clients
2598 pid0r = nid
2599 pid1r = nid + stride
2600 fid0r = nid / stride
2601 call blank(hdr,iHeaderSize)
2602
2603 call addfid(hname,fid0r)
2604 if(nid.eq.pid0r) write(6,*) ' FILE:',hname
2605 call byte_open(hname,ierr)
2606
2607 if(ierr.ne.0) goto 102
2608 call byte_read (hdr, iHeaderSize/4,ierr)
2609 if(ierr.ne.0) goto 102
2610 call byte_read (bytetest,1,ierr)
2611 if(ierr.ne.0) goto 102
2612 call mfi_parse_hdr (hdr,ierr) ! replace hdr with correct one
2613 if (nelr.gt.lelr) then
2614 write(6,*) 'ERROR: increase lelr in SIZE to ', nelr
2615 call exitt
2616 endif
2617 call byte_read (er,nelr,ierr) ! get element mapping
2618 if(if_byte_sw) call byte_reverse(er,nelr,ierr)
2619 else
2620 pid0r = 0
2621 pid1r = 0
2622 fid0r = 0
2623 endif
2624
2625 else
2626
2627 pid0r = nid
2628 pid1r = nid
2629 offs0 = iHeaderSize + 4
2630 nfiler = np
2631
2632 ! number of elements to read
2633 nelr = nelgr/np
2634 do i = 0,mod(nelgr,np)-1
2635 if(i.eq.nid) nelr = nelr + 1
2636 enddo
2637 nelBr = igl_running_sum(nelr) - nelr
2638 offs = offs0 + nelBr*isize
2639
2640 call addfid(hname,fid0r)
2641 if(nio.eq.0) write(6,*) ' FILE:',hname
2642 call byte_open_mpi(hname,ifh_mbyte,.true.,ierr)
2643
2644 if(ierr.ne.0) goto 102
2645 call byte_set_view(offs,ifh_mbyte)
2646 if (nelr.gt.lelr) then
2647 write(6,*) 'ERROR: increase lelr in SIZE to ', nelr
2648 call exitt
2649 endif
2650 call byte_read_mpi(er,nelr,-1,ifh_mbyte,ierr)
2651 if(ierr.ne.0) goto 102
2652 if(if_byte_sw) call byte_reverse(er,nelr,ierr)
2653
2654 endif
2655
2656 102 continue
2657 call err_chk(ierr,'Error reading header/element map.$')
2658
2659 return
2660 end
2661c-----------------------------------------------------------------------
2662 subroutine axis_interp_ic(pm1)
2663
2664 include 'SIZE'
2665 include 'TOTAL'
2666 include 'RESTART'
2667
2668 real pm1(lx1,ly1,lz1,lelv)
2669
2670 common /ctmp0/ axism1 (lx1,ly1)
2671 integer e
2672
2673 if (.not.ifaxis) return
2674
2675 do e=1,nelv
2676 if (ifrzer(e)) then
2677 if (ifgetx) then
2678 call mxm (xm1(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
2679 call copy (xm1(1,1,1,e),axism1,lx1*ly1)
2680 call mxm (ym1(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
2681 call copy (ym1(1,1,1,e),axism1,lx1*ly1)
2682 endif
2683 if (ifgetu) then
2684 call mxm (vx(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
2685 call copy (vx(1,1,1,e),axism1,lx1*ly1)
2686 call mxm (vy(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
2687 call copy (vy(1,1,1,e),axism1,lx1*ly1)
2688 endif
2689 if (ifgetw) then
2690 call mxm (vz(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
2691 call copy (vz(1,1,1,e),axism1,lx1*ly1)
2692 endif
2693 if (ifgetp) then
2694 call mxm (pm1(1,1,1,e),lx1,iatlj1,ly1,axism1,ly1)
2695 call copy (pm1(1,1,1,e),axism1,lx1*ly1)
2696 endif
2697 if (ifgett) then
2698 call mxm (t (1,1,1,e,1),lx1,iatlj1,ly1,axism1,ly1)
2699 call copy (t (1,1,1,e,1),axism1,lx1*ly1)
2700 endif
2701 do ips=1,npscal
2702 is1 = ips + 1
2703 if (ifgtps(ips)) then
2704 call mxm (t(1,1,1,e,is1),lx1,iatlj1,ly1,axism1,ly1)
2705 call copy(t(1,1,1,e,is1),axism1,lx1*ly1)
2706 endif
2707 enddo
2708 endif
2709 enddo
2710
2711 return
2712 end
2713c-----------------------------------------------------------------------
2714 subroutine map_pm1_to_pr(pm1,ifile)
2715
2716 include 'SIZE'
2717 include 'TOTAL'
2718 include 'RESTART'
2719
2720 real pm1(lx1*ly1*lz1,lelv)
2721 integer e
2722
2723 nxyz2 = lx2*ly2*lz2
2724
2725 if (ifmhd.and.ifile.eq.2) then
2726 do e=1,nelv
2727 if (if_full_pres) then
2728 call copy (pm(1,1,1,e),pm1(1,e),nxyz2)
2729 else
2730 call map12 (pm(1,1,1,e),pm1(1,e),e)
2731 endif
2732 enddo
2733 elseif (ifsplit) then
2734 call copy (pr,pm1,lx1*ly1*lz1*nelv)
2735 else
2736 do e=1,nelv
2737 if (if_full_pres) then
2738 call copy (pr(1,1,1,e),pm1(1,e),nxyz2)
2739 else
2740 call map12 (pr(1,1,1,e),pm1(1,e),e)
2741 endif
2742 enddo
2743 endif
2744
2745 return
2746 end
2747c-----------------------------------------------------------------------
2748 subroutine full_restart(fnames,n_restart)
2749 include 'SIZE'
2750 include 'TOTAL'
2751
2752 character *(*) fnames(*)
2753
2754 ifile = istep+1 ! istep=0,1,...
2755
2756 if (ifile.le.n_restart) then
2757 p67 = param(67)
2758 param(67) = 6.00
2759 call chcopy (initc,fnames(ifile),80)
2760 call bcast (initc,80)
2761 call restart(1)
2762 call setprop
2763 param(67)=p67
2764 endif
2765
2766 return
2767 end
2768c-----------------------------------------------------------------------
2769 subroutine projfld_c0()
2770
2771 include 'SIZE'
2772 include 'TOTAL'
2773
2774 nxyz1 = lx1*ly1*lz1
2775 ntott = nelt*nxyz1
2776 ntotv = nelv*nxyz1
2777
2778 if(nid.eq.0 .and. loglevel.gt.2) write(6,*) 'projfld_c0'
2779
2780c if (ifflow.and..not.ifdg) then ! Current dg is for scalars only
2781 if (ifflow) then
2782 ifield = 1
2783 call opdssum(vx,vy,vz)
2784 call opcolv (vx,vy,vz,vmult)
2785 if (ifsplit) call dsavg(pr) ! continuous pressure
2786 if (ifvcor) call ortho(pr) ! remove any mean
2787 endif
2788
2789c if (ifmhd.and..not.ifdg) then ! Current dg is for scalars only
2790 if (ifmhd) then
2791 ifield = ifldmhd
2792 call opdssum(bx,by,bz)
2793 call opcolv (bx,by,bz,vmult)
2794 endif
2795
2796 if (ifheat.and..not.ifdg) then ! Don't project if using DG
2797 ifield = 2
2798 call dssum(t ,lx1,ly1,lz1)
2799 call col2 (t ,tmult,ntott)
2800 do ifield=3,nfield
2801 if(gsh_fld(ifield).ge.0) then
2802 call dssum(t(1,1,1,1,ifield-1),lx1,ly1,lz1)
2803 if(iftmsh(ifield)) then
2804 call col2 (t(1,1,1,1,ifield-1),tmult,ntott)
2805 else
2806 call col2 (t(1,1,1,1,ifield-1),vmult,ntotv)
2807 endif
2808 endif
2809 enddo
2810 endif
2811
2812c if (ifpert.and..not.ifdg) then ! Still not DG
2813 if (ifpert) then
2814 do jp=1,npert
2815 ifield = 1
2816 call opdssum(vxp(1,jp),vyp(1,jp),vzp(1,jp))
2817 call opcolv (vxp(1,jp),vyp(1,jp),vzp(1,jp),vmult)
2818
2819 if (.not.ifdg) then
2820 do ifield=2,nfield
2821 call dssum(tp(1,ifield-1,jp),lx1,ly1,lz1)
2822 if(iftmsh(ifield)) then
2823 call col2 (tp(1,ifield-1,jp),tmult,ntott)
2824 else
2825 call col2 (tp(1,ifield-1,jp),vmult,ntotv)
2826 endif
2827 enddo
2828 endif
2829 enddo
2830 endif
2831 jp = 0
2832
2833 return
2834 end
Note: See TracBrowser for help on using the repository browser.