KTH framework for Nek5000 toolboxes; testing version  0.0.1
ssolv.f
Go to the documentation of this file.
1  SUBROUTINE sstest (ISSS)
2 C------------------------------------------------------------------------
3 C
4 C Test if Steady State Solver should be activated.
5 C
6 C------------------------------------------------------------------------
7  include 'SIZE'
8  include 'INPUT'
9  include 'TSTEP'
10  isss = 0
11  iadv = 0
12  DO 100 ifield=1,nfield
13  IF (ifadvc(ifield)) iadv = 1
14  100 CONTINUE
15  IF (.NOT.iftran.AND.(iadv.EQ.1)) isss = 1
16  IF (isss.EQ.1.AND.nfield.GT.4.AND.nid.EQ.0) THEN
17  WRITE (6,*) ' '
18  WRITE (6,*) 'Trying to activate the steady state solver'
19  WRITE (6,*) 'using NFIELD =',nfield
20  WRITE (6,*) 'Maximum number of fields is 4'
21  call exitt
22  ENDIF
23  RETURN
24  END
25 C
26  SUBROUTINE ssinit (KMAX)
27 C-----------------------------------------------------------
28 C
29 C Initialize steady state solver
30 C
31 C-----------------------------------------------------------
32  include 'SIZE'
33  include 'INPUT'
34  include 'EIGEN'
35  include 'TSTEP'
36  include 'STEADY'
37 C
38  iftran = .true.
39  ifchar = .true.
40  ctarg = 3.
41  IF (lx1.GE.10) THEN
42  ctarg = 5.
43  ENDIF
44  CALL setprop
45  taumin = 1.e20
46  mfield = 1
47  IF (.NOT.ifflow) mfield=2
48  DO 10 ifield=mfield,nfield
49  diffus = avdiff(ifield)/avtran(ifield)
50  tauss(ifield) = 1./(eigaa*diffus)
51  txnext(ifield) = tauss(ifield)
52  IF (tauss(ifield).LT.taumin) taumin = tauss(ifield)
53  10 CONTINUE
54 C
55  nbdinp = 1.
56  time = 0.
57  dt = 0.
58  dtinit = taumin/5.
59  nsteps = 10000
60  iostep = 10000
61 C
62  ifmodp = .true.
63  ifskip = .true.
64  nsskip = 1
65 C
66  prelax = 1.e-1
67  IF (ifsplit) prelax = 1.e-4
68 C
69  ifssvt = .false.
70  ifexvt = .false.
71 C
72  kmax = 5
73 C
74  CALL setchar
75 C
76 C
77  RETURN
78  END
79 C
80  SUBROUTINE chkext (IFACCX,Z,S)
81 C------------------------------------------------------------------
82 C
83 C Accept extrapolation?
84 C
85 C------------------------------------------------------------------
86  include 'SIZE'
87  include 'TSTEP'
88  include 'INPUT'
89  include 'STEADY'
90  LOGICAL IFACCX
91  real*8 z(1),s(1)
92  REAL H1NRM1 (LDIMT1), H1NRM2(LDIMT1)
93 C
94  CALL rzero (h1nrm1,nfield)
95  IF (ifflow) THEN
96  ifield = 1
97  CALL unorm
98  h1nrm1(ifield) = vnrmh1
99  ENDIF
100  DO 10 ifield=2,nfield
101  CALL unorm
102  h1nrm1(ifield) = tnrmh1(ifield-1)
103  10 CONTINUE
104 C
105  CALL mkvec (z)
106  CALL mkarr (s)
107 C
108  CALL rzero (h1nrm2,nfield)
109  IF (ifflow) THEN
110  ifield = 1
111  CALL unorm
112  h1nrm2(ifield) = vnrmh1
113  ENDIF
114  DO 20 ifield=2,nfield
115  CALL unorm
116  h1nrm2(ifield) = tnrmh1(ifield-1)
117  20 CONTINUE
118 C
119  xlim = .2
120  ifaccx = .true.
121  rdmax = 0.
122  rdlim = .5*tolrel+1.e-4
123  IF (ifflow) THEN
124  ifield = 1
125  rdiff = abs((h1nrm2(ifield)-h1nrm1(ifield))/h1nrm1(ifield))
126  IF (rdiff.GT.rdmax) rdmax = rdiff
127  IF (nio.EQ.0) WRITE (6,*) ' ifield, rdiff ',ifield,rdiff
128  IF (rdiff.GT.xlim) ifaccx = .false.
129  ENDIF
130  DO 100 ifield=2,nfield
131  rdiff = abs((h1nrm2(ifield)-h1nrm1(ifield))/h1nrm1(ifield))
132  IF (rdiff.GT.rdmax) rdmax = rdiff
133  IF (nio.EQ.0) WRITE (6,*) ' ifield, rdiff ',ifield,rdiff
134  IF (rdiff.GT.xlim) ifaccx = .false.
135  100 CONTINUE
136 C
137  IF (.NOT.ifaccx) THEN
138  IF (nio.EQ.0) THEN
139  WRITE (6,*) ' '
140  write (6,*) 'Extrapolation attempt rejected'
141  write (6,*) ' '
142  ENDIF
143  CALL mkarr (z)
144  ELSE
145  IF (nio.EQ.0) THEN
146  write (6,*) ' '
147  write (6,*) 'Extrapolation accepted'
148  write (6,*) ' '
149  ENDIF
150  IF (rdmax.LT.rdlim) ifssvt = .true.
151  CALL filllag
152  ENDIF
153 C
154  RETURN
155  END
156 C
157  SUBROUTINE filllag
158  include 'SIZE'
159  include 'SOLN'
160  include 'INPUT'
161  include 'TSTEP'
162  nbdinp = 3
163  IF (ifflow) THEN
164  CALL lagvel
165  CALL lagvel
166  ENDIF
167  IF (ifheat) THEN
168  DO 100 ifield=2,nfield
169  CALL lagscal
170  CALL lagscal
171  100 CONTINUE
172  ENDIF
173  RETURN
174  END
175 C
176  SUBROUTINE gonstep (N,ITEST)
177 C----------------------------------------------------------------
178 C
179 C Do N steps; return if steady state
180 C
181 C----------------------------------------------------------------
182  include 'SIZE'
183  include 'INPUT'
184  include 'STEADY'
185  EXTERNAL gostep
186 C
187  DO 1000 jstep=1,n
188  IF (itest.EQ.0.AND.ifssvt) GOTO 1001
189  IF (itest.EQ.1.AND.(ifssvt.OR.ifexvt)) GOTO 1001
190  CALL gostep
191  1000 CONTINUE
192  1001 CONTINUE
193 C
194  RETURN
195  END
196 C
197  SUBROUTINE go1step (X,Y,NVEC)
198 C----------------------------------------------------------------
199 C
200 C Advance one (or more) time step(s)
201 C
202 C----------------------------------------------------------------
203  include 'SIZE'
204  include 'TOTAL'
205  real*8 x(1), y(1)
206 C
207  CALL mkarr (x)
208  IF (.NOT.ifskip) njstep=1
209  IF ( ifskip) njstep=nsskip
210 C
211  DO 9000 jstep=1,njstep
212 C
213  istep = istep+1
214  CALL settime
215  CALL setprop
216  IF (ifmodp) CALL modprop
217  CALL setsolv
218  CALL comment
219  DO 100 igeom=1,2
220  IF (ifgeom) THEN
221  CALL gengeom (igeom)
222  CALL geneig (igeom)
223  ENDIF
224  IF (ifflow) CALL fluid (igeom)
225  IF (ifheat) CALL heat (igeom)
226  IF (ifmvbd) CALL meshv (igeom)
227  100 CONTINUE
228  CALL prepost(.false.)
229  CALL userchk
230 C
231  9000 CONTINUE
232 C
233  IF (istep.GT.1) CALL chkssvt
234  CALL mkvec (y)
235 C
236  RETURN
237  END
238 C
239  SUBROUTINE gostep
240 C----------------------------------------------------------------
241 C
242 C Advance one (or more) time step(s)
243 C
244 C----------------------------------------------------------------
245  include 'SIZE'
246  include 'TOTAL'
247 C
248  IF (.NOT.ifskip) njstep=1
249  IF ( ifskip) njstep=nsskip
250 C
251  DO 9000 jstep=1,njstep
252 C
253  istep = istep+1
254  CALL settime
255  CALL setprop
256  IF (ifmodp) CALL modprop
257  CALL setsolv
258  CALL comment
259  DO 100 igeom=1,2
260  IF (ifgeom) THEN
261  CALL gengeom (igeom)
262  CALL geneig (igeom)
263  ENDIF
264  IF (ifflow) CALL fluid (igeom)
265  IF (ifheat) CALL heat (igeom)
266  IF (ifmvbd) CALL meshv (igeom)
267  100 CONTINUE
268  CALL prepost(.false.)
269  CALL userchk
270 C
271  9000 CONTINUE
272 C
273  IF (istep.GT.1) CALL chkssvt
274 C
275  RETURN
276  END
277 C
278  SUBROUTINE modprop
279 C------------------------------------------------------------------
280 C
281 C Modify the properties
282 C
283 C------------------------------------------------------------------
284  include 'SIZE'
285  include 'SOLN'
286  include 'TSTEP'
287  include 'STEADY'
288  include 'INPUT'
289 C
290  mfield=1
291  IF (.NOT.ifflow) mfield=2
292  DO 100 ifield=mfield,nfield
293  ntot = lx1*ly1*lz1*nelfld(ifield)
294  tau = .02*tauss(ifield)
295  decay = 1.+99.*exp(-time/tau)
296  CALL cmult (vdiff(1,1,1,1,ifield),decay,ntot)
297 c if (nid.eq.0)
298 c $ write (6,*) '.......... diff = ',IFIELD,vdiff(1,1,1,1,IFIELD)
299  100 CONTINUE
300 C
301  RETURN
302  END
303 C
304  SUBROUTINE mkvec (X)
305 C-------------------------------------------------------------
306 C
307 C Fill up the vector X with VX, VY, ....
308 C
309 C-------------------------------------------------------------
310  include 'SIZE'
311  include 'SOLN'
312  include 'INPUT'
313  include 'TSTEP'
314  real*8 x(1)
315 C
316  ntotv = lx1*ly1*lz1*nelv
317 C
318  IF (ifflow) THEN
319  DO 100 i=1,ntotv
320  x(i) = vx(i,1,1,1)
321  100 CONTINUE
322  DO 200 i=1,ntotv
323  x(i+ntotv) = vy(i,1,1,1)
324  200 CONTINUE
325  IF (ldim.EQ.3) THEN
326  ioff = 2*ntotv
327  DO 300 i=1,ntotv
328  x(i+ioff) = vz(i,1,1,1)
329  300 CONTINUE
330  ENDIF
331  ENDIF
332 C
333  IF (ifheat) THEN
334  ioff = ldim*ntotv
335  DO 401 ifield=2,nfield
336  ntot = lx1*ly1*lz1*nelfld(ifield)
337  DO 400 i=1,ntot
338  x(i+ioff) = t(i,1,1,1,ifield-1)
339  400 CONTINUE
340  ioff = ioff+ntot
341  401 CONTINUE
342  ENDIF
343 C
344  RETURN
345  END
346 C
347  SUBROUTINE mkarr (X)
348 C------------------------------------------------------------------
349 C
350 C Split the vector X into VX, VY, .....
351 C
352 C------------------------------------------------------------------
353  include 'SIZE'
354  include 'SOLN'
355  include 'TSTEP'
356  include 'INPUT'
357  real*8 x(1)
358 C
359  ntotv = lx1*ly1*lz1*nelv
360 C
361  IF (ifflow) THEN
362  DO 10 i=1,ntotv
363  vx(i,1,1,1) = x(i)
364  10 CONTINUE
365  DO 20 i=1,ntotv
366  vy(i,1,1,1) = x(i+ntotv)
367  20 CONTINUE
368  IF (ldim.EQ.3) THEN
369  ioff = 2*ntotv
370  DO 30 i=1,ntotv
371  vz(i,1,1,1) = x(i+ioff)
372  30 CONTINUE
373  ENDIF
374  ENDIF
375 C
376  IF (ifheat) THEN
377  ioff = ldim*ntotv
378  DO 41 ifield=2,nfield
379  ntot = lx1*ly1*lz1*nelfld(ifield)
380  DO 40 i=1,ntot
381  t(i,1,1,1,ifield-1) = x(i+ioff)
382  40 CONTINUE
383  ioff = ioff+ntot
384  41 CONTINUE
385  ENDIF
386 C
387  RETURN
388  END
389 C
390  SUBROUTINE ssparam (KMAX,L)
391 C------------------------------------------------------------------------------
392 C
393 C Set steady state parameters
394 C
395 C------------------------------------------------------------------------------
396  include 'SIZE'
397  include 'TSTEP'
398  include 'INPUT'
399  include 'STEADY'
400 C
401  IF (l.EQ.0) THEN
402  CALL ssinit (kmax)
403  ELSEIF (l.EQ.1) THEN
404  istep = 0
405 C
406  prelax = 1.e-2
407  IF (ifsplit) prelax = 1.e-5
408 C
409  ctarg = 1.
410  IF (ifsplit) ctarg = 1.
411  IF (lx1.GE.10) THEN
412  ctarg = 2.
413  IF (ifsplit) ctarg = 2.
414  ENDIF
415 C
416  kmax = 5
417  nbdinp = 3
418  nsskip = 2
419  ifskip = .true.
420  ifmodp = .false.
421 C
422  ELSEIF (l.EQ.2) THEN
423 C
424  prelax = 1.e-3
425  IF (ifsplit) prelax = 1.e-5
426 C
427  ctarg = 1.
428  IF (ifsplit) ctarg = 1.
429  IF (lx1.GE.10) THEN
430  ctarg = 2.
431  IF (ifsplit) ctarg = 2.
432  ENDIF
433 C
434  kmax = 5
435  nbdinp = 3
436  nsskip = 2
437  ifskip = .true.
438  ifmodp = .false.
439 C
440  ELSE
441  ENDIF
442  CALL setchar
443  RETURN
444  END
445 C
446  SUBROUTINE chkssvt
447 C-----------------------------------------------------------------------
448 C
449 C Check for global steady state (velocity and temp/passive scalar)
450 C
451 C-----------------------------------------------------------------------
452  include 'SIZE'
453  include 'INPUT'
454  include 'TSTEP'
455  include 'STEADY'
456 C
457  IF (ifflow) THEN
458  imesh = 1
459  ifield = 1
460  CALL chkssv
461  ENDIF
462 C
463  imesh = 2
464  DO 100 ifield=2,nfield
465  CALL chksst
466  100 CONTINUE
467 C
468  ifssvt = .true.
469  ifexvt = .true.
470  mfield = 1
471  IF (.NOT.ifflow) mfield=2
472  DO 200 ifield=mfield,nfield
473  IF(.NOT.ifstst(ifield)) ifssvt = .false.
474  IF(.NOT.ifextr(ifield)) ifexvt = .false.
475  200 CONTINUE
476  RETURN
477  END
478 C
479  SUBROUTINE chkssv
480 C--------------------------------------------------------------------
481 C
482 C Check steady state for velocity
483 C
484 C--------------------------------------------------------------------
485  include 'SIZE'
486  include 'SOLN'
487  include 'MASS'
488  include 'INPUT'
489  include 'EIGEN'
490  include 'TSTEP'
491  include 'STEADY'
492  COMMON /ctolpr/ divex
493  COMMON /cprint/ ifprint
494  LOGICAL IFPRINT
495 C
496  COMMON /scrss2/ dv1(lx1,ly1,lz1,lelv)
497  $ , dv2(lx1,ly1,lz1,lelv)
498  $ , dv3(lx1,ly1,lz1,lelv)
499  COMMON /scruz/ w1(lx1,ly1,lz1,lelv)
500  $ , w2(lx1,ly1,lz1,lelv)
501  $ , w3(lx1,ly1,lz1,lelv)
502  $ , bdivv(lx2,ly2,lz2,lelv)
503  COMMON /scrmg/ t1(lx1,ly1,lz1,lelv)
504  $ , t2(lx1,ly1,lz1,lelv)
505  $ , t3(lx1,ly1,lz1,lelv)
506  $ , divv(lx2,ly2,lz2,lelv)
507  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
508  $ , h2(lx1,ly1,lz1,lelv)
509 C
510  CALL opsub3 (dv1,dv2,dv3,vx,vy,vz,vxlag,vylag,vzlag)
511  CALL normvc (dvnnh1,dvnnsm,dvnnl2,dvnnl8,dv1,dv2,dv3)
512  intype = -1
513  CALL sethlm (h1,h2,intype)
514  CALL ophx (w1,w2,w3,dv1,dv2,dv3,h1,h2)
515  CALL opdssum(w1,w2,w3)
516  CALL opmask (w1,w2,w3)
517  CALL opcolv3(t1,t2,t3,w1,w2,w3,binvm1)
518  CALL ophx (w1,w2,w3,t1,t2,t3,h1,h2)
519  CALL opcol2 (w1,w2,w3,dv1,dv2,dv3)
520  ntot1 = lx1*ly1*lz1*nelv
521  usnrm1 = glsum(w1,ntot1)
522  usnrm2 = glsum(w2,ntot1)
523  usnrm3 = 0.
524  IF (ldim.EQ.3) usnrm3 = glsum(w3,ntot1)
525  usnorm = sqrt( (usnrm1+usnrm2+usnrm3)/volvm1 )
526 C
527  ntot2 = lx2*ly2*lz2*nelv
528  CALL opdiv (bdivv,vx,vy,vz)
529  CALL col3 (divv,bdivv,bm2inv,ntot2)
530  dnorm = sqrt(glsc2(divv,bdivv,ntot2)/volvm2)
531 C
532  tolold = tolps
533  CALL settolv
534  tolhv3 = tolhv*(ldim)
535  IF (ifstrs) tolhv3 = tolhv
536  IF (nio.EQ.0 .AND. ifprint) THEN
537  WRITE (6,*) 'USNORM, TOLHV',usnorm,tolhv3
538  WRITE (6,*) 'DNORM, TOLPS',dnorm,tolps
539  ENDIF
540  IF (dnorm.GT.(1.1*divex).AND.divex.GT.0.
541  $ .AND.tolpdf.EQ.0.) tolpdf = 5.*dnorm
542  usrel = usnorm/tolhv3
543  drel = dnorm/tolps
544 C
545  IF (tolrel.GT.0.) THEN
546  exfac = .3/tolrel
547  ELSE
548  WRITE (6,*) 'WARNING: TOLREL=0. Please modify *.rea'
549  call exitt
550  ENDIF
551  ifextr(ifield) = .false.
552 c IF ((USREL.LT.EXFAC).OR.(TIME.GT.TXNEXT(IFIELD)))
553 c $ IFEXTR(IFIELD) = .TRUE.
554  IF (usrel.LT.exfac) ifextr(ifield) = .true.
555  if (nio.eq.0 .and. ifprint)
556  $WRITE (6,*) 'Tau, Txnext ',ifield,tauss(ifield),txnext(ifield)
557 C
558  ifstst(ifield) = .false.
559  uslim = 2.*tolhv3
560  dlim = 2.*tolps
561  IF (usnorm.LT.uslim.AND.dnorm.LT.dlim .AND. .NOT.ifsplit)
562  $ ifstst(ifield) = .true.
563  IF (usnorm.LT.uslim .AND. ifsplit)
564  $ ifstst(ifield) = .true.
565 C
566  RETURN
567  END
568 C
569  SUBROUTINE chksst
570 C----------------------------------------------------------------------
571 C
572 C Check for steady state for temperature/passive scalar
573 C
574 C----------------------------------------------------------------------
575  include 'SIZE'
576  include 'SOLN'
577  include 'MASS'
578  include 'TSTEP'
579  include 'STEADY'
580  COMMON /scruz/ deltat(lx1,ly1,lz1,lelt)
581  $ , wa(lx1,ly1,lz1,lelt)
582  $ , wb(lx1,ly1,lz1,lelt)
583  COMMON /scrvh/ h1(lx1,ly1,lz1,lelt)
584  $ , h2(lx1,ly1,lz1,lelt)
585  COMMON /cprint/ ifprint
586  LOGICAL IFPRINT
587 C
588  ntot = lx1*ly1*lz1*nelt
589  CALL sub3 (deltat(1,1,1,1),t(1,1,1,1,ifield-1),
590  $ tlag(1,1,1,1,1,ifield-1),ntot)
591  CALL normsc (dvnnh1,dvnnsm,dvnnl2,dvnnl8,deltat,imesh)
592  intype = -1
593  isd = 1
594  CALL sethlm (h1,h2,intype)
595  CALL axhelm (wa,deltat,h1,h2,imesh,isd)
596  CALL dssum (wa,lx1,ly1,lz1)
597  CALL col2 (wa,tmask(1,1,1,1,ifield-1),ntot)
598  CALL col3 (wb,wa,bintm1,ntot)
599  CALL axhelm (wa,wb,h1,h2,imesh,isd)
600  CALL col2 (wa,deltat,ntot)
601  usnorm = sqrt(glsum(wa,ntot)/voltm1)
602 C
603  CALL settolt
604  IF (nio.EQ.0 .AND. ifprint)
605  $WRITE (6,*) 'USNORM, TOLHT',usnorm,tolht(ifield)
606  usrel = usnorm/tolht(ifield)
607 C
608  IF (tolrel.GT.0.) THEN
609  exfac = .3/tolrel
610  ELSE
611  WRITE (6,*) 'WARNING: TOLREL=0. Please modify *.rea'
612  call exitt
613  ENDIF
614  ifextr(ifield) = .false.
615 c IF ((USREL.LT.EXFAC).OR.(TIME.GT.TXNEXT(IFIELD)))
616 c $ IFEXTR(IFIELD) = .TRUE.
617  IF (usrel.LT.exfac) ifextr(ifield) = .true.
618  IF(nio.EQ.0 .AND. ifprint)
619  $WRITE (6,*) 'Tau, Txnext ',ifield,tauss(ifield),txnext(ifield)
620 C
621  ifstst(ifield) = .false.
622  uslim = 2.*tolht(ifield)
623  IF (usnorm.LT.uslim) ifstst(ifield) = .true.
624 C
625  RETURN
626  END
627 C
628  SUBROUTINE ssnormd (DV1,DV2,DV3)
629  include 'SIZE'
630  include 'MASS'
631  include 'TSTEP'
632  include 'STEADY'
633  REAL DV1(1),DV2(1),DV3(1)
634  CALL normvc (dvdfh1,dvdfsm,dvdfl2,dvdfl8,dv1,dv2,dv3)
635  RETURN
636  END
637 C
638  SUBROUTINE ssnormp (DV1,DV2,DV3)
639  include 'SIZE'
640  include 'TSTEP'
641  include 'STEADY'
642  REAL DV1(1),DV2(1),DV3(1)
643  CALL normvc (dvprh1,dvprsm,dvprl2,dvprl8,dv1,dv2,dv3)
644  RETURN
645  END
646 C
647  SUBROUTINE settolv
648 C-------------------------------------------------------------------
649 C
650 C Set tolerances for velocity solver
651 C
652 C-------------------------------------------------------------------
653  include 'SIZE'
654  include 'INPUT'
655  include 'EIGEN'
656  include 'MASS'
657  include 'TSTEP'
658  include 'SOLN'
659  REAL LENGTH
660 C
661  ntot = lx1*ly1*lz1*nelfld(ifield)
662  avvisc = glmin(vdiff(1,1,1,1,ifield),ntot)
663  avdens = glmax(vtrans(1,1,1,1,ifield),ntot)
664 C
665  IF (iftran) THEN
666  IF (istep.EQ.1) vnorm = vnrml8
667  IF (istep.GT.1) vnorm = vnrmsm
668  IF (vnorm.EQ.0.) vnorm = tolabs
669  factor = 1.+(avdens/(eigaa*avvisc*dt))
670  ELSE
671  vnorm = vnrml8
672  IF (vnorm.EQ.0.) vnorm = tolabs
673  factor = 1.
674  ENDIF
675 C
676  tolps = tolrel*vnorm * sqrt(eigas)/(4.*factor)
677  tolhv = tolrel*vnorm * sqrt(eigaa)*avvisc/2.
678  tolhv = tolhv/3.
679  IF (.NOT.iftran .AND. .NOT.ifnav) tolhv = tolhv/10.
680  tolhr = tolhv
681  tolhs = tolhv
682 C
683 C Non-zero default pressure tolerance
684 C NOTE: This tolerance may change due to precision problems.
685 C See subroutine CHKSSV
686 C
687  IF (tolpdf.NE.0.) tolps = tolpdf
688 C
689  RETURN
690  END
691 C
692  SUBROUTINE settolt
693 C-------------------------------------------------------------------
694 C
695 C Set tolerances for temerature/passive scalar solver
696 C
697 C-------------------------------------------------------------------
698  include 'SIZE'
699  include 'INPUT'
700  include 'EIGEN'
701  include 'MASS'
702  include 'TSTEP'
703  include 'SOLN'
704  REAL LENGTH
705 C
706  ntot = lx1*ly1*lz1*nelfld(ifield)
707  avcond = glmin(vdiff(1,1,1,1,ifield),ntot)
708 C
709  IF (iftran) THEN
710  IF (istep.EQ.1) tnorm = tnrml8(ifield-1)
711  IF (istep.GT.1) tnorm = tnrmsm(ifield-1)
712  IF (tnorm.EQ.0.) tnorm = tolabs
713  ELSE
714  tnorm = tnrml8(ifield-1)
715  IF (tnorm.EQ.0.) tnorm = tolabs
716  ENDIF
717 C
718  tolht(ifield) = tolrel*tnorm * sqrt(eigaa)*avcond
719 C
720  RETURN
721  END
722 C
723  SUBROUTINE chktolp (TOLMIN)
724  include 'SIZE'
725  include 'SOLN'
726  include 'MASS'
727  include 'TSTEP'
728  COMMON /scrmg/ divfld(lx2,ly2,lz2,lelv)
729  $ , work(lx2,ly2,lz2,lelv)
730  ntot2 = lx2*ly2*lz2*nelv
731  CALL opdiv (divfld,vx,vy,vz)
732  CALL col3 (work,divfld,bm2inv,ntot2)
733  CALL col2 (work,divfld,ntot2)
734  divv = sqrt(glsum(work,ntot2)/volvm2)
735 C
736  ifield = 1
737  CALL settolv
738  tolmin = divv/100.
739  IF (tolmin.LT.tolps) tolmin = tolps
740  RETURN
741  END
742 C
743  SUBROUTINE setchar
744 C-----------------------------------------------------------------------
745 C
746 C If characteristics, need number of sub-timesteps (DT/DS).
747 C Current sub-timeintegration scheme: RK4.
748 C If not characteristics, i.e. standard semi-implicit scheme,
749 C check user-defined Courant number.
750 C
751 C----------------------------------------------------------------------
752  include 'SIZE'
753  include 'INPUT'
754  include 'TSTEP'
755 
756  REAL MAX_CFL_RK4
757  DATA max_cfl_rk4 /2.0/ ! stability limit for RK4 including safety factor
758 C
759  IF (ifchar) THEN
760  ict = max(int(ctarg/max_cfl_rk4),1)
761  ntaubd = ict
762  dct = ctarg - ict*max_cfl_rk4
763  IF (dct.GT.0.1) ntaubd = ntaubd+1
764  IF (nio.EQ.0) write(6,*) 'RK4 substeps:', ntaubd
765  ELSE
766  ntaubd = 0
767  IF (ctarg.GT.0.5) THEN
768  IF (nio.EQ.0)
769  $ WRITE (6,*) 'Reset the target Courant number to .5'
770  ctarg = 0.5
771  ENDIF
772  ENDIF
773 C
774  RETURN
775  END
776 C
777  SUBROUTINE project
778 C--------------------------------------------------------------------
779 C
780 C Project current solution onto the closest incompressible field
781 C
782 C--------------------------------------------------------------------
783  include 'SIZE'
784  include 'TOTAL'
785  COMMON /scrns/ w1(lx1,ly1,lz1,lelv)
786  $ , w2(lx1,ly1,lz1,lelv)
787  $ , w3(lx1,ly1,lz1,lelv)
788  $ , dv1(lx1,ly1,lz1,lelv)
789  $ , dv2(lx1,ly1,lz1,lelv)
790  $ , dv3(lx1,ly1,lz1,lelv)
791  $ , respr(lx2,ly2,lz2,lelv)
792  COMMON /scrvh/ h1(lx1,ly1,lz1,lelv)
793  $ , h2(lx1,ly1,lz1,lelv)
794 C
795  IF (nio.EQ.0) WRITE(6,5)
796  5 FORMAT(/,' Project',/)
797 C
798  ntot1 = lx1*ly1*lz1*nelv
799  ntot2 = lx2*ly2*lz2*nelv
800  intype = 1
801  CALL rzero (h1,ntot1)
802  CALL rone (h2,ntot1)
803  CALL opdiv (respr,vx,vy,vz)
804  CALL chsign (respr,ntot2)
805  CALL ortho (respr)
806  CALL uzawa (respr,h1,h2,intype,icg)
807  CALL opgradt (w1,w2,w3,respr)
808  CALL opbinv (dv1,dv2,dv3,w1,w2,w3,h2)
809  CALL opadd2 (vx,vy,vz,dv1,dv2,dv3)
810  RETURN
811  END
exitt
void exitt()
Definition: comm_mpi.f:604
lagscal
subroutine lagscal
Definition: conduct.f:295
settime
subroutine settime
Definition: drive2.f:540
geneig
subroutine geneig(igeom)
Definition: drive2.f:606
meshv
subroutine meshv(igeom)
Definition: drive2.f:812
fluid
subroutine fluid(igeom)
Definition: drive2.f:658
comment
subroutine comment
Definition: drive2.f:86
heat
subroutine heat(igeom)
Definition: drive2.f:744
gengeom
subroutine gengeom(igeom)
Definition: drive2.f:371
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
glmin
function glmin(a, n)
Definition: math.f:973
glsc2
function glsc2(x, y, n)
Definition: math.f:794
col2
subroutine col2(a, b, n)
Definition: math.f:564
glsum
function glsum(x, n)
Definition: math.f:861
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
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
opgradt
subroutine opgradt(outx, outy, outz, inpfld)
Definition: navier1.f:4096
uzawa
subroutine uzawa(rcg, h1, h2, h2inv, intype, iter)
Definition: navier1.f:2827
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
opcol2
subroutine opcol2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2454
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
opcolv3
subroutine opcolv3(a1, a2, a3, b1, b2, b3, c)
Definition: navier1.f:2380
opmask
subroutine opmask(res1, res2, res3)
Definition: navier1.f:2314
lagvel
subroutine lagvel
Definition: navier1.f:1930
normvc
subroutine normvc(h1, semi, l2, linf, x1, x2, x3)
Definition: navier1.f:2083
ophx
subroutine ophx(out1, out2, out3, inp1, inp2, inp3, h1, h2)
Definition: navier1.f:743
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
opsub3
subroutine opsub3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
Definition: navier1.f:2370
normsc
subroutine normsc(h1, semi, l2, linf, x, imesh)
Definition: navier1.f:2026
tnorm
function tnorm(temp)
Definition: pertsupport.f:58
prepost
subroutine prepost(ifdoin, prefin)
Definition: prepost.f:91
ssnormd
subroutine ssnormd(DV1, DV2, DV3)
Definition: ssolv.f:629
modprop
subroutine modprop
Definition: ssolv.f:279
ssnormp
subroutine ssnormp(DV1, DV2, DV3)
Definition: ssolv.f:639
project
subroutine project
Definition: ssolv.f:778
sstest
subroutine sstest(ISSS)
Definition: ssolv.f:2
settolv
subroutine settolv
Definition: ssolv.f:648
gostep
subroutine gostep
Definition: ssolv.f:240
chktolp
subroutine chktolp(TOLMIN)
Definition: ssolv.f:724
ssinit
subroutine ssinit(KMAX)
Definition: ssolv.f:27
filllag
subroutine filllag
Definition: ssolv.f:158
mkarr
subroutine mkarr(X)
Definition: ssolv.f:348
settolt
subroutine settolt
Definition: ssolv.f:693
go1step
subroutine go1step(X, Y, NVEC)
Definition: ssolv.f:198
gonstep
subroutine gonstep(N, ITEST)
Definition: ssolv.f:177
ssparam
subroutine ssparam(KMAX, L)
Definition: ssolv.f:391
chkext
subroutine chkext(IFACCX, Z, S)
Definition: ssolv.f:81
chkssvt
subroutine chkssvt
Definition: ssolv.f:447
chkssv
subroutine chkssv
Definition: ssolv.f:480
setchar
subroutine setchar
Definition: ssolv.f:744
mkvec
subroutine mkvec(X)
Definition: ssolv.f:305
chksst
subroutine chksst
Definition: ssolv.f:570
unorm
subroutine unorm
Definition: subs1.f:352
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
setsolv
subroutine setsolv
Definition: subs1.f:1083
setprop
subroutine setprop
Definition: subs1.f:2618