KTH framework for Nek5000 toolboxes; testing version  0.0.1
drive2.f
Go to the documentation of this file.
1  subroutine initdim
2 C-------------------------------------------------------------------
3 C
4 C Transfer array dimensions to common
5 C
6 C-------------------------------------------------------------------
7  include 'SIZE'
8  include 'INPUT'
9 C
10  nelt=lelt
11  nelv=lelv
12 
13  nx1=lx1
14  ny1=ly1
15  nz1=lz1
16 
17  nx2=lx2
18  ny2=ly2
19  nz2=lz2
20 
21  nx3=lx3
22  ny3=ly3
23  nz3=lz3
24 
25  nxd=lxd
26  nyd=lyd
27  nzd=lzd
28 
29  ndim=ldim
30 C
31  RETURN
32  END
33 C
34  subroutine initdat
35 C--------------------------------------------------------------------
36 C
37 C Initialize and set default values.
38 C
39 C--------------------------------------------------------------------
40  include 'SIZE'
41  include 'TOTAL'
42  include 'CTIMER'
43  COMMON /doit/ ifdoit
44  LOGICAL IFDOIT
45 
46 c Set default logicals
47 
48  ifdoit = .false.
49  ifcvode = .false.
50  ifexplvis = .false.
51  ifvvisp = .true.
52 
53  ifsplit = .false.
54  if (lx1.eq.lx2) ifsplit=.true.
55 
56  if_full_pres = .false.
57 
58  CALL rzero (param,200)
59 
60  CALL blank(ccurve ,12*lelt)
61  nel8 = 8*lelt
62  CALL rzero(xc,nel8)
63  CALL rzero(yc,nel8)
64  CALL rzero(zc,nel8)
65 
66  ntot=lx1*ly1*lz1*lelt
67  CALL rzero(abx1,ntot)
68  CALL rzero(abx2,ntot)
69  CALL rzero(aby1,ntot)
70  CALL rzero(aby2,ntot)
71  CALL rzero(abz1,ntot)
72  CALL rzero(abz2,ntot)
73  CALL rzero(vgradt1,ntot)
74  CALL rzero(vgradt2,ntot)
75 
76  ntot=lx2*ly2*lz2*lelt
77  CALL rzero(usrdiv,ntot)
78  CALL rzero(qtl,ntot)
79 
80  nsteps = 0
81 
82  RETURN
83  END
84 C
85  subroutine comment
86 C---------------------------------------------------------------------
87 C
88 C No need to comment !!
89 C
90 C---------------------------------------------------------------------
91  include 'SIZE'
92  include 'INPUT'
93  include 'GEOM'
94  include 'TSTEP'
95  include 'CTIMER'
96 
97  LOGICAL IFCOUR
98  SAVE ifcour
99  COMMON /cprint/ ifprint
100  LOGICAL IFPRINT
101  real*8 eetime0,eetime1,eetime2
102  SAVE eetime0,eetime1,eetime2
103  DATA eetime0,eetime1,eetime2 /0.0, 0.0, 0.0/
104 
105 C
106 C Only node zero makes comments.
107  IF (nio.NE.0) RETURN
108 C
109 C
110  IF (eetime0.EQ.0.0 .AND. istep.EQ.1) eetime0=dnekclock()
111  eetime1=eetime2
112  eetime2=dnekclock()
113 C
114  IF (istep.EQ.0) THEN
115  ifcour = .false.
116  DO 10 ifield=1,nfield
117  IF (ifadvc(ifield)) ifcour = .true.
118  10 CONTINUE
119  IF (ifwcno) ifcour = .true.
120  ELSEIF (istep.GT.0 .AND. lastep.EQ.0 .AND. iftran) THEN
121  ttime_stp = eetime2-eetime1 ! time per timestep
122  ttime = eetime2-eetime0 ! sum of all timesteps
123  IF(istep.EQ.1) THEN
124  ttime_stp = 0
125  ttime = 0
126  ENDIF
127  IF ( ifcour)
128  $ WRITE(6,100)istep,time,dt,courno,ttime,ttime_stp
129  IF (.NOT.ifcour) WRITE (6,101) istep,time,dt
130  ELSEIF (lastep.EQ.1) THEN
131  ttime_stp = eetime2-eetime1 ! time per timestep
132  ttime = eetime2-eetime0 ! sum of all timesteps
133  ENDIF
134  100 FORMAT('Step',i7,', t=',1pe14.7,', DT=',1pe14.7
135  $,', C=',0pf7.3,2(1pe11.4))
136  101 FORMAT('Step',i7,', time=',1pe12.5,', DT=',1pe11.3)
137 
138  RETURN
139  END
140 C
141  subroutine setvar
142 C------------------------------------------------------------------------
143 C
144 C Initialize variables
145 C
146 C------------------------------------------------------------------------
147  include 'SIZE'
148  include 'INPUT'
149  include 'GEOM'
150  include 'DEALIAS'
151  include 'TSTEP'
152  include 'NEKNEK'
153 
154  param(120) = 500 ! print runtime stats
155 
156 C
157 C Geometry on Mesh 3 or 1?
158 C
159  ifgmsh3 = .true.
160  IF ( ifstrs ) ifgmsh3 = .false.
161  IF (.NOT.ifflow) ifgmsh3 = .false.
162  IF ( ifsplit ) ifgmsh3 = .false.
163 
164  ngeom = 2
165 C
166  nfield = 1
167  IF (ifheat) THEN
168  nfield = 2 + npscal
169  nfldtm = 1 + npscal
170  ENDIF
171 c
172  nfldt = nfield
173  if (ifmhd) then
174  nfldt = nfield + 1
175  nfldtm = nfldtm + 1
176  endif
177 c
178  mfield = 1
179  IF (ifmvbd) mfield = 0
180 C
181  DO 100 ifield=mfield,nfldt+(ldimt-1 - npscal)
182  IF (iftmsh(ifield)) THEN
183  nelfld(ifield) = nelt
184  ELSE
185  nelfld(ifield) = nelv
186  ENDIF
187  100 CONTINUE
188 
189  ! Maximum iteration counts for linear solver
190  nmxv = 1000
191  if (iftran) nmxv = 200
192  nmxh = nmxv ! not used anymore
193  nmxp = 200
194  do ifield = 2,ldimt+1
195  nmxt(ifield-1) = 200
196  enddo
197  nmxe = 100
198  nmxnl = 10
199 C
200  param(86) = 0 ! No skew-symm. convection for now
201 C
202  dt = abs(param(12))
203  dtinit = dt
204  fintim = param(10)
205  nsteps = param(11)
206  iocomm = param(13)
207  timeio = param(14)
208  iostep = param(15)
209  lastep = 0
210  tolpdf = abs(param(21))
211  tolhdf = abs(param(22))
212  tolrel = abs(param(24))
213  tolabs = abs(param(25))
214  ctarg = param(26)
215  nbdinp = abs(param(27))
216  nabmsh = param(28)
217 
218  if (nbdinp.gt.lorder) then
219  if (nid.eq.0) then
220  write(6,*) 'ERROR: torder > lorder.',nbdinp,lorder
221  write(6,*) 'Change SIZE and recompile entire code.'
222  endif
223  call exitt
224  endif
225 
226 C Check accuracy requested.
227 C
228  IF (tolrel.LE.0.) tolrel = 0.01
229 C
230 C Relaxed pressure iteration; maximum decrease in the residual.
231 C
232  prelax = 0.1*tolrel
233  IF (.NOT.iftran .AND. .NOT.ifnav) prelax = 1.e-5
234 C
235 C Tolerance for nonlinear iteration
236 C
237  tolnl = 1.e-4
238 C
239 C Fintim overrides nsteps
240 C
241  IF (fintim.NE.0.) nsteps = 1000000000
242  IF (.NOT.iftran ) nsteps = 1
243 C
244 C Print interval defaults to 1
245 C
246  IF (iocomm.EQ.0) iocomm = nsteps+1
247 C
248 C Set default for mesh integration scheme
249 C
250  IF (nabmsh.LE.0 .OR. nabmsh.GT.3) THEN
251  nabmsh = nbdinp
252  param(28) = (nabmsh)
253  ENDIF
254 C
255 C Courant number only applicable if convection in ANY field.
256 C
257  iadv = 0
258  ifld1 = 1
259  IF (.NOT.ifflow) ifld1 = 2
260  DO 200 ifield=ifld1,nfldt
261  IF (ifadvc(ifield)) iadv = 1
262  200 CONTINUE
263 C
264 C If characteristics, need number of sub-timesteps (DT/DS).
265 C Current sub-timeintegration scheme: RK4.
266 C If not characteristics, i.e. standard semi-implicit scheme,
267 C check user-defined Courant number.
268 C
269  IF (iadv.EQ.1) CALL setchar
270 C
271 C Initialize order of time-stepping scheme (BD)
272 C Initialize time step array.
273 C
274  nbd = 0
275  CALL rzero (dtlag,10)
276 
277  ! neknek
278  ifneknekm = .false.
279  ninter = 1
280  nfld_neknek = ndim + nfield
281 
282  CALL blank(cbc_bmap,sizeof(cbc_bmap))
283 
284  one = 1.
285  pi = 4.*atan(one)
286 
287  RETURN
288  END
289 C
290  subroutine echopar
291 C
292 C Echo the nonzero parameters from the readfile to the logfile
293 C
294  include 'SIZE'
295  include 'INPUT'
296  CHARACTER*132 STRING
297  CHARACTER*1 STRING1(132)
298  equivalence(string,string1)
299 C
300  IF (nid.ne.0) RETURN
301 C
302  OPEN (unit=9,file=reafle,status='OLD')
303  rewind(unit=9)
304 C
305 C
306  READ(9,*,err=400)
307  READ(9,*,err=400) vnekton
308  nktonv=vnekton
309  vnekmin=2.5
310  IF(vnekton.LT.vnekmin)THEN
311  print*,' Error: This NEKTON Solver Requires a .rea file'
312  print*,' from prenek version ',vnekmin,' or higher'
313  print*,' Please run the session through the preprocessor'
314  print*,' to bring the .rea file up to date.'
315  call exitt
316  ENDIF
317  READ(9,*,err=400) ldimr
318 c error check
319  IF(ldimr.NE.ldim)THEN
320  WRITE(6,10) ldimr,ldim
321  10 FORMAT(//,2x,'Error: This NEKTON Solver has been compiled'
322  $ /,2x,' for spatial dimension equal to',i2,'.'
323  $ /,2x,' The data file has dimension',i2,'.')
324  CALL exitt
325  ENDIF
326 C
327  CALL blank(string,132)
328 c CALL CHCOPY(STRING,REAFLE,132)
329  ls=ltrunc(string,132)
330  READ(9,*,err=400) nparam
331  WRITE(6,82) nparam,(string1(j),j=1,ls)
332 C
333  DO 20 i=1,nparam
334  CALL blank(string,132)
335  READ(9,80,err=400) string
336  ls=ltrunc(string,132)
337  IF (param(i).ne.0.0) WRITE(6,81) i,(string1(j),j=1,ls)
338  20 CONTINUE
339  80 FORMAT(a132)
340  81 FORMAT(i4,3x,132a1)
341  82 FORMAT(i4,3x,'Parameters from file:',132a1)
342  CLOSE (unit=9)
343  write(6,*) ' '
344 
345 c if(param(2).ne.param(8).and.nio.eq.0) then
346 c write(6,*) 'Note VISCOS not equal to CONDUCT!'
347 c write(6,*) 'Note VISCOS =',PARAM(2)
348 c write(6,*) 'Note CONDUCT =',PARAM(8)
349 c endif
350 c
351  return
352 C
353 C Error handling:
354 C
355  400 CONTINUE
356  WRITE(6,401)
357  401 FORMAT(2x,'ERROR READING PARAMETER DATA'
358  $ ,/,2x,'ABORTING IN ROUTINE ECHOPAR.')
359  CALL exitt
360 C
361  500 CONTINUE
362  WRITE(6,501)
363  501 FORMAT(2x,'ERROR READING LOGICAL DATA'
364  $ ,/,2x,'ABORTING IN ROUTINE ECHOPAR.')
365  CALL exitt
366 C
367  RETURN
368  END
369 C
370  subroutine gengeom (igeom)
371 C----------------------------------------------------------------------
372 C
373 C Generate geometry data
374 C
375 C----------------------------------------------------------------------
376  include 'SIZE'
377  include 'INPUT'
378  include 'TSTEP'
379  include 'GEOM'
380  include 'WZ'
381 C
382  COMMON /scruz/ xm3(lx3,ly3,lz3,lelt)
383  $ , ym3(lx3,ly3,lz3,lelt)
384  $ , zm3(lx3,ly3,lz3,lelt)
385 C
386 
387  if (nio.eq.0.and.istep.le.1) write(6,*) 'generate geometry data'
388 
389  IF (igeom.EQ.1) THEN
390  RETURN
391  ELSEIF (igeom.EQ.2) THEN
392  CALL lagmass
393  IF (istep.EQ.0) CALL gencoor (xm3,ym3,zm3)
394  IF (istep.GE.1) CALL updcoor
395  CALL geom1 (xm3,ym3,zm3)
396  CALL geom2
397  CALL updmsys (1)
398  CALL volume
399  CALL setinvm
400  CALL setdef
401  CALL sfastax
402  ENDIF
403 
404  if (nio.eq.0.and.istep.le.1) then
405  write(6,*) 'done :: generate geometry data'
406  write(6,*) ' '
407  endif
408 
409  return
410  end
411 c-----------------------------------------------------------------------
412  subroutine files
413 C
414 C Defines machine specific input and output file names.
415 C
416  include 'SIZE'
417  include 'INPUT'
418  include 'PARALLEL'
419 C
420  CHARACTER*132 NAME
421  CHARACTER*1 SESS1(132),PATH1(132),NAM1(132)
422  equivalence(session,sess1)
423  equivalence(path,path1)
424  equivalence(name,nam1)
425  CHARACTER*1 DMP(4),FLD(4),REA(4),HIS(4),SCH(4) ,ORE(4), NRE(4)
426  CHARACTER*1 RE2(4),PAR(4)
427  CHARACTER*4 DMP4 ,FLD4 ,REA4 ,HIS4 ,SCH4 ,ORE4 , NRE4
428  CHARACTER*4 RE24 ,PAR4
429  equivalence(dmp,dmp4), (fld,fld4), (rea,rea4), (his,his4)
430  $ , (sch,sch4), (ore,ore4), (nre,nre4)
431  $ , (re2,re24), (par,par4)
432  DATA dmp4,fld4,rea4 /'.dmp','.fld','.rea'/
433  DATA his4,sch4 /'.his','.sch'/
434  DATA ore4,nre4 /'.ore','.nre'/
435  DATA re24 /'.re2' /
436  DATA par4 /'.par' /
437  CHARACTER*78 STRING
438 C
439 C Find out the session name:
440 C
441 c CALL BLANK(SESSION,132)
442 c CALL BLANK(PATH ,132)
443 
444 c ierr = 0
445 c IF(NID.EQ.0) THEN
446 c OPEN (UNIT=8,FILE='SESSION.NAME',STATUS='OLD',ERR=24)
447 c READ(8,10) SESSION
448 c READ(8,10) PATH
449 c 10 FORMAT(A132)
450 c CLOSE(UNIT=8)
451 c GOTO 23
452 c 24 ierr = 1
453 c 23 ENDIF
454 c call err_chk(ierr,' Cannot open SESSION.NAME!$')
455 
456  len = ltrunc(path,132)
457  if(indx1(path1(len),'/',1).lt.1) then
458  call chcopy(path1(len+1),'/',1)
459  endif
460 
461 c call bcast(SESSION,132*CSIZE)
462 c call bcast(PATH,132*CSIZE)
463 
464  CALL blank(parfle,132)
465  CALL blank(reafle,132)
466  CALL blank(re2fle,132)
467  CALL blank(fldfle,132)
468  CALL blank(hisfle,132)
469  CALL blank(schfle,132)
470  CALL blank(dmpfle,132)
471  CALL blank(orefle,132)
472  CALL blank(nrefle,132)
473  CALL blank(name ,132)
474 C
475 C Construct file names containing full path to host:
476 C
477  ls=ltrunc(session,132)
478  lpp=ltrunc(path,132)
479  lsp=ls+lpp
480 c
481  call chcopy(nam1( 1),path1,lpp)
482  call chcopy(nam1(lpp+1),sess1,ls )
483  l1 = lpp+ls+1
484  ln = lpp+ls+4
485 c
486 c
487 c .rea file
488  call chcopy(nam1(l1),rea , 4)
489  call chcopy(reafle ,nam1,ln)
490 c write(6,*) 'reafile:',reafle
491 c
492 c .par file
493  call chcopy(nam1(l1),par , 4)
494  call chcopy(parfle ,nam1,ln)
495 c
496 c .re2 file
497  call chcopy(nam1(l1),re2 , 4)
498  call chcopy(re2fle ,nam1,ln)
499 c
500 c .fld file
501  call chcopy(nam1(l1),fld , 4)
502  call chcopy(fldfle ,nam1,ln)
503 c
504 c .his file
505  call chcopy(nam1(l1),his , 4)
506  call chcopy(hisfle ,nam1,ln)
507 c
508 c .sch file
509  call chcopy(nam1(l1),sch , 4)
510  call chcopy(schfle ,nam1,ln)
511 c
512 c
513 c .dmp file
514  call chcopy(nam1(l1),dmp , 4)
515  call chcopy(dmpfle ,nam1,ln)
516 c
517 c .ore file
518  call chcopy(nam1(l1),ore , 4)
519  call chcopy(orefle ,nam1,ln)
520 c
521 c .nre file
522  call chcopy(nam1(l1),nre , 4)
523  call chcopy(nrefle ,nam1,ln)
524 c
525 C Write the name of the .rea file to the logfile.
526 C
527 C IF (NIO.EQ.0) THEN
528 C CALL CHCOPY(STRING,REAFLE,78)
529 C WRITE(6,1000) STRING
530 C WRITE(6,1001)
531 C 1000 FORMAT(//,1X,'Beginning session:',/,2X,A78)
532 C 1001 FORMAT(/,' ')
533 C ENDIF
534 C
535  RETURN
536 
537  END
538 C
539  subroutine settime
540 C----------------------------------------------------------------------
541 C
542 C Store old time steps and compute new time step, time and timef.
543 C Set time-dependent coefficients in time-stepping schemes.
544 C
545 C----------------------------------------------------------------------
546  include 'SIZE'
547  include 'GEOM'
548  include 'INPUT'
549  include 'TSTEP'
550  COMMON /cprint/ ifprint
551  LOGICAL IFPRINT
552  SAVE
553 C
554  irst = param(46)
555 C
556 C Set time step.
557 C
558  DO 10 ilag=10,2,-1
559  dtlag(ilag) = dtlag(ilag-1)
560  10 CONTINUE
561  CALL setdt
562  dtlag(1) = dt
563  IF (istep.EQ.1 .and. irst.le.0) dtlag(2) = dt
564 C
565 C Set time.
566 C
567  timef = time
568  time = time+dt
569 C
570 C Set coefficients in AB/BD-schemes.
571 C
572  CALL setordbd
573  if (irst.gt.0) nbd = nbdinp
574  CALL rzero (bd,10)
575  CALL setbd (bd,dtlag,nbd)
576  if (param(27).lt.0) then
577  nab = nbdinp
578  else
579  nab = 3
580  endif
581  IF (istep.lt.nab.and.irst.le.0) nab = istep
582  CALL rzero (ab,10)
583  CALL setabbd (ab,dtlag,nab,nbd)
584  IF (ifmvbd) THEN
585  nbdmsh = 1
586  nabmsh = param(28)
587  IF (nabmsh.GT.istep .and. irst.le.0) nabmsh = istep
588  IF (ifsurt) nabmsh = nbd
589  CALL rzero (abmsh,10)
590  CALL setabbd (abmsh,dtlag,nabmsh,nbdmsh)
591  ENDIF
592 
593 C
594 C Set logical for printout to screen/log-file
595 C
596  ifprint = .false.
597  IF (iocomm.GT.0.AND.mod(istep,iocomm).EQ.0) ifprint=.true.
598  IF (istep.eq.1 .or. istep.eq.0 ) ifprint=.true.
599  IF (nio.eq.-1) ifprint=.false.
600 
601  RETURN
602  END
603 
604 
605  subroutine geneig (igeom)
606 C-----------------------------------------------------------------------
607 C
608 C Compute eigenvalues.
609 C Used for automatic setting of tolerances and to find critical
610 C time step for explicit mode.
611 C Currently eigenvalues are computed only for the velocity mesh.
612 C
613 C-----------------------------------------------------------------------
614  include 'SIZE'
615  include 'EIGEN'
616  include 'INPUT'
617  include 'TSTEP'
618 C
619  IF (igeom.EQ.1) RETURN
620 C
621 C Decide which eigenvalues to be computed.
622 C
623  IF (ifflow) THEN
624 C
625  ifaa = .false.
626  ifae = .false.
627  ifas = .false.
628  ifast = .false.
629  ifga = .true.
630  ifge = .false.
631  ifgs = .false.
632  ifgst = .false.
633 C
634 C For now, only compute eigenvalues during initialization.
635 C For deforming geometries the eigenvalues should be
636 C computed every time step (based on old eigenvectors => more memory)
637 C
638  imesh = 1
639  ifield = 1
640  tolev = 1.e-3
641  tolhe = tolhdf
642  tolhr = tolhdf
643  tolhs = tolhdf
644  tolps = tolpdf
645  CALL eigenv
646  CALL esteig
647 C
648  ELSEIF (ifheat.AND..NOT.ifflow) THEN
649 C
650  CALL esteig
651 C
652  ENDIF
653 C
654  RETURN
655  END
656 C-----------------------------------------------------------------------
657  subroutine fluid (igeom)
658 C
659 C Driver for solving the incompressible Navier-Stokes equations.
660 C
661 C Current version:
662 C (1) Velocity/stress formulation.
663 C (2) Constant/variable properties.
664 C (3) Implicit/explicit time stepping.
665 C (4) Automatic setting of tolerances .
666 C (5) Lagrangian/"Eulerian"(operator splitting) modes
667 C
668 C-----------------------------------------------------------------------
669  include 'SIZE'
670  include 'DEALIAS'
671  include 'INPUT'
672  include 'SOLN'
673  include 'TSTEP'
674 
675  real*8 ts, dnekclock
676 
677  ifield = 1
678  imesh = 1
679  call unorm
680  call settolv
681 
682  ts = dnekclock()
683 
684  if(nio.eq.0 .and. igeom.eq.2)
685  & write(*,'(13x,a)') 'Solving for fluid'
686 
687  if (ifsplit) then
688 
689 c PLAN 4: TOMBO SPLITTING
690 c - Time-dependent Navier-Stokes calculation (Re>>1).
691 c - Same approximation spaces for pressure and velocity.
692 c - Incompressibe or Weakly compressible (div u .ne. 0).
693 
694  call plan4 (igeom)
695  if (igeom.ge.2) call chkptol ! check pressure tolerance
696  if (igeom.eq.ngeom) then
697  if (ifneknekc) then
698  call vol_flow_ms ! check for fixed flow rate
699  else
700  call vol_flow ! check for fixed flow rate
701  endif
702  endif
703  if (igeom.ge.2) call printdiverr
704 
705  elseif (iftran) then
706 
707 c call plan1 (igeom) ! Orig. NEKTON time stepper
708 
709  if (ifrich) then
710  call plan5(igeom)
711  else
712  call plan3 (igeom) ! Same as PLAN 1 w/o nested iteration
713  ! Std. NEKTON time stepper !
714  endif
715 
716  if (igeom.ge.2) call chkptol ! check pressure tolerance
717  if (igeom.eq.ngeom) then
718  if (ifneknekc) then
719  call vol_flow_ms ! check for fixed flow rate
720  else
721  call vol_flow ! check for fixed flow rate
722  endif
723  endif
724 
725  else ! steady Stokes, non-split
726 
727 c - Steady/Unsteady Stokes/Navier-Stokes calculation.
728 c - Consistent approximation spaces for velocity and pressure.
729 c - Explicit treatment of the convection term.
730 c - Velocity/stress formulation.
731 
732  call plan1 (igeom) ! The NEKTON "Classic".
733 
734  endif
735 
736  if(nio.eq.0 .and. igeom.ge.2)
737  & write(*,'(4x,i7,a,1p2e12.4)')
738  & istep,' Fluid done',time,dnekclock()-ts
739 
740  return
741  end
742 c-----------------------------------------------------------------------
743  subroutine heat (igeom)
744 C
745 C Driver for temperature or passive scalar.
746 C
747 C Current version:
748 C (1) Varaiable properties.
749 C (2) Implicit time stepping.
750 C (3) User specified tolerance for the Helmholtz solver
751 C (not based on eigenvalues).
752 C (4) A passive scalar can be defined on either the
753 C temperatur or the velocity mesh.
754 C (5) A passive scalar has its own multiplicity (B.C.).
755 C
756  include 'SIZE'
757  include 'INPUT'
758  include 'TSTEP'
759  include 'DEALIAS'
760 
761  real*8 ts, dnekclock
762 
763  ts = dnekclock()
764 
765  if (nio.eq.0 .and. igeom.eq.2)
766  & write(*,'(13x,a)') 'Solving for Hmholtz scalars'
767 
768  do ifield = 2,nfield
769  if (idpss(ifield-1).eq.0) then ! helmholtz
770  intype = -1
771  if (.not.iftmsh(ifield)) imesh = 1
772  if ( iftmsh(ifield)) imesh = 2
773  call unorm
774  call settolt
775  call cdscal(igeom)
776  endif
777  enddo
778 
779  if (nio.eq.0 .and. igeom.eq.2)
780  & write(*,'(4x,i7,a,1p2e12.4)')
781  & istep,' Scalars done',time,dnekclock()-ts
782 
783  return
784  end
785 c-----------------------------------------------------------------------
786  subroutine heat_cvode (igeom)
787 C
788  include 'SIZE'
789  include 'INPUT'
790  include 'TSTEP'
791  include 'DEALIAS'
792 
793  real*8 ts, dnekclock
794 
795  ts = dnekclock()
796 
797  if (igeom.ne.2) return
798 
799  if (nio.eq.0)
800  & write(*,'(13x,a)') 'Solving for CVODE scalars'
801 
802  call cdscal_cvode(igeom)
803 
804  if (nio.eq.0)
805  & write(*,'(4x,i7,a,1p2e12.4)')
806  & istep,' CVODE scalars done',time,dnekclock()-ts
807 
808  return
809  end
810 c-----------------------------------------------------------------------
811  subroutine meshv (igeom)
812 
813 C Driver for mesh velocity used in conjunction with moving geometry.
814 C
815 C-----------------------------------------------------------------------
816  include 'SIZE'
817  include 'INPUT'
818  include 'TSTEP'
819 C
820  IF (igeom.EQ.1) RETURN
821 C
822  ifield = 0
823  nel = nelfld(ifield)
824  imesh = 1
825  IF (iftmsh(ifield)) imesh = 2
826 C
827  CALL updmsys (0)
828  CALL mvbdry (nel)
829  CALL elasolv (nel)
830 C
831  RETURN
832  END
833 c-----------------------------------------------------------------------
834  subroutine time00
835 c
836  include 'SIZE'
837  include 'TOTAL'
838  include 'CTIMER'
839 C
840  nmxmf=0
841  nmxms=0
842  ndsum=0
843  nvdss=0
844  nsett=0
845  ncdtp=0
846  npres=0
847  nmltd=0
848  ngsum=0
849  nprep=0
850  ndsnd=0
851  ndadd=0
852  nhmhz=0
853  naxhm=0
854  ngop =0
855  nusbc=0
856  ncopy=0
857  ninvc=0
858  ninv3=0
859  nsolv=0
860  nslvb=0
861  nddsl=0
862  ncrsl=0
863  ndott=0
864  nbsol=0
865  nadvc=0
866  nspro=0
867  ncvf =0
868 c
869  tmxmf=0.0
870  tmxms=0.0
871  tdsum=0.0
872  tvdss=0.0
873  tvdss=0.0
874  tdsmn=9.9e9
875  tdsmx=0.0
876  tsett=0.0
877  tcdtp=0.0
878  tpres=0.0
879  teslv=0.0
880  tmltd=0.0
881  tgsum=0.0
882  tgsmn=9.9e9
883  tgsmx=0.0
884  tprep=0.0
885  tdsnd=0.0
886  tdadd=0.0
887  thmhz=0.0
888  taxhm=0.0
889  tgop =0.0
890  tusbc=0.0
891  tcopy=0.0
892  tinvc=0.0
893  tinv3=0.0
894  tsolv=0.0
895  tslvb=0.0
896  tddsl=0.0
897  tcrsl=0.0
898  tdott=0.0
899  tbsol=0.0
900  tbso2=0.0
901  tspro=0.0
902  tadvc=0.0
903  ttime=0.0
904  tcvf =0.0
905  tproj=0.0
906  tuchk=0.0
907  tmakf=0.0
908  tmakq=0.0
909 C
910  return
911  end
912 C
913 c-----------------------------------------------------------------------
914  subroutine runstat
915 
916 #ifdef TIMER
917 
918  include 'SIZE'
919  include 'TOTAL'
920  include 'CTIMER'
921 
922  real min_dsum, max_dsum, avg_dsum
923  real min_vdss, max_vdss, avg_vdss
924  real min_gop, max_gop, avg_gop
925  real min_gop_sync, max_gop_sync, avg_gop_sync
926  real min_crsl, max_crsl, avg_crsl
927  real min_usbc, max_usbc, avg_usbc
928  real min_syc, max_syc, avg_syc
929  real min_wal, max_wal, avg_wal
930  real min_irc, max_irc, avg_irc
931  real min_isd, max_isd, avg_isd
932  real min_comm, max_comm, avg_comm
933 
934  real comm_timers(8)
935  integer comm_counters(8)
936  character*132 s132
937 
938  tstop=dnekclock()
939  tttstp=ttime ! sum over all timesteps
940 
941 c call opcount(3) ! print op-counters
942 
943  tcomm = tisd + tirc + tsyc + tgp2+ twal + trc + tsd
944  min_comm = tcomm
945  call gop(min_comm,wwork,'m ',1)
946  max_comm = tcomm
947  call gop(max_comm,wwork,'M ',1)
948  avg_comm = tcomm
949  call gop(avg_comm,wwork,'+ ',1)
950  avg_comm = avg_comm/np
951 c
952  min_isd = tisd
953  call gop(min_isd,wwork,'m ',1)
954  max_isd = tisd
955  call gop(max_isd,wwork,'M ',1)
956  avg_isd = tisd
957  call gop(avg_isd,wwork,'+ ',1)
958  avg_isd = avg_isd/np
959 c
960  min_irc = tirc
961  call gop(min_irc,wwork,'m ',1)
962  max_irc = tirc
963  call gop(max_irc,wwork,'M ',1)
964  avg_irc = tirc
965  call gop(avg_irc,wwork,'+ ',1)
966  avg_irc = avg_irc/np
967 c
968  min_syc = tsyc
969  call gop(min_syc,wwork,'m ',1)
970  max_syc = tsyc
971  call gop(max_syc,wwork,'M ',1)
972  avg_syc = tsyc
973  call gop(avg_syc,wwork,'+ ',1)
974  avg_syc = avg_syc/np
975 c
976  min_wal = twal
977  call gop(min_wal,wwork,'m ',1)
978  max_wal = twal
979  call gop(max_wal,wwork,'M ',1)
980  avg_wal = twal
981  call gop(avg_wal,wwork,'+ ',1)
982  avg_wal = avg_wal/np
983 c
984  min_gop = tgp2
985  call gop(min_gop,wwork,'m ',1)
986  max_gop = tgp2
987  call gop(max_gop,wwork,'M ',1)
988  avg_gop = tgp2
989  call gop(avg_gop,wwork,'+ ',1)
990  avg_gop = avg_gop/np
991 c
992  min_gop_sync = tgop_sync
993  call gop(min_gop_sync,wwork,'m ',1)
994  max_gop_sync = tgop_sync
995  call gop(max_gop_sync,wwork,'M ',1)
996  avg_gop_sync = tgop_sync
997  call gop(avg_gop_sync,wwork,'+ ',1)
998  avg_gop_sync = avg_gop_sync/np
999 c
1000  min_vdss = tvdss
1001  call gop(min_vdss,wwork,'m ',1)
1002  max_vdss = tvdss
1003  call gop(max_vdss,wwork,'M ',1)
1004  avg_vdss = tvdss
1005  call gop(avg_vdss,wwork,'+ ',1)
1006  avg_vdss = avg_vdss/np
1007 c
1008  min_dsum = tdsum
1009  call gop(min_dsum,wwork,'m ',1)
1010  max_dsum = tdsum
1011  call gop(max_dsum,wwork,'M ',1)
1012  avg_dsum = tdsum
1013  call gop(avg_dsum,wwork,'+ ',1)
1014  avg_dsum = avg_dsum/np
1015 c
1016 
1017  min_crsl = tcrsl
1018  call gop(min_crsl,wwork,'m ',1)
1019  max_crsl = tcrsl
1020  call gop(max_crsl,wwork,'M ',1)
1021  avg_crsl = tcrsl
1022  call gop(avg_crsl,wwork,'+ ',1)
1023  avg_crsl = avg_crsl/np
1024 c
1025  min_usbc = tusbc
1026  call gop(min_usbc,wwork,'m ',1)
1027  max_usbc = tusbc
1028  call gop(max_usbc,wwork,'M ',1)
1029  avg_usbc = tusbc
1030  call gop(avg_usbc,wwork,'+ ',1)
1031  avg_usbc = avg_usbc/np
1032 c
1033  tttstp = tttstp + 1e-7
1034  if (nio.eq.0) then
1035  write(6,*) ''
1036  write(6,'(A)') 'runtime statistics:'
1037 
1038  pinit=tinit/tttstp
1039  write(6,*) 'init time',tinit,pinit
1040  pprep=tprep/tttstp
1041  write(6,*) 'prep time',nprep,tprep,pprep
1042 
1043 c Pressure solver timings
1044  ppres=tpres/tttstp
1045  write(6,*) 'pres time',npres,tpres,ppres
1046 
1047 c Coarse grid solver timings
1048  pcrsl=tcrsl/tttstp
1049  write(6,*) 'crsl time',ncrsl,tcrsl,pcrsl
1050  write(6,*) 'crsl min ',min_crsl
1051  write(6,*) 'crsl max ',max_crsl
1052  write(6,*) 'crsl avg ',avg_crsl
1053 
1054 c Helmholz solver timings
1055  phmhz=thmhz/tttstp
1056  write(6,*) 'hmhz time',nhmhz,thmhz,phmhz
1057 
1058 c E solver timings
1059  peslv=teslv/tttstp
1060  write(6,*) 'eslv time',neslv,teslv,peslv
1061 
1062 c makef timings
1063  pmakf=tmakf/tttstp
1064  write(6,*) 'makf time',tmakf,pmakf
1065 
1066 c makeq timings
1067  pmakq=tmakq/tttstp
1068  write(6,*) 'makq time',tmakq,pmakq
1069 
1070 c CVODE RHS timings
1071  pcvf=tcvf/tttstp
1072  if(ifcvode) write(6,*) 'cfun time',ncvf,tcvf,pcvf
1073 
1074 c Resiual projection timings
1075  pproj=tproj/tttstp
1076  write(6,*) 'proj time',tproj,pproj
1077 
1078 c Variable properties timings
1079  pspro=tspro/tttstp
1080  write(6,*) 'usvp time',nspro,tspro,pspro
1081 
1082 c User q and f timings
1083  pusfq=tusfq/tttstp
1084  write(6,*) 'usfq time',0,tusfq,pusfq
1085 
1086 c USERBC timings
1087  pusbc=tusbc/tttstp
1088  write(6,*) 'usbc time',nusbc,tusbc,pusbc
1089  write(6,*) 'usbc min ',min_usbc
1090  write(6,*) 'usbc max ',max_usbc
1091  write(6,*) 'usb avg ',avg_usbc
1092 
1093 c User check timings
1094  puchk=tuchk/tttstp
1095  write(6,*) 'uchk time',tuchk,puchk
1096 
1097 c Operator timings
1098  pmltd=tmltd/tttstp
1099  write(6,*) 'mltd time',nmltd,tmltd,pmltd
1100  pcdtp=tcdtp/tttstp
1101  write(6,*) 'cdtp time',ncdtp,tcdtp,pcdtp
1102  paxhm=taxhm/tttstp
1103  write(6,*) 'axhm time',naxhm,taxhm,paxhm
1104  padvc=tadvc/tttstp
1105  write(6,*) 'advc time',nadvc,tadvc,padvc
1106 
1107 c Low-level routines
1108  pmxmf=tmxmf/tttstp
1109  write(6,*) 'mxmf time',tmxmf,pmxmf
1110  padc3=tadc3/tttstp
1111  write(6,*) 'adc3 time',tadc3,padc3
1112  pcol2=tcol2/tttstp
1113  write(6,*) 'col2 time',tcol2,pcol2
1114  pcol3=tcol3/tttstp
1115  write(6,*) 'col3 time',tcol3,pcol3
1116  pa2s2=ta2s2/tttstp
1117  write(6,*) 'a2s2 time',ta2s2,pa2s2
1118  padd2=tadd2/tttstp
1119  write(6,*) 'add2 time',tadd2,padd2
1120  pinvc=tinvc/tttstp
1121  write(6,*) 'invc time',tinvc,pinvc
1122 
1123 c pinv3=tinv3/tttstp
1124 c write(6,*) 'inv3 time',ninv3,tinv3,pinv3
1125 
1126  pgop=tgop/tttstp
1127  write(6,*) 'tgop time',ngop,tgop,pgop
1128 
1129  pdadd=tdadd/tttstp
1130  write(6,*) 'dadd time',ndadd,tdadd,pdadd
1131 
1132 c Vector direct stiffness summuation timings
1133  pvdss=tvdss/tttstp
1134  write(6,*) 'vdss time',nvdss,tvdss,pvdss
1135  write(6,*) 'vdss min ',min_vdss
1136  write(6,*) 'vdss max ',max_vdss
1137  write(6,*) 'vdss avg ',avg_vdss
1138 
1139 c Direct stiffness summuation timings
1140  pdsum=tdsum/tttstp
1141  write(6,*) 'dsum time',ndsum,tdsum,pdsum
1142  write(6,*) 'dsum min ',min_dsum
1143  write(6,*) 'dsum max ',max_dsum
1144  write(6,*) 'dsum avg ',avg_dsum
1145 
1146 c pgsum=tgsum/tttstp
1147 c write(6,*) 'gsum time',ngsum,tgsum,pgsum
1148 
1149 c pdsnd=tdsnd/tttstp
1150 c write(6,*) 'dsnd time',ndsnd,tdsnd,pdsnd
1151 
1152 c pdsmx=tdsmx/tttstp
1153 c write(6,*) 'dsmx time',ndsmx,tdsmx,pdsmx
1154 c pdsmn=tdsmn/tttstp
1155 c write(6,*) 'dsmn time',ndsmn,tdsmn,pdsmn
1156 c pslvb=tslvb/tttstp
1157 c write(6,*) 'slvb time',nslvb,tslvb,pslvb
1158 
1159  pddsl=tddsl/tttstp
1160  write(6,*) 'ddsl time',nddsl,tddsl,pddsl
1161 
1162 c pbsol=tbsol/tttstp
1163 c write(6,*) 'bsol time',nbsol,tbsol,pbsol
1164 c pbso2=tbso2/tttstp
1165 c write(6,*) 'bso2 time',nbso2,tbso2,pbso2
1166 
1167  write(6,*) ''
1168  endif
1169 
1170  if (lastep.eq.1) then
1171  if (nio.eq.0) ! header for timing
1172  $ write(6,1) 'tusbc','tdadd','tcrsl','tvdss','tdsum',
1173  $ ' tgop',ifsync
1174  1 format(/,'#',2x,'nid',6(7x,a5),4x,'qqq',1x,l4)
1175 
1176  call blank(s132,132)
1177  write(s132,132) nid,tusbc,tdadd,tcrsl,tvdss,tdsum,tgop
1178  132 format(i12,1p6e12.4,' qqq')
1179  call pprint_all(s132,132,6)
1180  endif
1181 #endif
1182 
1183  return
1184  end
1185 c-----------------------------------------------------------------------
1186  subroutine pprint_all(s,n_in,io)
1187  character*1 s(n_in)
1188  character*1 w(132)
1189 
1190  include 'SIZE'
1191  include 'PARALLEL'
1192 
1193  n = min(132,n_in)
1194 
1195  mtag = 999
1196  m = 1
1197  call nekgsync()
1198 
1199  if (nid.eq.0) then
1200  l = ltrunc(s,n)
1201  write(io,1) (s(k),k=1,l)
1202  1 format(132a1)
1203 
1204  do i=1,np-1
1205  call csend(mtag,s,1,i,0) ! send handshake
1206  m = 132
1207  call blank(w,m)
1208  call crecv(i,w,m)
1209  if (m.le.132) then
1210  l = ltrunc(w,m)
1211  write(io,1) (w(k),k=1,l)
1212  else
1213  write(io,*) 'pprint long message: ',i,m
1214  l = ltrunc(w,132)
1215  write(io,1) (w(k),k=1,l)
1216  endif
1217  enddo
1218  else
1219  call crecv(mtag,w,m) ! wait for handshake
1220  l = ltrunc(s,n)
1221  call csend(nid,s,l,0,0) ! send data to node 0
1222  endif
1223  return
1224  end
1225 c-----------------------------------------------------------------------
1226 
1227  subroutine opcount(ICALL)
1228 C
1229  include 'SIZE'
1230  include 'PARALLEL'
1231  include 'OPCTR'
1232 
1233  character*6 sname(maxrts)
1234  integer ind (maxrts)
1235  integer idum (maxrts)
1236 C
1237  if (icall.eq.1) then
1238  nrout=0
1239  endif
1240  if (icall.eq.1.or.icall.eq.2) then
1241  dcount = 0.0
1242  do 100 i=1,maxrts
1243  ncall(i) = 0
1244  dct(i) = 0.0
1245  100 continue
1246  endif
1247  if (icall.eq.3) then
1248 C
1249 C Sort and print out diagnostics
1250 C
1251  if (nid.eq.0) then
1252  write(6,*) nid,' opcount',dcount
1253  do i = 1,np-1
1254  call csend(i,idum,4,i,0)
1255  call crecv(i,ddcount,wdsize)
1256  write(6,*) i,' opcount',ddcount
1257  enddo
1258  else
1259  call crecv (nid,idum,4)
1260  call csend (nid,dcount,wdsize,0,0)
1261  endif
1262 
1263  dhc = dcount
1264  call gop(dhc,dwork,'+ ',1)
1265  if (nio.eq.0) then
1266  write(6,*) ' TOTAL OPCOUNT',dhc
1267  endif
1268 C
1269  CALL drcopy(rct,dct,nrout)
1270  CALL sort(rct,ind,nrout)
1271  CALL chswapr(rname,6,ind,nrout,sname)
1272  call iswap(ncall,ind,nrout,idum)
1273 C
1274  if (nio.eq.0) then
1275  do 200 i=1,nrout
1276  write(6,201) nid,rname(i),rct(i),ncall(i)
1277  200 continue
1278  201 format(2x,' opnode',i4,2x,a6,g18.7,i12)
1279  endif
1280  endif
1281  return
1282  end
1283 C
1284 c-----------------------------------------------------------------------
1285  subroutine dofcnt
1286  include 'SIZE'
1287  include 'TOTAL'
1288  COMMON /scrns/ work(lctmp1)
1289 
1290  integer*8 ntot,ntotp,ntotv
1291 
1292  nxyz = nx1*ny1*nz1
1293  nel = nelv
1294 
1295  ! unique points on v-mesh
1296  vpts = glsum(vmult,nel*nxyz) + .1
1297  nvtot=vpts
1298 
1299  ! unique points on pressure mesh
1300  work(1)=nel*nxyz
1301  ppts = glsum(work,1) + .1
1302  ntot=ppts
1303 C
1304  if (nio.eq.0) write(6,'(A,2i13)')
1305  & 'gridpoints unique/tot: ',nvtot,ntot
1306 
1307  ntot1=nx1*ny1*nz1*nelv
1308  ntot2=nx2*ny2*nz2*nelv
1309 
1310  ntotv = glsc2(tmult,tmask,ntot1)
1311  ntotp = i8glsum(ntot2,1)
1312 
1313  if (ifflow) ntotv = glsc2(vmult,v1mask,ntot1) + .1
1314  if (ifsplit) ntotp = glsc2(vmult,pmask ,ntot1) + .1
1315  if (nio.eq.0) write(6,'(A,2i13)')
1316  $ 'dofs vel/pr: ',ntotv,ntotp
1317 
1318  return
1319  end
1320 c-----------------------------------------------------------------------
1321  subroutine vol_flow
1322 c
1323 c
1324 c Adust flow volume at end of time step to keep flow rate fixed by
1325 c adding an appropriate multiple of the linear solution to the Stokes
1326 c problem arising from a unit forcing in the X-direction. This assumes
1327 c that the flow rate in the X-direction is to be fixed (as opposed to Y-
1328 c or Z-) *and* that the periodic boundary conditions in the X-direction
1329 c occur at the extreme left and right ends of the mesh.
1330 c
1331 c pff 6/28/98
1332 c
1333  include 'SIZE'
1334  include 'TOTAL'
1335 c
1336 c Swap the comments on these two lines if you don't want to fix the
1337 c flow rate for periodic-in-X (or Z) flow problems.
1338 c
1339  parameter(kx1=lx1,ky1=ly1,kz1=lz1,kx2=lx2,ky2=ly2,kz2=lz2)
1340 c
1341  common /cvflow_a/ vxc(kx1,ky1,kz1,lelv)
1342  $ , vyc(kx1,ky1,kz1,lelv)
1343  $ , vzc(kx1,ky1,kz1,lelv)
1344  $ , prc(kx2,ky2,kz2,lelv)
1345  $ , vdc(kx1*ky1*kz1*lelv,2)
1346  common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
1347  $ , scale_vf(3)
1348  common /cvflow_i/ icvflow,iavflow
1349  common /cvflow_c/ chv(3)
1350  character*1 chv
1351 c
1352  real bd_vflow,dt_vflow
1353  save bd_vflow,dt_vflow
1354  data bd_vflow,dt_vflow /-99.,-99./
1355 
1356  logical ifcomp
1357 
1358 c Check list:
1359 
1360 c param (55) -- volume flow rate, if nonzero
1361 c forcing in X? or in Z?
1362 
1363 
1364  ntot1 = lx1*ly1*lz1*nelv
1365  ntot2 = lx2*ly2*lz2*nelv
1366 
1367  if (param(55).eq.0.) return
1368  if (kx1.eq.1) then
1369  write(6,*) 'ABORT. Recompile vol_flow with kx1=lx1, etc.'
1370  call exitt
1371  endif
1372 
1373  icvflow = 1 ! Default flow dir. = X
1374  if (param(54).ne.0) icvflow = abs(param(54))
1375  iavflow = 0 ! Determine flow rate
1376  if (param(54).lt.0) iavflow = 1 ! from mean velocity
1377  flow_rate = param(55)
1378 
1379  chv(1) = 'X'
1380  chv(2) = 'Y'
1381  chv(3) = 'Z'
1382 
1383 c If either dt or the backwards difference coefficient change,
1384 c then recompute base flow solution corresponding to unit forcing:
1385 
1386  ifcomp = .false.
1387  if (dt.ne.dt_vflow.or.bd(1).ne.bd_vflow.or.ifmvbd) ifcomp=.true.
1388  if (.not.ifcomp) then
1389  ifcomp=.true.
1390  do i=1,ntot1
1391  if (vdiff(i,1,1,1,1).ne.vdc(i,1)) goto 20
1392  if (vtrans(i,1,1,1,1).ne.vdc(i,2)) goto 20
1393  enddo
1394  ifcomp=.false. ! If here, then vdiff/vtrans unchanged.
1395  20 continue
1396  endif
1397  call gllog(ifcomp,.true.)
1398 
1399  call copy(vdc(1,1),vdiff(1,1,1,1,1),ntot1)
1400  call copy(vdc(1,2),vtrans(1,1,1,1,1),ntot1)
1401  dt_vflow = dt
1402  bd_vflow = bd(1)
1403 
1404  if (ifcomp) call compute_vol_soln(vxc,vyc,vzc,prc)
1405 
1406  if (icvflow.eq.1) current_flow=glsc2(vx,bm1,ntot1)/domain_length ! for X
1407  if (icvflow.eq.2) current_flow=glsc2(vy,bm1,ntot1)/domain_length ! for Y
1408  if (icvflow.eq.3) current_flow=glsc2(vz,bm1,ntot1)/domain_length ! for Z
1409 
1410  if (iavflow.eq.1) then
1411  xsec = volvm1 / domain_length
1412  flow_rate = param(55)*xsec
1413  endif
1414 
1415  delta_flow = flow_rate-current_flow
1416 
1417 c Note, this scale factor corresponds to FFX, provided FFX has
1418 c not also been specified in userf. If ffx is also specified
1419 c in userf then the true FFX is given by ffx_userf + scale.
1420 
1421  scale = delta_flow/base_flow
1422  scale_vf(icvflow) = scale
1423  if (nio.eq.0) write(6,1) istep,chv(icvflow)
1424  $ ,time,scale,delta_flow,current_flow,flow_rate
1425  1 format(i11,' Volflow ',a1,11x,1p5e13.4)
1426 
1427  call add2s2(vx,vxc,scale,ntot1)
1428  call add2s2(vy,vyc,scale,ntot1)
1429  call add2s2(vz,vzc,scale,ntot1)
1430  call add2s2(pr,prc,scale,ntot2)
1431 
1432  return
1433  end
1434 c-----------------------------------------------------------------------
1435  subroutine compute_vol_soln(vxc,vyc,vzc,prc)
1436 c
1437 c Compute the solution to the time-dependent Stokes problem
1438 c with unit forcing, and find associated flow rate.
1439 c
1440 c pff 2/28/98
1441 c
1442  include 'SIZE'
1443  include 'TOTAL'
1444 c
1445  real vxc(lx1,ly1,lz1,lelv)
1446  $ , vyc(lx1,ly1,lz1,lelv)
1447  $ , vzc(lx1,ly1,lz1,lelv)
1448  $ , prc(lx2,ly2,lz2,lelv)
1449 c
1450  common /cvflow_r/ flow_rate,base_flow,domain_length,xsec
1451  $ , scale_vf(3)
1452  common /cvflow_i/ icvflow,iavflow
1453  common /cvflow_c/ chv(3)
1454  character*1 chv
1455 c
1456  integer icalld
1457  save icalld
1458  data icalld/0/
1459 c
1460 c
1461  ntot1 = lx1*ly1*lz1*nelv
1462  if (icalld.eq.0) then
1463  icalld=icalld+1
1464  xlmin = glmin(xm1,ntot1)
1465  xlmax = glmax(xm1,ntot1)
1466  ylmin = glmin(ym1,ntot1) ! for Y!
1467  ylmax = glmax(ym1,ntot1)
1468  zlmin = glmin(zm1,ntot1) ! for Z!
1469  zlmax = glmax(zm1,ntot1)
1470 c
1471  if (icvflow.eq.1) domain_length = xlmax - xlmin
1472  if (icvflow.eq.2) domain_length = ylmax - ylmin
1473  if (icvflow.eq.3) domain_length = zlmax - zlmin
1474 c
1475  endif
1476 c
1477  if (ifsplit) then
1478 c call plan2_vol(vxc,vyc,vzc,prc)
1479  call plan4_vol(vxc,vyc,vzc,prc)
1480  else
1481  call plan3_vol(vxc,vyc,vzc,prc)
1482  endif
1483 c
1484 c Compute base flow rate
1485 c
1486  if (icvflow.eq.1) base_flow = glsc2(vxc,bm1,ntot1)/domain_length
1487  if (icvflow.eq.2) base_flow = glsc2(vyc,bm1,ntot1)/domain_length
1488  if (icvflow.eq.3) base_flow = glsc2(vzc,bm1,ntot1)/domain_length
1489 c
1490  if (nio.eq.0 .and. loglevel.gt.2) write(6,1)
1491  $ istep,chv(icvflow),base_flow,domain_length,flow_rate
1492  1 format(i11,' basflow ',a1,11x,1p3e13.4)
1493 c
1494  return
1495  end
1496 c-----------------------------------------------------------------------
1497  subroutine plan2_vol(vxc,vyc,vzc,prc)
1498 c
1499 c Compute pressure and velocity using fractional step method.
1500 c (classical splitting scheme).
1501 c
1502 c
1503  include 'SIZE'
1504  include 'TOTAL'
1505 c
1506  real vxc(lx1,ly1,lz1,lelv)
1507  $ , vyc(lx1,ly1,lz1,lelv)
1508  $ , vzc(lx1,ly1,lz1,lelv)
1509  $ , prc(lx2,ly2,lz2,lelv)
1510 C
1511  COMMON /scrns/ resv1(lx1,ly1,lz1,lelv)
1512  $ , resv2(lx1,ly1,lz1,lelv)
1513  $ , resv3(lx1,ly1,lz1,lelv)
1514  $ , respr(lx2,ly2,lz2,lelv)
1515  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
1516  $ , h2(lx1,ly1,lz1,lelv)
1517 c
1518  common /cvflow_i/ icvflow,iavflow
1519 C
1520 C
1521 C Compute pressure
1522 C
1523  ntot1 = lx1*ly1*lz1*nelv
1524 c
1525  if (icvflow.eq.1) then
1526  call cdtp (respr,v1mask,rxm2,sxm2,txm2,1)
1527  elseif (icvflow.eq.2) then
1528  call cdtp (respr,v2mask,rxm2,sxm2,txm2,1)
1529  else
1530  call cdtp (respr,v3mask,rxm2,sxm2,txm2,1)
1531  endif
1532 c
1533  call ortho (respr)
1534 c
1535  call ctolspl (tolspl,respr)
1536  call rone (h1,ntot1)
1537  call rzero (h2,ntot1)
1538 c
1539  call hmholtz ('PRES',prc,respr,h1,h2,pmask,vmult,
1540  $ imesh,tolspl,nmxp,1)
1541  call ortho (prc)
1542 C
1543 C Compute velocity
1544 C
1545  call opgrad (resv1,resv2,resv3,prc)
1546  call opchsgn (resv1,resv2,resv3)
1547  call add2col2 (resv1,bm1,v1mask,ntot1)
1548 c
1549  intype = -1
1550  call sethlm (h1,h2,intype)
1551  call ophinv (vxc,vyc,vzc,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
1552 C
1553  return
1554  end
1555 c-----------------------------------------------------------------------
1556  subroutine plan3_vol(vxc,vyc,vzc,prc)
1557 c
1558 c Compute pressure and velocity using fractional step method.
1559 c (PLAN3).
1560 c
1561 c
1562  include 'SIZE'
1563  include 'TOTAL'
1564 c
1565  real vxc(lx1,ly1,lz1,lelv)
1566  $ , vyc(lx1,ly1,lz1,lelv)
1567  $ , vzc(lx1,ly1,lz1,lelv)
1568  $ , prc(lx2,ly2,lz2,lelv)
1569 C
1570  COMMON /scrns/ rw1(lx1,ly1,lz1,lelv)
1571  $ , rw2(lx1,ly1,lz1,lelv)
1572  $ , rw3(lx1,ly1,lz1,lelv)
1573  $ , dv1(lx1,ly1,lz1,lelv)
1574  $ , dv2(lx1,ly1,lz1,lelv)
1575  $ , dv3(lx1,ly1,lz1,lelv)
1576  $ , respr(lx2,ly2,lz2,lelv)
1577  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
1578  $ , h2(lx1,ly1,lz1,lelv)
1579  COMMON /scrhi/ h2inv(lx1,ly1,lz1,lelv)
1580  common /cvflow_i/ icvflow,iavflow
1581 c
1582 c
1583 c Compute velocity, 1st part
1584 c
1585  ntot1 = lx1*ly1*lz1*nelv
1586  ntot2 = lx2*ly2*lz2*nelv
1587  ifield = 1
1588 c
1589  if (icvflow.eq.1) then
1590  call copy (rw1,bm1,ntot1)
1591  call rzero (rw2,ntot1)
1592  call rzero (rw3,ntot1)
1593  elseif (icvflow.eq.2) then
1594  call rzero (rw1,ntot1)
1595  call copy (rw2,bm1,ntot1)
1596  call rzero (rw3,ntot1)
1597  else
1598  call rzero (rw1,ntot1) ! Z-flow!
1599  call rzero (rw2,ntot1) ! Z-flow!
1600  call copy (rw3,bm1,ntot1) ! Z-flow!
1601  endif
1602  intype = -1
1603  call sethlm (h1,h2,intype)
1604  call ophinv (vxc,vyc,vzc,rw1,rw2,rw3,h1,h2,tolhv,nmxv)
1605  call ssnormd (vxc,vyc,vzc)
1606 c
1607 c Compute pressure (from "incompr")
1608 c
1609  intype = 1
1610  dtinv = 1./dt
1611 c
1612  call rzero (h1,ntot1)
1613  call copy (h2,vtrans(1,1,1,1,ifield),ntot1)
1614  call cmult (h2,dtinv,ntot1)
1615  call invers2 (h2inv,h2,ntot1)
1616  call opdiv (respr,vxc,vyc,vzc)
1617  call chsign (respr,ntot2)
1618  call ortho (respr)
1619 c
1620 c
1621 c Set istep=0 so that h1/h2 will be re-initialized in eprec
1622  i_tmp = istep
1623  istep = 0
1624  call esolver (respr,h1,h2,h2inv,intype)
1625  istep = i_tmp
1626 c
1627  call opgradt (rw1,rw2,rw3,respr)
1628  call opbinv (dv1,dv2,dv3,rw1,rw2,rw3,h2inv)
1629  call opadd2 (vxc,vyc,vzc,dv1,dv2,dv3)
1630 c
1631  call cmult2 (prc,respr,bd(1),ntot2)
1632 c
1633  return
1634  end
1635 c-----------------------------------------------------------------------
1636  subroutine plan4_vol(vxc,vyc,vzc,prc)
1637 
1638 c Compute pressure and velocity using fractional step method.
1639 c (Tombo splitting scheme).
1640 
1641 
1642 
1643  include 'SIZE'
1644  include 'TOTAL'
1645 
1646  real vxc(lx1,ly1,lz1,lelv)
1647  $ , vyc(lx1,ly1,lz1,lelv)
1648  $ , vzc(lx1,ly1,lz1,lelv)
1649  $ , prc(lx1,ly1,lz1,lelv)
1650 
1651  common /scrns/ resv1(lx1,ly1,lz1,lelv)
1652  $ , resv2(lx1,ly1,lz1,lelv)
1653  $ , resv3(lx1,ly1,lz1,lelv)
1654  $ , respr(lx1,ly1,lz1,lelv)
1655  common /scrvh/ h1(lx1,ly1,lz1,lelv)
1656  $ , h2(lx1,ly1,lz1,lelv)
1657 
1658  common /cvflow_i/ icvflow,iavflow
1659 
1660  n = lx1*ly1*lz1*nelv
1661  call invers2 (h1,vtrans,n)
1662  call rzero (h2, n)
1663 
1664 c Compute pressure
1665 
1666  if (icvflow.eq.1) call cdtp(respr,h1,rxm2,sxm2,txm2,1)
1667  if (icvflow.eq.2) call cdtp(respr,h1,rym2,sym2,tym2,1)
1668  if (icvflow.eq.3) call cdtp(respr,h1,rzm2,szm2,tzm2,1)
1669 
1670  call ortho (respr)
1671  call ctolspl (tolspl,respr)
1672 
1673  call hmholtz ('PRES',prc,respr,h1,h2,pmask,vmult,
1674  $ imesh,tolspl,nmxp,1)
1675  call ortho (prc)
1676 
1677 C Compute velocity
1678 
1679  call opgrad (resv1,resv2,resv3,prc)
1680  if (ifaxis) call col2 (resv2,omask,n)
1681  call opchsgn (resv1,resv2,resv3)
1682 
1683  if (icvflow.eq.1) call add2col2(resv1,v1mask,bm1,n) ! add forcing
1684  if (icvflow.eq.2) call add2col2(resv2,v2mask,bm1,n)
1685  if (icvflow.eq.3) call add2col2(resv3,v3mask,bm1,n)
1686 
1687 
1688  if (ifexplvis) call split_vis ! split viscosity into exp/imp part
1689 
1690  intype = -1
1691  call sethlm (h1,h2,intype)
1692  call ophinv (vxc,vyc,vzc,resv1,resv2,resv3,h1,h2,tolhv,nmxv)
1693 
1694  if (ifexplvis) call redo_split_vis ! restore vdiff
1695 
1696  end
1697 c-----------------------------------------------------------------------
1698  subroutine a_dmp
1699 c
1700  include 'SIZE'
1701  include 'TOTAL'
1702  COMMON /scrns/ w(lx1,ly1,lz1,lelt)
1703  COMMON /scruz/ v(lx1,ly1,lz1,lelt)
1704  $ , h1(lx1,ly1,lz1,lelt)
1705  $ , h2(lx1,ly1,lz1,lelt)
1706 c
1707  ntot = lx1*ly1*lz1*nelv
1708  call rone (h1,ntot)
1709  call rzero(h2,ntot)
1710  do i=1,ntot
1711  call rzero(v,ntot)
1712  v(i,1,1,1) = 1.
1713  call axhelm (w,v,h1,h2,1,1)
1714  call outrio (w,ntot,55)
1715  enddo
1716 c write(6,*) 'quit in a_dmp'
1717 c call exitt
1718  return
1719  end
1720 c-----------------------------------------------------------------------
1721  subroutine outrio (v,n,io)
1722 c
1723  real v(1)
1724 c
1725  write(6,*) 'outrio:',n,io,v(1)
1726  write(io,6) (v(k),k=1,n)
1727  6 format(1pe19.11)
1728 c
1729 c nr = min(12,n)
1730 c write(io,6) (v(k),k=1,nr)
1731 c 6 format(1p12e11.3)
1732  return
1733  end
1734 c-----------------------------------------------------------------------
1735  subroutine reset_prop
1736 C------------------------------------------------------------------------
1737 C
1738 C Set variable property arrays
1739 C
1740 C------------------------------------------------------------------------
1741  include 'SIZE'
1742  include 'TOTAL'
1743 C
1744 C Caution: 2nd and 3rd strainrate invariants residing in scratch
1745 C common /SCREV/ are used in STNRINV and NEKASGN
1746 C
1747  COMMON /screv/ sii(lx1,ly1,lz1,lelt)
1748  $ , siii(lx1,ly1,lz1,lelt)
1749  COMMON /scruz/ ta(lx1,ly1,lz1,lelt)
1750 C
1751  real rstart
1752  save rstart
1753  data rstart /1/
1754 c
1755  rfinal = 1./param(2) ! Target Re
1756 c
1757  ntot = lx1*ly1*lz1*nelv
1758  iramp = 200
1759  istpp = istep
1760 c istpp = istep+2033+1250
1761  if (istpp.ge.iramp) then
1762  vfinal=1./rfinal
1763  call cfill(vdiff,vfinal,ntot)
1764  else
1765  one = 1.
1766  pi2 = 2.*atan(one)
1767  sarg = (pi2*istpp)/iramp
1768  sarg = sin(sarg)
1769  rnew = rstart + (rfinal-rstart)*sarg
1770  vnew = 1./rnew
1771  call cfill(vdiff,vnew,ntot)
1772  if (nio.eq.0) write(6,*) istep,' New Re:',rnew,sarg,istpp
1773  endif
1774  return
1775  end
1776 C-----------------------------------------------------------------------
1777  subroutine prinit
1778 
1779  include 'SIZE'
1780  include 'TOTAL'
1781 
1782  if(nio.eq.0) write(6,*) 'initialize pressure solver'
1783  isolver = param(40)
1784 
1785  if (isolver.eq.0) then ! semg_xxt
1786  if (nelgt.gt.350000)
1787  $ call exitti('problem size too large for XXT solver!$',0)
1788  call set_overlap
1789  else if (isolver.eq.1) then ! semg_amg
1790  call set_overlap
1791  else if (isolver.eq.2) then ! semg_amg_hypre
1792  call set_overlap
1793  else if (isolver.eq.3) then ! fem_amg_hypre
1794  null_space = 0
1795  if (ifvcor) null_space = 1
1796  call fem_amg_setup(nx1,ny1,nz1,
1797  $ nelv,ndim,
1798  $ xm1,ym1,zm1,
1799  $ pmask,binvm1,null_space,
1800  $ gsh_fld(1),fem_amg_param)
1801  endif
1802 
1803  return
1804  end
exitt
void exitt()
Definition: comm_mpi.f:604
setinvm
subroutine setinvm
Definition: coef.f:1334
geom1
subroutine geom1(xm3, ym3, zm3)
Definition: coef.f:361
volume
subroutine volume
Definition: coef.f:1003
geom2
subroutine geom2
Definition: coef.f:823
lagmass
subroutine lagmass
Definition: coef.f:1315
crecv
subroutine crecv(mtype, buf, lenm)
Definition: comm_mpi.f:313
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
csend
subroutine csend(mtype, buf, len, jnid, jpid)
Definition: comm_mpi.f:303
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
cdscal
subroutine cdscal(igeom)
Definition: conduct.f:3
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
cdscal_cvode
subroutine cdscal_cvode
Definition: cvode_driver.f:21
time00
subroutine time00
Definition: drive2.f:835
files
subroutine files
Definition: drive2.f:413
initdim
subroutine initdim
Definition: drive2.f:2
compute_vol_soln
subroutine compute_vol_soln(vxc, vyc, vzc, prc)
Definition: drive2.f:1436
pprint_all
subroutine pprint_all(s, n_in, io)
Definition: drive2.f:1187
runstat
subroutine runstat
Definition: drive2.f:915
settime
subroutine settime
Definition: drive2.f:540
prinit
subroutine prinit
Definition: drive2.f:1778
geneig
subroutine geneig(igeom)
Definition: drive2.f:606
outrio
subroutine outrio(v, n, io)
Definition: drive2.f:1722
meshv
subroutine meshv(igeom)
Definition: drive2.f:812
fluid
subroutine fluid(igeom)
Definition: drive2.f:658
plan2_vol
subroutine plan2_vol(vxc, vyc, vzc, prc)
Definition: drive2.f:1498
comment
subroutine comment
Definition: drive2.f:86
reset_prop
subroutine reset_prop
Definition: drive2.f:1736
setvar
subroutine setvar
Definition: drive2.f:142
heat
subroutine heat(igeom)
Definition: drive2.f:744
initdat
subroutine initdat
Definition: drive2.f:35
echopar
subroutine echopar
Definition: drive2.f:291
dofcnt
subroutine dofcnt
Definition: drive2.f:1286
opcount
subroutine opcount(ICALL)
Definition: drive2.f:1228
vol_flow
subroutine vol_flow
Definition: drive2.f:1322
heat_cvode
subroutine heat_cvode(igeom)
Definition: drive2.f:787
gengeom
subroutine gengeom(igeom)
Definition: drive2.f:371
plan3_vol
subroutine plan3_vol(vxc, vyc, vzc, prc)
Definition: drive2.f:1557
a_dmp
subroutine a_dmp
Definition: drive2.f:1699
plan4_vol
subroutine plan4_vol(vxc, vyc, vzc, prc)
Definition: drive2.f:1637
esteig
subroutine esteig
Definition: eigsolv.f:9
eigenv
subroutine eigenv
Definition: eigsolv.f:91
fem_amg_setup
#define fem_amg_setup
Definition: fem_amg_preco.c:8
setdef
subroutine setdef
Definition: genxyz.f:352
gencoor
subroutine gencoor(xm3, ym3, zm3)
Definition: genxyz.f:534
sfastax
subroutine sfastax
Definition: hmholtz.f:311
hmholtz
subroutine hmholtz(name, u, rhs, h1, h2, mask, mult, imsh, tli, maxit, isd)
Definition: hmholtz.f:3
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
indx1
integer function indx1(S1, S2, L2)
Definition: ic.f:1357
ophinv
subroutine ophinv(o1, o2, o3, i1, i2, i3, h1, h2, tolh, nmxhi)
Definition: induct.f:1023
iswap
subroutine iswap(b, ind, n, temp)
Definition: math.f:549
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
glmin
function glmin(a, n)
Definition: math.f:973
add2col2
subroutine add2col2(a, b, c, n)
Definition: math.f:1565
glsc2
function glsc2(x, y, n)
Definition: math.f:794
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
copy
subroutine copy(a, b, n)
Definition: math.f:260
blank
subroutine blank(A, N)
Definition: math.f:19
i8glsum
integer *8 function i8glsum(a, n)
Definition: math.f:947
chswapr
subroutine chswapr(b, L, ind, n, temp)
Definition: math.f:1210
ltrunc
function ltrunc(string, l)
Definition: math.f:494
gllog
subroutine gllog(la, lb)
Definition: math.f:986
glsum
function glsum(x, n)
Definition: math.f:861
drcopy
subroutine drcopy(r, d, N)
Definition: math.f:1226
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rone
subroutine rone(a, n)
Definition: math.f:230
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
chsign
subroutine chsign(a, n)
Definition: math.f:305
sort
subroutine sort(a, ind, n)
Definition: math.f:1325
vol_flow_ms
subroutine vol_flow_ms
Definition: multimesh.f:655
mvbdry
subroutine mvbdry(nel)
Definition: mvmesh.f:250
updcoor
subroutine updcoor
Definition: mvmesh.f:226
elasolv
subroutine elasolv(nel)
Definition: mvmesh.f:729
esolver
subroutine esolver(RES, H1, H2, H2INV, INTYPE)
Definition: navier0.f:2
opgradt
subroutine opgradt(outx, outy, outz, inpfld)
Definition: navier1.f:4096
cdtp
subroutine cdtp(dtx, x, rm2, sm2, tm2, isd)
Definition: navier1.f:331
plan1
subroutine plan1(igeom)
Definition: navier1.f:2
ortho
subroutine ortho(respr)
Definition: navier1.f:224
opbinv
subroutine opbinv(out1, out2, out3, inp1, inp2, inp3, h2inv)
Definition: navier1.f:776
opdiv
subroutine opdiv(outfld, inpx, inpy, inpz)
Definition: navier1.f:4065
setordbd
subroutine setordbd
Definition: navier1.f:2003
setbd
subroutine setbd(bd, dtbd, nbd)
Definition: navier1.f:1736
setabbd
subroutine setabbd(ab, dtlag, nab, nbd)
Definition: navier1.f:1663
ctolspl
subroutine ctolspl(tolspl, respr)
Definition: navier1.f:194
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
opchsgn
subroutine opchsgn(a, b, c)
Definition: navier1.f:2464
opgrad
subroutine opgrad(out1, out2, out3, inp)
Definition: navier1.f:298
chkptol
subroutine chkptol
Definition: navier4.f:269
set_overlap
subroutine set_overlap
Definition: navier6.f:31
plan4
subroutine plan4(igeom)
Definition: plan4.f:3
split_vis
subroutine split_vis
Definition: plan4.f:459
redo_split_vis
subroutine redo_split_vis
Definition: plan4.f:477
printdiverr
subroutine printdiverr
Definition: plan4.f:652
plan5
subroutine plan5(igeom)
Definition: plan5.f:3
plan3
subroutine plan3(IGEOM)
Definition: planx.f:2
ssnormd
subroutine ssnormd(DV1, DV2, DV3)
Definition: ssolv.f:629
settolv
subroutine settolv
Definition: ssolv.f:648
settolt
subroutine settolt
Definition: ssolv.f:693
setchar
subroutine setchar
Definition: ssolv.f:744
setdt
subroutine setdt
Definition: subs1.f:179
unorm
subroutine unorm
Definition: subs1.f:352
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341
updmsys
subroutine updmsys(IFLD)
Definition: subs2.f:712