|
tick_geo.f.html |
 |
|
Source file: tick_geo.f
|
|
Directory: /home/rjl/git/claworg/clawpack-4.x/geoclaw/2d/lib
|
|
Converted: Sat Aug 6 2011 at 21:59:28
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c -------------------------------------------------------------
c
subroutine tick(nvar,iout,nstart,nstop,cut,vtime,time,ichkpt,
& naux,nout,tout,tchk,t0,rest)
c
use geoclaw_module
implicit double precision (a-h,o-z)
include "call.i"
logical vtime, dumpout, dumpchk, rest
dimension dtnew(maxlv), ntogo(maxlv), tlevel(maxlv)
dimension tout(nout)
dimension tchk(maxout)
c
c :::::::::::::::::::::::::::: TICK :::::::::::::::::::::::::::::
c main driver routine. controls:
c integration of all grids.
c error estimation / regridding
c output counting
c updating of fine to coarse grids
c parameters:
c nstop = # of coarse grid time steps to be taken
c iout = output interval every 'iout' coarse time steps
c (if 0, not used - set to inf.)
c vtime = true for variable timestep, calculated each coarse step
c
c integration strategy is to advance a fine grid until it catches
c up to the coarse grid. this strategy is applied recursively.
c coarse grid goes first.
c
c nsteps: used to count how number steps left for a level to be
c integrated before it catches up with the next coarser level.
c ncycle: counts number of coarse grid steps = # cycles.
c
c icheck: counts the number of steps (incrementing by 1
c each step) to keep track of when that level should
c have its error estimated and finer levels should be regridded.
c ::::::::::::::::::::::::::::::::::::;::::::::::::::::::::::::::
c
integer verbosity_regrid
c # Eventually add this to new input params?
c # For now use verbosity, or set to 0 to suppress printing regrid info:
verbosity_regrid = method(4)
ncycle = nstart
call setbestsrc() ! need at very start of run, including restart
if ( iout .eq. 0) iout = iinfinity
if (ichkpt .eq. 0) ichkpt = iinfinity
nextout = 1
if (nout .gt. 0) then
tfinal = tout(nout)
c if this is a restart, make sure output times start after restart time
if (nstart .gt. 0) then
do ii = 1, nout
if (tout(ii) .gt. time) then
nextout = ii
go to 2
endif
end do
2 continue
endif
else
tfinal = rinfinity
endif
nextchk = 1
if (nstart .gt. 0 .and. ichkpt.lt.0) then
c if this is a restart, make sure chkpt times start after restart time
do ii = 1, -ichkpt
if (tchk(ii) .gt. time) then
nextchk = ii
go to 3
endif
enddo
3 continue
endif
tlevel(1) = time
do 5 i = 2, mxnest
tlevel(i) = tlevel(1)
5 continue
c
c ------ start of coarse grid integration loop. ------------------
c
20 if (ncycle .ge. nstop .or. time .ge. tfinal) goto 999
if (nextout .le. nout) then
outtime = tout(nextout)
else
outtime = rinfinity
endif
if (nextchk .le. -ichkpt) then
chktime = tchk(nextchk)
else
chktime = rinfinity
endif
c if (time.lt.outtime .and. time + possk(1) .ge. outtime) then
if (time.lt.outtime .and. time+1.001*possk(1) .ge. outtime) then
c apr 2010 mjb: modified so allow slightly larger timestep to
c hit output time exactly, instead of taking minuscule timestep
c should still be stable since increase dt in only 3rd digit.
c ## adjust time step to hit outtime exactly, and make output
oldposs = possk(1)
possk(1) = outtime - time
c write(*,*)" old possk is ", possk(1)
diffdt = oldposs - possk(1) ! if positive new step is smaller
if (method(4).ge.2) then ! verbosity test to
write(*,122) diffdt,outtime ! notify of change
122 format(" Adjusting timestep by ",e10.3,
. " to hit output time of ",e12.6)
c write(*,*)" new possk is ", possk(1)
if (diffdt .lt. 0.) then ! new step is slightly larger
pctIncrease = -100.*diffdt/oldposs ! minus sign to make whole expr. positive
write(*,123) pctIncrease
123 format(" New step is ",e8.2," % larger.",
. " Should still be stable")
endif
endif
do 12 i = 2, mxnest
12 possk(i) = possk(i-1) / kratio(i-1)
nextout = nextout + 1
dumpout = .true.
else
dumpout = .false.
endif
if (time.lt.chktime .and. time + possk(1) .ge. chktime) then
c ## adjust time step to hit chktime exactly, and do checkpointing
possk(1) = chktime - time
do 13 i = 2, mxnest
13 possk(i) = possk(i-1) / kratio(i-1)
nextchk = nextchk + 1
dumpchk = .true.
else
dumpchk = .false.
endif
c take output stuff from here - put it back at end.
c
level = 1
ntogo(level) = 1
do 10 i = 1, maxlv
10 dtnew(i) = rinfinity
c
c ------------- regridding time? ---------
c
c check if either
c (i) this level should have its error estimated before being advanced
c (ii) this level needs to provide boundary values for either of
c next 2 finer levels to have their error estimated.
c this only affects two grid levels higher, occurs because
c previous time step needs boundary vals for giant step.
c no error estimation on finest possible grid level
c
60 continue
if (icheck(level) .ge. kcheck) then
lbase = level
else if (level+1 .ge. mxnest) then
go to 90
else if (icheck(level+1) .ge. kcheck) then
lbase = level+1
else if (level+2 .ge. mxnest) then
go to 90
else if (icheck(level+2) .ge. kcheck) then
lbase = level+2
else
go to 90
endif
if (lbase .eq. mxnest .or. lbase .gt. lfine) go to 70
c
c regrid level 'lbase+1' up to finest level.
c level 'lbase' stays fixed.
c
if (rprint) write(outunit,101) lbase
101 format(8h level ,i5,32h stays fixed during regridding )
call regrid(nvar,lbase,cut,naux,t0)
call setbestsrc() ! need at every grid change
c call conck(1,nvar,time,rest)
c call outtre(lstart(lbase+1),.true.,nvar,naux)
c note negative time to signal regridding output in plots
c call valout(lbase,lfine,-tlevel(lbase),nvar,naux)
c
c maybe finest level in existence has changed. reset counters.
c
if (rprint .and. lbase .lt. lfine) then
call outtre(lstart(lbase+1),.false.,nvar,naux)
endif
70 continue
do 80 i = lbase, lfine
80 icheck(i) = 0
do 81 i = lbase+1, lfine
81 tlevel(i) = tlevel(lbase)
c
c MJB: modified to check level where new grids start, which is lbase+1
if (verbosity_regrid.ge.lbase+1) then
do levnew = lbase+1,lfine
write(6,1006) intratx(levnew-1),intraty(levnew-1),
& kratio(levnew-1),levnew
1006 format(' Refinement ratios... in x:', i3,
& ' in y:',i3,' in t:',i3,' for level ',i4)
end do
endif
c ------- done regridding --------------------
c
c integrate all grids at level 'level'.
c
90 continue
call advanc(level,nvar,dtlevnew,vtime,naux)
c # rjl modified 6/17/05 to print out *after* advanc and print cfl
c # rjl & mjb changed to cfl_level, 3/17/10
timenew = tlevel(level)+possk(level)
if (tprint) then
write(outunit,100)level,cfl_level,possk(level),timenew
endif
if (method(4).ge.level) then
write(6,100)level,cfl_level,possk(level),timenew
endif
100 format(' AMRCLAW: level ',i2,' CFL = ',e8.3,
& ' dt = ',e10.4, ' final t = ',e12.6)
c # to debug individual grid updates...
c call valout(level,level,time,nvar,naux)
c
c done with a level of integration. update counts, decide who next.
c
ntogo(level) = ntogo(level) - 1
dtnew(level) = dmin1(dtnew(level),dtlevnew)
tlevel(level) = tlevel(level) + possk(level)
icheck(level) = icheck(level) + 1
c
if (level .lt. lfine) then
level = level + 1
c # check if should adjust finer grid time step to start wtih
if (((possk(level-1) - dtnew(level-1))/dtnew(level-1)) .gt.
. .05) then
dttemp = dtnew(level-1)/kratio(level-1)
ntogo(level) = (tlevel(level-1)-tlevel(level))/dttemp+.9
else
ntogo(level) = kratio(level-1)
endif
possk(level) = possk(level-1)/ntogo(level)
go to 60
endif
c
105 if (level .eq. 1) go to 110
if (ntogo(level) .gt. 0) then
c same level goes again. check for ok time step
106 if ((possk(level)-dtnew(level))/dtnew(level)
. .gt. .05) then
write(6,*)" ***adjusting timestep for level ",level
write(6,*)" old ntogo dt ",ntogo(level),possk(level)
c adjust time steps for this and finer levels
ntogo(level) = ntogo(level) + 1
possk(level) = (tlevel(level-1)-tlevel(level))/
. ntogo(level)
if (varRefTime) then
kratio(level-1) = ceiling(possk(level-1) /
. possk(level))
endif
write(6,*)" new ntogo dt ",ntogo(level),possk(level)
go to 106
endif
go to 60
else
level = level - 1
call update(level,nvar)
endif
go to 105
c
c --------------one complete coarse grid integration cycle done. -----
c
c time for output? done with the whole thing?
c
110 continue
time = time + possk(1)
ncycle = ncycle + 1
call conck(1,nvar,time,rest)
if ((ichkpt.gt.0 .and. mod(ncycle,ichkpt).eq.0)
& .or. dumpchk) then
call check(ncycle,time,nvar,naux)
endif
if ((mod(ncycle,iout).eq.0) .or. dumpout) then
call valout(1,lfine,time,nvar,naux)
if (printout) call outtre(mstart,.true.,nvar,naux)
endif
if ( .not. vtime) go to 201
c
c ## adjust time steps if variable time step and/or variable refinement ratios in time
c
if (.not. varRefTime) then
c
c find new dt for next cycle (passed back from integration routine).
do 115 i = 2, lfine
ii = lfine+1-i
dtnew(ii) = min(dtnew(ii),dtnew(ii+1)*kratio(ii))
115 continue
possk(1) = dtnew(1)
do 120 i = 2, mxnest
120 possk(i) = possk(i-1) / kratio(i-1)
else ! since refinement ratio in time can change need to set new timesteps in different order
c ! use same alg. as when setting refinement when first make new fine grids
possk(1) = dtnew(1)
do 125 i = 2, lfine
if (dtnew(i) .gt. possk(i-1)) then
kratio(i-1) = 1 ! cant have larger timestep than parent level
possk(i) = possk(i-1)
else
kratio(i-1) = ceiling(possk(i-1)/dtnew(i)) ! round up for stable integer ratio
possk(i) = possk(i-1)/kratio(i-1) ! set exact timestep on this level
endif
125 continue
endif
201 go to 20
c
999 continue
c
c # computation is complete to final time or requested number of steps
c
if (ncycle .ge. nstop .and. tfinal .lt. rinfinity) then
c # warn the user that calculation finished prematurely
write(outunit,102) nstop
write(6,102) nstop
102 format('*** Computation halted after nv(1) = ',i8,
& ' steps on coarse grid')
endif
c
c # final output (unless we just did it above)
c
if (.not. dumpout) then
if (mod(ncycle,iout).ne.0) then
call valout(1,lfine,time,nvar,naux)
if (printout) call outtre(mstart,.true.,nvar,naux)
endif
endif
c # checkpoint everything for possible future restart
c # (unless we just did it based on ichkpt or dumpchk)
c
c # don't checkpoint at all if user set ichkpt=0
if (ichkpt .gt. 0) then
if ((ncycle/ichkpt)*ichkpt .ne. ncycle) then
call check(ncycle,time,nvar,naux)
endif
endif
if (ichkpt .lt. 0) then
if (.not. dumpchk) then
call check(ncycle,time,nvar,naux)
endif
endif
return
end