|
step1fw.f.html |
 |
|
Source file: step1fw.f
|
|
Directory: /Users/rjl/clawpack_src/clawpack_master/apps/fvmbook/chap16/vctraffic
|
|
Converted: Sat Mar 23 2024 at 11:04:24
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c
c ===================================================================
subroutine step1(maxmx,meqn,mwaves,mbc,mx,q,aux,dx,dt,
& method,mthlim,cfl,f,fwave,s,amdq,apdq,dtdx,rp1)
c ===================================================================
c
c # Take one time step, updating q.
c
c # MODIFIED to use fwave instead of wave.
c
c method(1) = 1 ==> Godunov method
c method(1) = 2 ==> Slope limiter method
c mthlim(p) controls what limiter is used in the pth family
c
c
c amdq, apdq, fwave, s, and f are used locally:
c
c amdq(1-mbc:maxmx+mbc, meqn) = left-going flux-differences
c apdq(1-mbc:maxmx+mbc, meqn) = right-going flux-differences
c e.g. amdq(i,m) = m'th component of A^- \Delta q from i'th Riemann
c problem (between cells i-1 and i).
c
c fwave(1-mbc:maxmx+mbc, meqn, mwaves) = waves from solution of
c Riemann problems,
c fwave(i,m,mw) = mth component of jump in f across
c wave in family mw in Riemann problem between
c states i-1 and i.
c
c s(1-mbc:maxmx+mbc, mwaves) = wave speeds,
c s(i,mw) = speed of wave in family mw in Riemann problem between
c states i-1 and i.
c
c f(1-mbc:maxmx+mbc, meqn) = correction fluxes for second order method
c f(i,m) = mth component of flux at left edge of ith cell
c --------------------------------------------------------------------
c
implicit double precision (a-h,o-z)
dimension q(meqn,1-mbc:maxmx+mbc)
dimension aux(*,1-mbc:maxmx+mbc)
dimension f(1-mbc:maxmx+mbc, meqn)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension fwave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension dtdx(1-mbc:maxmx+mbc)
dimension method(7),mthlim(mwaves)
logical limit
c
c # check if any limiters are used:
limit = .false.
do 5 mw=1,mwaves
if (mthlim(mw) .gt. 0) limit = .true.
5 continue
c
mcapa = method(6)
do 10 i=1-mbc,mx+mbc
if (mcapa.gt.0) then
if (aux(mcapa,i) .le. 0.d0) then
write(6,*) 'Error -- capa must be positive'
stop
endif
dtdx(i) = dt / (dx*aux(mcapa,i))
else
dtdx(i) = dt/dx
endif
10 continue
c
c
c
c # solve Riemann problem at each interface
c -----------------------------------------
c
call rp1(maxmx,meqn,mwaves,mbc,mx,q,q,aux,aux,fwave,s,amdq,apdq)
c
c # Modify q for Godunov update:
c # Note this may not correspond to a conservative flux-differencing
c # for equations not in conservation form. It is conservative if
c # amdq + apdq = f(q(i)) - f(q(i-1)).
c
do 40 i=1,mx+1
do 39 m=1,meqn
q(m,i) = q(m,i) - dtdx(i)*apdq(i,m)
q(m,i-1) = q(m,i-1) - dtdx(i-1)*amdq(i,m)
39 continue
40 continue
c
c # compute maximum wave speed:
cfl = 0.d0
do 50 mw=1,mwaves
do 45 i=1,mx+1
c # if s>0 use dtdx(i) to compute CFL,
c # if s<0 use dtdx(i-1) to compute CFL:
cfl = dmax1(cfl, dtdx(i)*s(i,mw), -dtdx(i-1)*s(i,mw))
45 continue
50 continue
c
if (method(2) .eq. 1) go to 900
c
c # compute correction fluxes for second order q_{xx} terms:
c ----------------------------------------------------------
c
do 100 m = 1, meqn
do 100 i = 1-mbc, mx+mbc
f(i,m) = 0.d0
100 continue
c
c # apply limiter to waves:
if (limit) then
call limiter(maxmx,meqn,mwaves,mbc,mx,fwave,s,mthlim)
endif
c
do 120 i=1,mx+1
do 120 m=1,meqn
do 110 mw=1,mwaves
dtdxave = 0.5d0 * (dtdx(i-1) + dtdx(i))
f(i,m) = f(i,m) + 0.5d0 * dsign(1.d0,s(i,mw))
& * (1.d0 - dabs(s(i,mw))*dtdxave) * fwave(i,m,mw)
c
110 continue
120 continue
c
c
140 continue
c
c # update q by differencing correction fluxes
c ============================================
c
c # (Note: Godunov update has already been performed above)
c
do 150 m=1,meqn
do 150 i=1,mx
q(m,i) = q(m,i) - dtdx(i) * (f(i+1,m) - f(i,m))
150 continue
c
900 continue
return
end