c
c =========================================================
subroutine rp1euznd(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the 1D ZND-Euler equations using Roe's
c # approximate Riemann solver.
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 # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
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)
c
c # local storage
c ---------------
parameter (max2 = 100002) !# assumes at most 100000 grid points with mbc=2
dimension u(-1:max2), enth(-1:max2),a(-1:max2)
dimension fr(-1:max2,4), fl(-1:max2,4), ihll(-1:max2)
logical efix, pfix, hll, roe, hllfix
common /param/ gamma,gamma1,q0
c
c define local arrays
c
dimension delta(4), Y(2,-1:max2)
c
data efix /.true./ !# use entropy fix for transonic rarefactions
data pfix /.true./ !# use Larrouturou's positivity fix for species
data hll /.true./ !# use HLL solver if unphysical values occur
data roe /.true./ !# use Roe solver
c
c # Riemann solver returns fluxes
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
c
c # Compute Roe-averaged quantities:
c
do 20 i=2-mbc,mx+mbc
c
pl = gamma1*(qr(i-1,4) - qr(i-1,2)*q0 -
& 0.5d0*qr(i-1,3)**2/(qr(i-1,1) + qr(i-1,2)))
pr = gamma1*(ql(i, 4) - ql(i, 2)*q0 -
& 0.5d0*ql(i, 3)**2/(ql(i, 1) + ql(i, 2)))
rhsqrtl = dsqrt(qr(i-1,1) + qr(i-1,2))
rhsqrtr = dsqrt(ql(i, 1) + ql(i, 2))
rhsq2 = rhsqrtl + rhsqrtr
u(i) = (qr(i-1,3)/rhsqrtl + ql(i,3)/rhsqrtr) / rhsq2
enth(i) = (((qr(i-1,4)+pl)/rhsqrtl
& + (ql(i ,4)+pr)/rhsqrtr)) / rhsq2
Y(1,i) = (qr(i-1,1)/rhsqrtl + ql(i,1)/rhsqrtr) / rhsq2
Y(2,i) = (qr(i-1,2)/rhsqrtl + ql(i,2)/rhsqrtr) / rhsq2
c # speed of sound
a2 = gamma1*(enth(i) - 0.5d0*u(i)**2 - Y(2,i)*q0)
a(i) = dsqrt(a2)
c
20 continue
c
do 30 i=2-mbc,mx+mbc
c
c # find a1 thru a3, the coefficients of the 4 eigenvectors:
c
do k = 1, 4
delta(k) = ql(i,k) - qr(i-1,k)
enddo
drho = delta(1) + delta(2)
c
a2 = gamma1/a(i)**2 * (drho*0.5d0*u(i)**2 - delta(2)*q0
& - u(i)*delta(3) + delta(4))
a3 = 0.5d0*( a2 - ( u(i)*drho - delta(3) )/a(i) )
a1 = a2 - a3
c
c # Compute the waves.
c
c # 1-wave
wave(i,1,1) = a1*Y(1,i)
wave(i,2,1) = a1*Y(2,i)
wave(i,3,1) = a1*(u(i)-a(i))
wave(i,4,1) = a1*(enth(i) - u(i)*a(i))
s(i,1) = u(i)-a(i)
c
c # 2-wave and 3-wave
wave(i,1,2) = delta(1) - Y(1,i)*a2
wave(i,2,2) = delta(2) - Y(2,i)*a2
wave(i,3,2) = u(i)*(drho - a2)
wave(i,4,2) = 0.5d0*u(i)**2*(drho - a2) +
& q0*(delta(2) - Y(2,i)*a2)
s(i,2) = u(i)
c
c # 4-wave
wave(i,1,3) = a3*Y(1,i)
wave(i,2,3) = a3*Y(2,i)
wave(i,3,3) = a3*(u(i)+a(i))
wave(i,4,3) = a3*(enth(i)+u(i)*a(i))
s(i,3) = u(i)+a(i)
c
c aa2 = 1.d0
c if (a2.ne.0.d0) aa2 = drho/a2-1.d0
c write (6,*) i,a1,aa2,a3
c
ihll(i) = 0
if (dabs(a3).gt.1.d0) ihll(i) = 1
c
30 continue
c
call flx1(maxmx,meqn,mbc,mx,qr,maux,auxr,apdq)
call flx1(maxmx,meqn,mbc,mx,ql,maux,auxl,amdq)
c
do 35 i = 1-mbc, mx+mbc
do 35 m=1,meqn
fl(i,m) = amdq(i,m)
fr(i,m) = apdq(i,m)
35 continue
c
if (.not.roe) go to 900
c
c # compute Godunov flux f0:
c --------------------------
c
if (efix) go to 110
c
c # no entropy fix
c ----------------
c
c # amdq = SUM s*wave over left-going waves
c
do 100 m=1,meqn
do 100 i=2-mbc, mx+mbc
amdq(i,m) = 0.d0
do 90 mw=1,mwaves
if (s(i,mw) .lt. 0.d0)
& amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw)
90 continue
100 continue
go to 900
110 continue
c
c # With entropy fix
c ------------------
c
c # compute flux differences amdq.
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
rk1 = qr(i-1,1)
rk2 = qr(i-1,2)
rhou = qr(i-1,3)
rhoE = qr(i-1,4)
rho = rk1 + rk2
p = gamma1*(rhoE - rk2*q0 - 0.5d0*rhou**2/rho)
if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200
c = dsqrt(gamma*p/rho)
s0 = rhou/rho - c !# u-c in left state (cell i-1)
* write(6,*) 'left state 0', a(i), c, T
c
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,meqn
amdq(i,m) = 0.d0
60 continue
go to 200
endif
c
rk1 = rk1 + wave(i,1,1)
rk2 = rk2 + wave(i,2,1)
rhou = rhou + wave(i,3,1)
rhoE = rhoE + wave(i,4,1)
rho = rk1 + rk2
p = gamma1*(rhoE - rk2*q0 - 0.5d0*rhou**2/rho)
if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200
c = dsqrt(gamma*p/rho)
s1 = rhou/rho - c !# u-c to right of 1-wave
* write(6,*) 'left state 1', a(i), c, T
c
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,meqn
amdq(i,m) = sfract*wave(i,m,1)
120 continue
c
c # check 2-wave:
c ---------------
c
if (s(i,2) .ge. 0.d0) go to 200 !# 2-wave is rightgoing
do 140 m=1,meqn
amdq(i,m) = amdq(i,m) + s(i,2)*wave(i,m,2)
140 continue
c
c # check 3-wave:
c ---------------
c
rk1 = ql(i,1)
rk2 = ql(i,2)
rhou = ql(i,3)
rhoE = ql(i,4)
rho = rk1 + rk2
p = gamma1*(rhoE - rk2*q0 - 0.5d0*rhou**2/rho)
if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200
c = dsqrt(gamma*p/rho)
s3 = rhou/rho + c !# u+c in right state (cell i)
* write(6,*) 'right state 1', a(i), c, T
c
rk1 = rk1 - wave(i,1,3)
rk2 = rk2 - wave(i,2,3)
rhou = rhou - wave(i,3,3)
rhoE = rhoE - wave(i,4,3)
rho = rk1 + rk2
p = gamma1*(rhoE - rk2*q0 - 0.5d0*rhou**2/rho)
if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200
c = dsqrt(gamma*p/rho)
s2 = rhou/rho + c !# u+c to left of 3-wave
* write(6,*) 'right state 0', a(i), c, T
c
if (s2 .lt. 0.d0 .and. s3.gt.0.d0) then
c # transonic rarefaction in the 3-wave
sfract = s2 * (s3-s(i,3)) / (s3-s2)
else if (s(i,3) .lt. 0.d0) then
c # 3-wave is leftgoing
sfract = s(i,3)
else
c # 3-wave is rightgoing
go to 200
endif
c
do 160 m=1,meqn
amdq(i,m) = amdq(i,m) + sfract*wave(i,m,3)
160 continue
200 continue
c
900 continue
c
if (hll) then
do 250 i = 2-mbc, mx+mbc
hllfix = .false.
if (.not.roe) hllfix = .true.
c
rho1l = qr(i-1,1) + wave(i,1,1)
rho2l = qr(i-1,2) + wave(i,2,1)
rhoul = qr(i-1,3) + wave(i,3,1)
El = qr(i-1,4) + wave(i,4,1)
pl = gamma1*(El - rho2l*q0 - 0.5d0*rhoul**2/(rho1l+rho2l))
if (rho1l+rho2l.le.0.d0.or.pl.le.0.d0)
& hllfix = .true.
c
rho1r = ql(i,1) - wave(i,1,3)
rho2r = ql(i,2) - wave(i,2,3)
rhour = ql(i,3) - wave(i,3,3)
Er = ql(i,4) - wave(i,4,3)
pr = gamma1*(Er - rho2r*q0 - 0.5d0*rhour**2/(rho1r+rho2r))
if (rho1r+rho2r.le.0.d0.or.pr.le.0.d0.or.ihll(i).eq.1)
& hllfix = .true.
c
if (hllfix) then
c if (roe) write (6,*) 'Switching to HLL in',i
c
rl = qr(i-1,1) + qr(i-1,2)
ul = qr(i-1,3)/rl
pl = gamma1*(qr(i-1,4) - qr(i-1,2)*q0 -
& 0.5d0*qr(i-1,3)**2/rl)
al = dsqrt(gamma*pl/rl)
c
rr = ql(i ,1) + ql(i ,2)
ur = ql(i ,3)/rr
pr = gamma1*(ql(i ,4) - ql(i ,2)*q0 -
& 0.5d0*ql(i ,3)**2/rr)
ar = dsqrt(gamma*pr/rr)
c
sl = dmin1(ul-al,ur-ar)
sr = dmax1(ul+al,ur+ar)
c
do m=1,meqn
if (sl.ge.0.d0) fg = fr(i-1,m)
if (sr.le.0.d0) fg = fl(i,m)
if (sl.lt.0.d0.and.sr.gt.0.d0)
& fg = (sr*fr(i-1,m) - sl*fl(i,m) +
& sl*sr*(ql(i,m)-qr(i-1,m)))/ (sr-sl)
amdq(i,m) = fg-fr(i-1,m)
enddo
s(i,1) = sl
s(i,2) = 0.d0
s(i,3) = sr
endif
250 continue
endif
c
if (pfix) then
do 300 i=2-mbc,mx+mbc
amdr = amdq(i,1)+amdq(i,2)
rhol = qr(i-1,1)+qr(i-1,2)
rhor = ql(i ,1)+ql(i ,2)
do 300 m=1,2
if (qr(i-1,3)+amdr.gt.0.d0) then
Z = qr(i-1,m)/rhol
else
Z = ql(i ,m)/rhor
endif
amdq(i,m) = Z*amdr + (Z-qr(i-1,m)/rhol)*qr(i-1,3)
300 continue
endif
c
do 400 i = 2-mbc, mx+mbc
do 400 m=1,meqn
amdq(i,m) = fr(i-1,m) + amdq(i,m)
apdq(i,m) = -amdq(i,m)
400 continue
c
return
end
c