|
rp1sw.f.html |
 |
|
Source file: rp1sw.f
|
|
Directory: /home/rjl/git/rjleveque/clawpack-4.x/book/chap21/radialdam/1drad
|
|
Converted: Tue Jul 26 2011 at 12:59:05
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c =========================================================
subroutine rp1(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
& wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the 1D shallow water equations
c # (h)_t + (u h)_x = 0
c # (uh)_t + ( uuh + .5*gh^2 )_x = 0
c # using Roe's approximate Riemann solver with entropy fix for
c # transonic rarefractions.
c
c # On input, ql contains the state vector at the left edge of each cell
c # qr contains the state vector at the right edge of each cell
c # On output, wave contains the waves,
c # s the speeds,
c # amdq the left-going flux difference A^- \Delta q
c # apdq the right-going flux difference A^+ \Delta q
c
c # Note that the i'th Riemann problem has left state qr(i-1,:)
c # and right state ql(i,:)
c # From the basic clawpack routine step1, rp is called with ql = qr = q.
c
c Here meqn=mwaves=2 should be passed from the calling routine
c
implicit double precision (a-h,o-z)
dimension ql(1-mbc:maxmx+mbc, meqn)
dimension qr(1-mbc:maxmx+mbc, meqn)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension auxl(1-mbc:maxmx+mbc, *)
dimension auxr(1-mbc:maxmx+mbc, *)
c
c # local storage
c ---------------
dimension delta(2)
logical efix
common /comrp/ grav
data efix /.true./ !# use entropy fix for transonic rarefactions
c
do 30 i=2-mbc,mx+mbc
c
c
c # compute Roe-averaged quantities:
ubar = (qr(i-1,2)/dsqrt(qr(i-1,1)) + ql(i,2)/dsqrt(ql(i,1)))/
. ( dsqrt(qr(i-1,1)) + dsqrt(ql(i,1)) )
cbar=dsqrt(0.5d0*grav*(qr(i-1,1) + ql(i,1)))
c # delta(1)=h(i)-h(i-1) and delta(2)=hu(i)-hu(i-1)
delta(1) = ql(i,1) - qr(i-1,1)
delta(2) = ql(i,2) - qr(i-1,2)
c # compute coeffs in the evector expansion of delta(1),delta(2)
a1 = 0.5d0*(-delta(2) + (ubar + cbar) * delta(1))/cbar
a2 = 0.5d0*( delta(2) - (ubar - cbar) * delta(1))/cbar
c # finally, compute the waves.
wave(i,1,1) = a1
wave(i,2,1) = a1*(ubar - cbar)
s(i,1) = ubar - cbar
wave(i,1,2) = a2
wave(i,2,2) = a2*(ubar + cbar)
s(i,2) = ubar + cbar
30 continue
c # compute Godunov flux f0:
c --------------------------
if (efix) go to 110
c # no entropy fix
c ----------------------------------------------
c # amdq = SUM s*wave over left-going waves
c # apdq = SUM s*wave over right-going waves
do 100 m=1,2
do 100 i=2-mbc, mx+mbc
amdq(i,m) = 0.d0
apdq(i,m) = 0.d0
do 90 mw=1,mwaves
if (s(i,mw) .lt. 0.d0) then
amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw)
else
apdq(i,m) = apdq(i,m) + s(i,mw)*wave(i,m,mw)
endif
90 continue
100 continue
go to 900
c
c -----------------------------------------------
c
110 continue
c # With entropy fix
c ------------------
c
c # compute flux differences amdq and apdq.
c # First compute amdq as sum of s*wave for left going waves.
c # Incorporate entropy fix by adding a modified fraction of wave
c # if s should change sign.
c
do 200 i=2-mbc,mx+mbc
c ------------------------------------------------------
c # check 1-wave:
c ---------------
c
c # u-c in left state (cell i-1)
s0 = qr(i-1,2)/qr(i-1,1) - dsqrt(grav*qr(i-1,1))
c # check for fully supersonic case:
if (s0.ge.0.d0 .and. s(i,1).gt.0.d0) then
c # everything is right-going
do 60 m=1,2
amdq(i,m) = 0.d0
60 continue
go to 200
endif
c
c # u-c to right of 1-wave
hr1 = qr(i-1,1) + wave(i,1,1)
uhr1 = qr(i-1,2) + wave(i,2,1)
s1 = uhr1/hr1 - dsqrt(grav*hr1)
if (s0.lt.0.d0 .and. s1.gt.0.d0) then
c # transonic rarefaction in the 1-wave
sfract = s0 * (s1-s(i,1)) / (s1-s0)
else if (s(i,1) .lt. 0.d0) then
c # 1-wave is leftgoing
sfract = s(i,1)
else
c # 1-wave is rightgoing
sfract = 0.d0 !# this shouldn't happen since s0 < 0
endif
do 120 m=1,2
amdq(i,m) = sfract*wave(i,m,1)
120 continue
c -------------------------------------------------------
c # check 2-wave:
c ---------------
c # u+c in right state (cell i)
s3 = ql(i,2)/ql(i,1) + dsqrt(grav*ql(i,1))
c # u+c to left of 2-wave
hl2 = ql(i,1) - wave(i,1,2)
uhl2 = ql(i,2) - wave(i,2,2)
s2 = uhl2/hl2 + dsqrt(grav*hl2)
if (s2 .lt. 0.d0 .and. s3.gt.0.d0) then
c # transonic rarefaction in the 2-wave
sfract = s2 * (s3-s(i,2)) / (s3-s2)
else if (s(i,2) .lt. 0.d0) then
c # 2-wave is leftgoing
sfract = s(i,2)
else
c # 2-wave is rightgoing
go to 200
endif
c
do 160 m=1,2
amdq(i,m) = amdq(i,m) + sfract*wave(i,m,2)
160 continue
200 continue
c
c
c # compute the rightgoing flux differences:
c # df = SUM s*wave is the total flux difference and apdq = df - amdq
c
do 220 m=1,2
do 220 i = 2-mbc, mx+mbc
df = 0.d0
do 210 mw=1,mwaves
df = df + s(i,mw)*wave(i,m,mw)
210 continue
apdq(i,m) = df - amdq(i,m)
220 continue
900 continue
return
end