source: CIVL/examples/fortran/nek5000/core/ssolv.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: 21.2 KB
Line 
1 SUBROUTINE SSTEST (ISSS)
2C------------------------------------------------------------------------
3C
4C Test if Steady State Solver should be activated.
5C
6C------------------------------------------------------------------------
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
25C
26 SUBROUTINE SSINIT (KMAX)
27C-----------------------------------------------------------
28C
29C Initialize steady state solver
30C
31C-----------------------------------------------------------
32 INCLUDE 'SIZE'
33 INCLUDE 'INPUT'
34 INCLUDE 'EIGEN'
35 INCLUDE 'TSTEP'
36 INCLUDE 'STEADY'
37C
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
54C
55 NBDINP = 1.
56 TIME = 0.
57 DT = 0.
58 DTINIT = TAUMIN/5.
59 NSTEPS = 10000
60 IOSTEP = 10000
61C
62 IFMODP = .TRUE.
63 IFSKIP = .TRUE.
64 NSSKIP = 1
65C
66 PRELAX = 1.E-1
67 IF (IFSPLIT) PRELAX = 1.E-4
68C
69 IFSSVT = .FALSE.
70 IFEXVT = .FALSE.
71C
72 KMAX = 5
73C
74 CALL SETCHAR
75C
76C
77 RETURN
78 END
79C
80 SUBROUTINE CHKEXT (IFACCX,Z,S)
81C------------------------------------------------------------------
82C
83C Accept extrapolation?
84C
85C------------------------------------------------------------------
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)
93C
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
104C
105 CALL MKVEC (Z)
106 CALL MKARR (S)
107C
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
118C
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
136C
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
153C
154 RETURN
155 END
156C
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
175C
176 SUBROUTINE GONSTEP (N,ITEST)
177C----------------------------------------------------------------
178C
179C Do N steps; return if steady state
180C
181C----------------------------------------------------------------
182 INCLUDE 'SIZE'
183 INCLUDE 'INPUT'
184 INCLUDE 'STEADY'
185 EXTERNAL GOSTEP
186C
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
193C
194 RETURN
195 END
196C
197 SUBROUTINE GO1STEP (X,Y,NVEC)
198C----------------------------------------------------------------
199C
200C Advance one (or more) time step(s)
201C
202C----------------------------------------------------------------
203 INCLUDE 'SIZE'
204 INCLUDE 'TOTAL'
205 real*8 X(1), Y(1)
206C
207 CALL MKARR (X)
208 IF (.NOT.IFSKIP) NJSTEP=1
209 IF ( IFSKIP) NJSTEP=NSSKIP
210C
211 DO 9000 JSTEP=1,NJSTEP
212C
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
230C
231 9000 CONTINUE
232C
233 IF (ISTEP.GT.1) CALL CHKSSVT
234 CALL MKVEC (Y)
235C
236 RETURN
237 END
238C
239 SUBROUTINE GOSTEP
240C----------------------------------------------------------------
241C
242C Advance one (or more) time step(s)
243C
244C----------------------------------------------------------------
245 INCLUDE 'SIZE'
246 INCLUDE 'TOTAL'
247C
248 IF (.NOT.IFSKIP) NJSTEP=1
249 IF ( IFSKIP) NJSTEP=NSSKIP
250C
251 DO 9000 JSTEP=1,NJSTEP
252C
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
270C
271 9000 CONTINUE
272C
273 IF (ISTEP.GT.1) CALL CHKSSVT
274C
275 RETURN
276 END
277C
278 SUBROUTINE MODPROP
279C------------------------------------------------------------------
280C
281C Modify the properties
282C
283C------------------------------------------------------------------
284 INCLUDE 'SIZE'
285 INCLUDE 'SOLN'
286 INCLUDE 'TSTEP'
287 INCLUDE 'STEADY'
288 INCLUDE 'INPUT'
289C
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)
297c if (nid.eq.0)
298c $ write (6,*) '.......... diff = ',IFIELD,vdiff(1,1,1,1,IFIELD)
299 100 CONTINUE
300C
301 RETURN
302 END
303C
304 SUBROUTINE MKVEC (X)
305C-------------------------------------------------------------
306C
307C Fill up the vector X with VX, VY, ....
308C
309C-------------------------------------------------------------
310 INCLUDE 'SIZE'
311 INCLUDE 'SOLN'
312 INCLUDE 'INPUT'
313 INCLUDE 'TSTEP'
314 real*8 X(1)
315C
316 NTOTV = lx1*ly1*lz1*NELV
317C
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
332C
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
343C
344 RETURN
345 END
346C
347 SUBROUTINE MKARR (X)
348C------------------------------------------------------------------
349C
350C Split the vector X into VX, VY, .....
351C
352C------------------------------------------------------------------
353 INCLUDE 'SIZE'
354 INCLUDE 'SOLN'
355 INCLUDE 'TSTEP'
356 INCLUDE 'INPUT'
357 real*8 X(1)
358C
359 NTOTV = lx1*ly1*lz1*NELV
360C
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
375C
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
386C
387 RETURN
388 END
389C
390 SUBROUTINE SSPARAM (KMAX,L)
391C------------------------------------------------------------------------------
392C
393C Set steady state parameters
394C
395C------------------------------------------------------------------------------
396 INCLUDE 'SIZE'
397 INCLUDE 'TSTEP'
398 INCLUDE 'INPUT'
399 INCLUDE 'STEADY'
400C
401 IF (L.EQ.0) THEN
402 CALL SSINIT (KMAX)
403 ELSEIF (L.EQ.1) THEN
404 ISTEP = 0
405C
406 PRELAX = 1.E-2
407 IF (IFSPLIT) PRELAX = 1.E-5
408C
409 CTARG = 1.
410 IF (IFSPLIT) CTARG = 1.
411 IF (lx1.GE.10) THEN
412 CTARG = 2.
413 IF (IFSPLIT) CTARG = 2.
414 ENDIF
415C
416 KMAX = 5
417 NBDINP = 3
418 NSSKIP = 2
419 IFSKIP = .TRUE.
420 IFMODP = .FALSE.
421C
422 ELSEIF (L.EQ.2) THEN
423C
424 PRELAX = 1.E-3
425 IF (IFSPLIT) PRELAX = 1.E-5
426C
427 CTARG = 1.
428 IF (IFSPLIT) CTARG = 1.
429 IF (lx1.GE.10) THEN
430 CTARG = 2.
431 IF (IFSPLIT) CTARG = 2.
432 ENDIF
433C
434 KMAX = 5
435 NBDINP = 3
436 NSSKIP = 2
437 IFSKIP = .TRUE.
438 IFMODP = .FALSE.
439C
440 ELSE
441 ENDIF
442 CALL SETCHAR
443 RETURN
444 END
445C
446 SUBROUTINE CHKSSVT
447C-----------------------------------------------------------------------
448C
449C Check for global steady state (velocity and temp/passive scalar)
450C
451C-----------------------------------------------------------------------
452 INCLUDE 'SIZE'
453 INCLUDE 'INPUT'
454 INCLUDE 'TSTEP'
455 INCLUDE 'STEADY'
456C
457 IF (IFFLOW) THEN
458 IMESH = 1
459 IFIELD = 1
460 CALL CHKSSV
461 ENDIF
462C
463 IMESH = 2
464 DO 100 IFIELD=2,NFIELD
465 CALL CHKSST
466 100 CONTINUE
467C
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
478C
479 SUBROUTINE CHKSSV
480C--------------------------------------------------------------------
481C
482C Check steady state for velocity
483C
484C--------------------------------------------------------------------
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
495C
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)
509C
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 )
526C
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)
531C
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
544C
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.
552c IF ((USREL.LT.EXFAC).OR.(TIME.GT.TXNEXT(IFIELD)))
553c $ 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)
557C
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.
565C
566 RETURN
567 END
568C
569 SUBROUTINE CHKSST
570C----------------------------------------------------------------------
571C
572C Check for steady state for temperature/passive scalar
573C
574C----------------------------------------------------------------------
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
587C
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)
602C
603 CALL SETTOLT
604 IF (NIO.EQ.0 .AND. IFPRINT)
605 $WRITE (6,*) 'USNORM, TOLHT',USNORM,TOLHT(IFIELD)
606 USREL = USNORM/TOLHT(IFIELD)
607C
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.
615c IF ((USREL.LT.EXFAC).OR.(TIME.GT.TXNEXT(IFIELD)))
616c $ 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)
620C
621 IFSTST(IFIELD) = .FALSE.
622 USLIM = 2.*TOLHT(IFIELD)
623 IF (USNORM.LT.USLIM) IFSTST(IFIELD) = .TRUE.
624C
625 RETURN
626 END
627C
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
637C
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
646C
647 SUBROUTINE SETTOLV
648C-------------------------------------------------------------------
649C
650C Set tolerances for velocity solver
651C
652C-------------------------------------------------------------------
653 INCLUDE 'SIZE'
654 INCLUDE 'INPUT'
655 INCLUDE 'EIGEN'
656 INCLUDE 'MASS'
657 INCLUDE 'TSTEP'
658 INCLUDE 'SOLN'
659 REAL LENGTH
660C
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)
664C
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
675C
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
682C
683C Non-zero default pressure tolerance
684C NOTE: This tolerance may change due to precision problems.
685C See subroutine CHKSSV
686C
687 IF (TOLPDF.NE.0.) TOLPS = TOLPDF
688C
689 RETURN
690 END
691C
692 SUBROUTINE SETTOLT
693C-------------------------------------------------------------------
694C
695C Set tolerances for temerature/passive scalar solver
696C
697C-------------------------------------------------------------------
698 INCLUDE 'SIZE'
699 INCLUDE 'INPUT'
700 INCLUDE 'EIGEN'
701 INCLUDE 'MASS'
702 INCLUDE 'TSTEP'
703 INCLUDE 'SOLN'
704 REAL LENGTH
705C
706 NTOT = lx1*ly1*lz1*NELFLD(IFIELD)
707 AVCOND = GLMIN (VDIFF(1,1,1,1,IFIELD),NTOT)
708C
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
717C
718 TOLHT(IFIELD) = TOLREL*TNORM * SQRT(EIGAA)*AVCOND
719C
720 RETURN
721 END
722C
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)
735C
736 IFIELD = 1
737 CALL SETTOLV
738 TOLMIN = DIVV/100.
739 IF (TOLMIN.LT.TOLPS) TOLMIN = TOLPS
740 RETURN
741 END
742C
743 SUBROUTINE SETCHAR
744C-----------------------------------------------------------------------
745C
746C If characteristics, need number of sub-timesteps (DT/DS).
747C Current sub-timeintegration scheme: RK4.
748C If not characteristics, i.e. standard semi-implicit scheme,
749C check user-defined Courant number.
750C
751C----------------------------------------------------------------------
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
758C
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
773C
774 RETURN
775 END
776C
777 SUBROUTINE PROJECT
778C--------------------------------------------------------------------
779C
780C Project current solution onto the closest incompressible field
781C
782C--------------------------------------------------------------------
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)
794C
795 IF (NIO.EQ.0) WRITE(6,5)
796 5 FORMAT(/,' Project',/)
797C
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
Note: See TracBrowser for help on using the repository browser.