c
c
c =========================================================
subroutine rp1eu(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,fl,fr)
c =========================================================
c
c # Riemann solver for the 1D Euler equations
c # The waves are computed using the Roe approximation.
c
c # This is quite a bit slower than the Roe solver,
c # but may give more accurate solutions for some problems.
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, s the speeds,
c # fl and fr the positive and negative Godunov flux.
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 routine step1, rp is called with ql = qr = q.
c
c Author: Randall J. LeVeque
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 fr(1-mbc:maxmx+mbc, meqn)
dimension fl(1-mbc:maxmx+mbc, meqn)
c
c # local storage
c ---------------
parameter (max2 = 100002) !# assumes at most 100000 grid points with mbc=2
dimension delta(3)
dimension u(-1:max2),enth(-1:max2),a(-1:max2)
common /param/ gamma,gamma1
dimension sl(2), sr(2)
c
c # Riemann solver returns flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
c
c # Compute Roe-averaged quantities:
c
do 20 i=2-mbc,mx+mbc
rhsqrtl = dsqrt(qr(i-1,1))
rhsqrtr = dsqrt(ql(i,1))
pl = gamma1*(qr(i-1,3) - 0.5d0*(qr(i-1,2)**2)/qr(i-1,1))
pr = gamma1*(ql(i,3) - 0.5d0*(ql(i,2)**2)/ql(i,1))
rhsq2 = rhsqrtl + rhsqrtr
u(i) = (qr(i-1,2)/rhsqrtl + ql(i,2)/rhsqrtr) / rhsq2
enth(i) = (((qr(i-1,3)+pl)/rhsqrtl
& + (ql(i,3)+pr)/rhsqrtr)) / rhsq2
a2 = gamma1*(enth(i) - .5d0*u(i)**2)
a(i) = dsqrt(a2)
20 continue
c
c
do 30 i=2-mbc,mx+mbc
c
c # find a1 thru a3, the coefficients of the 3 eigenvectors:
c
delta(1) = ql(i,1) - qr(i-1,1)
delta(2) = ql(i,2) - qr(i-1,2)
delta(3) = ql(i,3) - qr(i-1,3)
a2 = gamma1/a(i)**2 * ((enth(i)-u(i)**2)*delta(1)
& + u(i)*delta(2) - delta(3))
a3 = (delta(2) + (a(i)-u(i))*delta(1) - a(i)*a2) / (2.d0*a(i))
a1 = delta(1) - a2 - a3
c
c # Compute the waves.
c
wave(i,1,1) = a1
wave(i,2,1) = a1*(u(i)-a(i))
wave(i,3,1) = a1*(enth(i) - u(i)*a(i))
s(i,1) = u(i)-a(i)
c
wave(i,1,2) = a2
wave(i,2,2) = a2*u(i)
wave(i,3,2) = a2*0.5d0*u(i)**2
s(i,2) = u(i)
c
wave(i,1,3) = a3
wave(i,2,3) = a3*(u(i)+a(i))
wave(i,3,3) = a3*(enth(i)+u(i)*a(i))
s(i,3) = u(i)+a(i)
30 continue
c
c # compute Godunov flux f0 at each interface.
c # Uses exact Riemann solver
c
do 200 i = 2-mbc, mx+mbc
c
rhol = qr(i-1,1)
rhor = ql(i ,1)
ul = qr(i-1,2)/qr(i-1,1)
ur = ql(i ,2)/ql(i ,1)
pl = gamma1*(qr(i-1,3) - 0.5d0*ul**2*qr(i-1,1))
pr = gamma1*(ql(i, 3) - 0.5d0*ur**2*ql(i ,1))
c
c # iterate to find pstar, ustar:
c
alpha = 1.
pstar = 0.5*(pl+pr)
wr = dsqrt(pr*rhor) * phi(pstar/pr)
wl = dsqrt(pl*rhol) * phi(pstar/pl)
c if (pl.eq.pr .and. rhol.eq.rhor) go to 60
c
40 do 50 iter=1,20
p1 = (ul-ur+pr/wr+pl/wl) / (1./wr + 1./wl)
pstar = dmax1(p1,1d-6)*alpha + (1.-alpha)*pstar
wr1 = wr
wl1 = wl
wr = dsqrt(pr*rhor) * phi(pstar/pr)
wl = dsqrt(pl*rhol) * phi(pstar/pl)
if (dmax1(abs(wr1-wr),dabs(wl1-wl)) .lt. 1d-6)
& go to 60
50 continue
c
c # nonconvergence:
alpha = alpha/2.
if (alpha .gt. 0.1) go to 40
c write(6,*) 'no convergence',wr1,wr,wl1,wl
wr = .5*(wr+wr1)
wl = .5*(wl+wl1)
c
60 continue
ustar = (pl-pr+wr*ur+wl*ul) / (wr+wl)
c
c # left wave:
c ============
c
if (pstar .gt. pl) then
c
c # shock:
sl(1) = ul - wl/rhol
sr(1) = sl(1)
rho1 = wl/(ustar-sl(1))
c
else
c
c # rarefaction:
cl = dsqrt(gamma*pl/rhol)
cstar = cl + 0.5*gamma1*(ul-ustar)
sl(1) = ul-cl
sr(1) = ustar-cstar
rho1 = (pstar/pl)**(1./gamma) * rhol
endif
c
c # right wave:
c =============
c
if (pstar .ge. pr) then
c
c # shock
sl(2) = ur + wr/rhor
sr(2) = sl(2)
rho2 = wr/(sl(2)-ustar)
c
else
c
c # rarefaction:
cr = dsqrt(gamma*pr/rhor)
cstar = cr + 0.5*gamma1*(ustar-ur)
sr(2) = ur+cr
sl(2) = ustar+cstar
rho2 = (pstar/pr)**(1./gamma)*rhor
endif
c
c # compute flux:
c ===============
c
c # compute state (rhos,us,ps) at x/t = 0:
c
if (sl(1).gt.0) then
rhos = rhol
us = ul
ps = pl
else if (sr(1).le.0. .and. ustar.ge. 0.) then
rhos = rho1
us = ustar
ps = pstar
else if (ustar.lt.0. .and. sl(2).ge. 0.) then
rhos = rho2
us = ustar
ps = pstar
else if (sr(2).lt.0) then
rhos = rhor
us = ur
ps = pr
else if (sl(1).le.0. .and. sr(1).ge.0.) then
c # transonic 1-rarefaction
us = (gamma1*ul + 2.*cl)/(gamma+1.)
e0 = pl/(rhol**gamma)
rhos = (us**2/(gamma*e0))**(1./gamma1)
ps = e0*rhos**gamma
else if (sl(2).le.0. .and. sr(2).ge.0.) then
c # transonic 3-rarefaction
us = (gamma1*ur - 2.*cr)/(gamma+1.)
e0 = pr/(rhor**gamma)
rhos = (us**2/(gamma*e0))**(1./gamma1)
ps = e0*rhos**gamma
endif
c
fl(i,1) = rhos*us
fl(i,2) = rhos*us**2 + ps
fl(i,3) = us*(gamma*ps/gamma1 + 0.5*rhos*us**2)
200 continue
c
do 220 m=1,3
do 220 i = 2-mbc, mx+mbc
fr(i,m) = -fl(i,m)
220 continue
c
return
end
c
c
double precision function phi(w)
implicit double precision (a-h,o-z)
common /param/ gamma,gamma1
c
sqg = dsqrt(gamma)
if (w .gt. 1.) then
phi = dsqrt(w*(gamma+1.)/2. + gamma1/2.)
else if (w .gt. 0.99999) then
phi = sqg
else if (w .gt. .999) then
phi = sqg + (2*gamma**2 - 3.*gamma + 1)
& *(w-1.) / (4.*sqg)
else
phi = gamma1*(1.-w) / (2.*sqg*(1.-w**(gamma1/(2.*gamma))))
endif
return
end
c