|
|
rp1eu_geneos_wb.f.html
| |
|
Source file: rp1eu_geneos_wb.f
|
|
Directory: /var/www/html/clawpack/links/shockvacuum10/shockvacuum
|
|
Converted: Mon Aug 16 2010 at 14:13:31
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c
c
c =========================================================
subroutine rp1(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
& fwave,s,amdq,apdq)
c =========================================================
c
c f-wave solver for the Riemann problems for the 1D gas dynamics
c equations with a general equation of state (EOS) using a variant of
c the "Suliciu relaxation solver" described in Sec. 2.4.5 of
c Francois Bouchut's book "Nonlinear Stability of Finite
c Volume methods for Hyperbolic Conservation Laws and Well-Balanced
c Schemes for Sources" (Birkh\"auser 2004). This solver works well
c near vacuum states. The EOS should be specified in the subroutine
c pres_psi, which returns a vector pres containing the pressure
c and cs, the sound speed, in each cell.
c
c Note that in this Riemann solver it is assumed that ql=qr and only
c ql is used.
c
c The solver described in Bouchut has been modified by RJL and
c Jeremiah Murphy to work in f-wave form so that source terms can also
c be incorporated. [cite paper, to appear].
c
c The subroutine pres_psi should also return
c the source term psi*dx at the left edge of each cell: psidx(i) is
c used to modify dp at the interface between cell i-1 and i, where dp
c is the pressure difference between these cells.
c
c -------------------------------------------------------
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 fwave(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,*)
common /comxt/ dtcom,dxcom,tcom
double precision npoly, nstar
common /ceos/ grav, gamma, gamma1, sK, npoly, nstar, Crho, ieos
c
c # local storage
c ---------------
parameter (max2 = 10002) !# assumes at most 10000 grid points with mbc=2
dimension pres(-1:max2)
dimension cs(-1:max2)
dimension psidx(-1:max2)
if (mx+mbc .gt. max2) then
write(6,*) '*** Increase max2 in rp1'
stop
endif
c
c
c # compute the pressure pres, sound speed cs, and source term psi*dx:
call pres_psi(ql,pres,cs,psidx,mx,maxmx,mbc,meqn)
write(60,601) 1200,pres(1200),cs(1200),psidx(1200)
write(60,601) 1201,pres(1201),cs(1201),psidx(1201)
601 format(i5,3e16.6)
c
c
do i=2-mbc,mx+mbc
c find primative values to the left and right of the contact
c discontinuity.
pl = pres(i-1)
pr = pres(i)
rhol = ql(i-1,1)
rhor = ql(i,1)
ul = ql(i-1,2)/rhol
ur = ql(i,2)/rhor
alpha = (gamma+1.d0)*0.5d0
c # Jump in pressure is modified by source term correction for
c # well-balanced method:
dp = pr - pl - psidx(i)
c # Useful for debugging:
c write(40,401) i,rhol,rhor,ul,ur
c write(40,401) i,pl,pr,psidx(i),dp
c write(40,*) ' '
c 401 format(i5,4e15.6)
c # Formulas (2.145, 2.146) from Bouchut:
c # ccl and ccr are the $\tilde c_l$ and $\tilde c_r$ values there,
c # which are sound speeds multiplied by rho values
if (dp.gt.0.d0) then
cl = cs(i-1) + alpha*max((dp/rhor/cs(i) +ul-ur),0.d0)
cr = cs(i) + alpha*max((-dp/rhol/cl + ul-ur),0.d0)
else
cr = cs(i) + alpha*max((-dp/rhol/cs(i-1) + ul-ur),0.d0)
cl = cs(i-1) + alpha*max((dp/rhor/cr + ul-ur),0.d0)
endif
ccl = rhol*cl
ccr = rhor*cr
c # Internal energy:
eccl = ql(i-1,3)-0.5d0*rhol*ul**2
eccr = ql(i,3)-0.5d0*rhor*ur**2
c
c # Enthalpies, not used:
c hl = (eccl+pl)/rhol
c hr = (eccr+pr)/rhor
c # Equation (2.135) from Bouchut:
ustar = (ul*ccl + ur*ccr - dp)/(ccr + ccl)
pstar = (pl*ccr + pr*ccl + (ul - ur)*ccr*ccl)
& /(ccr+ccl)
rstarl = 1.d0/
& (1.d0/rhol + (ccr*(ur-ul) - dp)/ccl/(ccl+ccr))
rstarr = 1.d0/
& (1.d0/rhor + (ccl*(ur-ul) + dp)/ccr/(ccl+ccr))
c
c # From Jeremiah Murphy:
c pstarl = pstar - ccl*psidx/(ccl+ccr)
c pstarr = pstar + ccr*psidx/(ccl+ccr)
c # Equivalent expressions:
pstarl = pl + (ccl*dp + (ul-ur)*ccr*ccl)/(ccl+ccr)
pstarr = pr - (ccr*dp - (ul-ur)*ccr*ccl)/(ccl+ccr)
estarl = rstarl*(eccl/rhol + (pstarl**2 - pl**2)/2.d0/ccl**2)
estarr = rstarr*(eccr/rhor + (pstarr**2 - pr**2)/2.d0/ccr**2)
c
c # Compute the waves.
c
c f(qstarl) - f(ql)
fwave(i,1,1) = rstarl*ustar - rhol*ul
fwave(i,2,1) = rstarl*ustar**2 + pstarl - rhol*ul**2 - pl
fwave(i,3,1) = ustar*(estarl + 0.5d0*rstarl*ustar**2 + pstarl)
& - ul*(eccl + 0.5d0*rhol*ul**2 + pl)
s(i,1) = ul - cl
c
c f(qstarr) - f(qstarl)
fwave(i,1,2) = (rstarr - rstarl)*ustar
fwave(i,2,2) = (rstarr - rstarl)*ustar**2
fwave(i,3,2) = ustar*(estarr - estarl
& + 0.5d0*(rstarr - rstarl)*ustar**2)
s(i,2) = ustar
c
c f(qr) - f(qstarr)
fwave(i,1,3) = rhor*ur - rstarr*ustar
fwave(i,2,3) = rhor*ur**2 + pr - rstarr*ustar**2 - pstarr
fwave(i,3,3) = ur*(eccr + 0.5d0*rhor*ur**2 + pr)
& - ustar*(estarr + 0.5d0*rstarr*ustar**2 + pstarr)
s(i,3) = ur + cr
c write(40,401) i,fwave(i,2,1),fwave(i,2,2),fwave(i,2,3)
c write(40,*) ' '
enddo
c
c
c
c # amdq = SUM s*wave over left-going waves
c # apdq = SUM s*wave over right-going waves
c
stol = 0.d0
do m=1,meqn
do i=2-mbc, mx+mbc
amdq(i,m) = 0.d0
apdq(i,m) = 0.d0
do mw=1,mwaves
if (s(i,mw) .lt. stol) then
amdq(i,m) = amdq(i,m) + fwave(i,m,mw)
else if (s(i,mw) .gt. stol) then
apdq(i,m) = apdq(i,m) + fwave(i,m,mw)
else
amdq(i,m) = amdq(i,m) + 0.5d0*fwave(i,m,mw)
apdq(i,m) = apdq(i,m) + 0.5d0*fwave(i,m,mw)
endif
enddo
enddo
enddo
return
end