|
setgrd_geo.f.html |
 |
|
Source file: setgrd_geo.f
|
|
Directory: /home/rjl/git/rjleveque/clawpack-4.x/geoclaw/2d/lib
|
|
Converted: Sun May 15 2011 at 19:15:41
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c -----------------------------------------------------------------
c
subroutine setgrd (nvar,cut,naux,dtinit,t0)
c
use geoclaw_module
c
implicit double precision (a-h,o-z)
include "call.i"
dimension spoh(maxlv)
integer verbosity_regrid
logical vtime
data vtime/.false./
c # may as well not bother to calculate time step for error est.
c
c :::::::::::::::::::::::::::: SETGRD :::::::::::::::::::::::::::::::;
c set up the entire tree/grid structure. only at this time t = 0
c can we take advantage of initialization routines.
c remember that regridding/error estimation needs to have two
c time steps of soln. values.
c 6/21/05: added dtinit arg. to allow for better choice of initial timestep
c as discovered by advance/setgrd in first step.
c ::::::::::::::::::::::::::::::::::::::;::::::::::::::::::::::::::::;
c
dtinit = possk(1)
if (mxnest .eq. 1) go to 99
c
levnew = 2
time = t0
verbosity_regrid = method(4)
c
10 if (levnew .gt. mxnest) go to 30
levold = levnew - 1
if (lstart(levold) .eq. 0) go to 30
lbase = levold
lfnew = lbase
c
c set up level to be flagged. need a solution t=0,and t=dt.
c error estimation makes next one at t=2*dt for Richardson.
c
call advanc(levold,nvar,dtlev,vtime,naux)
evol = evol + rvol
rvol = 0.d0
kfac = 1
do i = 1, levold-1
kfac = kfac * kratio(i)
end do
dtinit = min(dtinit, dtlev*kfac)
c dont count it in real integration stats
do 20 level=1,mxnest
20 rvoll(level) = 0.d0
c
c flag, cluster, and make new grids
c
call grdfit(lbase,levold,nvar,naux,cut,time,t0)
if (newstl(levnew) .ne. 0) lfnew = levnew
c
c init new level. after each iteration. fix the data structure
c also reinitalize coarser grids so fine grids can be advanced
c and interpolate correctly for their bndry vals from coarser grids.
c
call ginit(newstl(levnew),.true., nvar, naux, t0)
lstart(levnew) = newstl(levnew)
lfine = lfnew
call ginit(lstart(levold),.false., nvar, naux, t0)
c
c count number of grids on newly created levels (needed for openmp
c parallelization). this is also done in regridding.
c set up numgrids now for new level, since advanc will use it for parallel execution
c
mptr = lstart(levnew)
ngrids = 0
ncells = 0
do while (mptr .gt. 0)
ngrids = ngrids + 1
ncells = ncells + (node(ndihi,mptr)-node(ndilo,mptr)+1)
. * (node(ndjhi,mptr)-node(ndjlo,mptr)+1)
mptr = node(levelptr, mptr)
end do
numgrids(levnew) = ngrids
numcells(levnew) = ncells
if (verbosity_regrid .ge. levnew) then
write(*,100) ngrids,ncells,levnew
100 format("there are ",i4," grids with ",i8,
& " cells at level ", i3)
endif
c
levnew = levnew + 1
go to 10
30 continue
c
c switch location of old and new storage for soln. vals,
c and reset time to 0.0 (or initial time t0)
c
if (mxnest .eq. 1) go to 99
c
lev = 1
40 if ((lev .eq. mxnest) .or. (lev .gt. lfine)) go to 60
mptr = lstart(lev)
50 itemp = node(store1,mptr)
node(store1,mptr) = node(store2,mptr)
node(store2,mptr) = itemp
rnode(timemult,mptr) = t0
mptr = node(levelptr,mptr)
if (mptr .ne. 0) go to 50
lev = lev + 1
go to 40
60 continue
c
c initial updating so can do conservation check. can do before
c bndry flux arrays set, since dont need them for this
c
do 65 level = 1, lfine-1
call update(lfine-level,nvar)
65 continue
c
c set up boundary flux conservation arrays
c
do 70 level = 1, lfine-1
call prepf(level+1,nvar,naux)
call prepc(level,nvar)
70 continue
c
if (.not. varRefTime) go to 99 ! keep consistent with gfixup and filval
do level = 1, lfine ! compute max speed for all grids at each level
spoh(level) = 0.d0 ! to keep track of max wave speed for all new grids
mptr = lstart(level)
75 continue
loc = node(store1,mptr)
locaux = node(storeaux,mptr)
nx = node(ndihi,mptr) - node(ndilo,mptr) + 1
ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1
mitot = nx + 2*nghost
mjtot = ny + 2*nghost
sp_over_h = get_max_speed(alloc(loc),mitot,mjtot,nvar,
& alloc(locaux),naux,nghost,
& hxposs(level),hyposs(level))
write(*,*)"max wave speed for level ",level," grid, ",mptr,
. " is", sp_over_h
spoh(level) = max(spoh(level),sp_over_h)
mptr = node(levelptr,mptr)
if (mptr .ne. 0) go to 75
end do
c
c check that level 1 has a time step that doesnt exceed cfl based on wave speed just calculated
c
dtc = possk(1)
if (dtc*spoh(1) .gt. cflv1) then
write(*,*)" coarse level time step too big: should reset"
write(*,*) " dtc ", dtc," gives cfl = ",dtc*spoh(1)
endif
c
c even if initial timestep was good could be too small
c try automatically resetting
c RANDY: do you like this option
c
possk(1) = cfl/spoh(1)
write(*,*)" automatically setting initial dt to ",possk(1)
c
c set new time step and refinement ratios in time
c at this initial time just setting initial conditions
c important to do this from coarser to finer levels to make sure
c subcycling lines up right.
c this sets ratios, and possk array.
c note that code below does not change coarsest level time step,
c this is set during timestepping or by user initially
c
do level = 2, lfine
dtc = possk(level-1)
dtf = cfl/spoh(level)
if (dtf .gt. dtc) then
kratio(level-1) = 1 ! cant have larger timestep than parent level
possk(level) = dtc ! cant have larger timestep than parent level
else
kratio(level-1) = ceiling(dtc/dtf)
possk(level) = possk(level-1)/kratio(level-1)
endif
write(6,*)" setting ref. ratio in time for level ",level," to ",
. kratio(level-1)
end do
c
99 continue
c
return
end