| 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
|
|---|
