c
c =========================================================
subroutine rp1eurhok(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the thermally perfect 1D multi-component
c # Euler equations using a Roe-type approximate Riemann solver.
c # Scheme is blended with HLL for robustness.
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)
dimension auxl(1-mbc:maxmx+mbc, maux)
dimension auxr(1-mbc:maxmx+mbc, maux)
c
c # local storage
c ---------------
parameter (max2 = 10002) !# assumes at most 10000 grid points with mbc=2
dimension u(-1:max2),enth(-1:max2),a(-1:max2),smax(-1:max2)
dimension g1a2(-1:max2)
logical efix, pfix, hll, roe, hllfix
c
c define local arrays
c
include "ck.i"
c
dimension delta(LeNsp+2)
dimension rkl(LeNsp), rkr(LeNsp)
dimension hkl(LeNsp), hkr(LeNsp)
dimension Y(LeNsp,-1:max2), pk(LeNsp,-1:max2)
dimension fl(-1:max2,LeNsp+3), fr(-1:max2,LeNsp+3)
c
data efix /.true./ !# use entropy fix
data pfix /.true./ !# use Larrouturou's positivity fix for species
data hll /.true./ !# use HLL instead of Roe solver, if unphysical values occur
data roe /.true./ !# turn off Roe solver when debugging HLL
c
c # Riemann solver returns fluxes
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
c
mu = Nsp+1
mE = Nsp+2
mT = Nsp+3
c
c # Compute Roe-averaged quantities:
c
do 20 i=2-mbc,mx+mbc
rhol = 0.d0
rhor = 0.d0
do k = 1, Nsp
rkl(k) = qr(i-1,k)
rkr(k) = ql(i ,k)
rhol = rhol + rkl(k)
rhor = rhor + rkr(k)
enddo
if( rhol.le.1.d-10 ) then
write(6,*) 'negative total density, left', rhol
stop
endif
if( rhor.le.1.d-10 ) then
write(6,*) 'negative total density, right', rhor
stop
endif
c
c # compute left/right rho/W and rho*Cp
c
rhoWl = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rhoWl = rhoWl + rkl(k)/Wk(k)
rhoWr = rhoWr + rkr(k)/Wk(k)
enddo
c
c # calculate left/right Temperatures
c
rhoel = qr(i-1,mE)-0.5d0*qr(i-1,mu)**2/rhol
call SolveTrhok(qr(i-1,mT),rhoel,rhoWl,rkl,Nsp,ifail)
rhoer = ql(i ,mE)-0.5d0*ql(i ,mu)**2/rhor
call SolveTrhok(ql(i ,mT),rhoer,rhoWr,rkr,Nsp,ifail)
c
Tl = qr(i-1,mT)
Tr = ql(i ,mT)
pl = rhoWl*RU*Tl
pr = rhoWr*RU*Tr
c
c # compute quantities for rho-average
c
rhsqrtl = dsqrt(rhol)
rhsqrtr = dsqrt(rhor)
rhsq2 = rhsqrtl + rhsqrtr
c
c # find rho-averaged specific velocity and enthalpy
c
u(i) = (qr(i-1,mu)/rhsqrtl +
& ql(i ,mu)/rhsqrtr) / rhsq2
enth(i) = (((qr(i-1,mE)+pl)/rhsqrtl
& + (ql(i ,mE)+pr)/rhsqrtr)) / rhsq2
c
c # compute rho-averages for T, cp, and W
c
T = (Tl * rhsqrtl + Tr * rhsqrtr) / rhsq2
W = rhsq2 / (rhoWl/rhsqrtl + rhoWr/rhsqrtr)
c
c # evaluate left/right entropies and mean cp
c
call tabintp( Tl, hkl, hms, Nsp )
call tabintp( Tr, hkr, hms, Nsp )
do k = 1, Nsp
Y(k,i) = (rkl(k)/rhsqrtl + rkr(k)/rhsqrtr) / rhsq2
enddo
Cp = Cpmix( Tl, Tr, hkl, hkr, Y(1,i) )
gamma1 = RU / ( W*Cp - RU )
gamma = gamma1 + 1.d0
c
c # find rho-averaged specific enthalpies,
c # compute rho-averaged mass fractions and
c # compute partial pressure derivatives
c
tmp = gamma * RU * T / gamma1
* ht = 0.d0
do k = 1, Nsp
hk = (hkl(k)*rhsqrtl + hkr(k)*rhsqrtr) / rhsq2
* ht = ht + Y(k,i)*(hkl(k)*rhsqrtl + hkr(k)*rhsqrtr) / rhsq2
pk(k,i) = 0.5d0*u(i)**2 - hk + tmp / Wk(k)
enddo
c
* write (6,4) qr(i-1,mE)+pl, ql(i,mE)+pr,
* & ht+0.5d0*u(i)**2, enth(i), ht+0.5d0*u(i)**2-enth(i)
* 4 format(e16.8,e16.8,e16.8,e16.8,e24.16)
c
c # compute speed of sound
c
a2 = enth(i)-u(i)**2
do k = 1, Nsp
a2 = a2 + Y(k,i) * pk(k,i)
enddo
g1a2(i) = 1.d0 / a2
a(i) = dsqrt(gamma1*a2)
c
rhoCpl = avgtabip( Tl, rkl, cpk, Nsp )
rhoCpr = avgtabip( Tr, rkr, cpk, Nsp )
gammal = RU / ( rhoCpl/rhoWl - RU ) + 1.d0
gammar = RU / ( rhoCpr/rhoWr - RU ) + 1.d0
ul = qr(i-1,mu)/rhol
ur = ql(i ,mu)/rhor
al = dsqrt(gammal*pl/rhol)
ar = dsqrt(gammar*pr/rhor)
smax(i) = dmax1(dmax1(dabs(ur-ar-(ul-al)),dabs(ur-ul)),
& dabs(ur+ar-(ul+al)))
c
20 continue
c
c
do 30 i=2-mbc,mx+mbc
c
c # find a1 thru a3, the coefficients of the mE eigenvectors:
c
dpdr = 0.d0
dpY = 0.d0
drho = 0.d0
do k = 1, Nsp
delta(k) = ql(i,k) - qr(i-1,k)
drho = drho + delta(k)
dpdr = dpdr + pk(k,i) * delta(k)
dpY = dpY + pk(k,i) * Y(k,i)
enddo
delta(mu) = ql(i,mu) - qr(i-1,mu)
delta(mE) = ql(i,mE) - qr(i-1,mE)
c
a2 = g1a2(i)*(dpdr - u(i)*delta(mu) + delta(mE))
a3 = 0.5d0*( a2 - ( u(i)*drho - delta(mu) )/a(i) )
a1 = a2 - a3
c
c
c # Compute the waves.
c # Note that the 1+k-waves, for 1 .le. k .le. Nsp travel at
c # the same speed and are lumped together in wave(.,.,2).
c # The 3-wave is then stored in wave(.,.,3).
c
do k = 1, Nsp
c # 1-wave
wave(i,k,1) = a1*Y(k,i)
c # 2-wave
wave(i,k,2) = delta(k) - Y(k,i)*a2
c # 3-wave
wave(i,k,3) = a3*Y(k,i)
enddo
c # 1-wave
wave(i,mu,1) = a1*(u(i)-a(i))
wave(i,mE,1) = a1*(enth(i) - u(i)*a(i))
wave(i,mT,1) = 0.d0
s(i,1) = u(i)-a(i)
c
c # 2-wave
wave(i,mu,2) = u(i)*(drho - a2)
wave(i,mE,2) = u(i)**2*(drho - a2) - dpdr + dpY*a2
wave(i,mT,2) = 0.d0
s(i,2) = u(i)
c
c # 3-wave
wave(i,mu,3) = a3*(u(i)+a(i))
wave(i,mE,3) = a3*(enth(i)+u(i)*a(i))
wave(i,mT,3) = 0.d0
s(i,3) = u(i)+a(i)
c
30 continue
c
c
c # compute fluxes as
c # F(Ur,Ul) = 0.5*( f(Ur)+f(Ul) - |A|(Ur-Ul) )
c --------------------------
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 (roe) then
do 40 i = 2-mbc, mx+mbc
do 40 m=1,meqn
amdq(i,m) = 0.5d0*(fl(i,m)+fr(i-1,m))
40 continue
c
do 50 i = 2-mbc, mx+mbc
do 50 m=1,meqn
sw = 0.d0
do 60 mw=1,mwaves
sl = dabs(s(i,mw))
if (efix.and.dabs(s(i,mw)).lt.smax(i).and.mw.ne.2)
& sl = s(i,mw)**2/(2.d0*smax(i))+0.5d0*smax(i)
sw = sw + sl*wave(i,m,mw)
60 continue
amdq(i,m) = amdq(i,m) - 0.5d0*sw
50 continue
endif
c
if (hll) then
do 55 i = 2-mbc, mx+mbc
c # set this to hllfix = .true. when debugging HLL
hllfix = .false.
if (.not.roe) hllfix = .true.
c
rhol = 0.d0
rhoWl = 0.d0
do k = 1, Nsp
rkl(k) = qr(i-1,k) + wave(i,k,1)
rhol = rhol + rkl(k)
rhoWl = rhoWl + rkl(k)/Wk(k)
enddo
rhoul = qr(i-1,mu) + wave(i,mu,1)
ul = rhoul/rhol
rhoEl = qr(i-1,mE) + wave(i,mE,1)
Tl = qr(i-1,mT)
rhoel = rhoEl - 0.5d0*rhoul**2/rhol
call SolveTrhok( Tl, rhoel, rhoWl, rkl, Nsp, ifail)
rhoCpl = avgtabip( Tl, rkl, cpk, Nsp )
gammal = RU / ( rhoCpl/rhoWl - RU ) + 1.d0
pl = rhoWl*RU*Tl
al = dsqrt(gammal*pl/rhol)
if (rhol.le.0.d0.or.pl.le.0.d0) hllfix = .true.
c
rhor = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rkr(k) = ql(i ,k) - wave(i,k,3)
rhor = rhor + rkr(k)
rhoWr = rhoWr + rkr(k)/Wk(k)
enddo
rhour = ql(i ,mu) - wave(i,mu,3)
ur = rhoul/rhol
rhoEr = ql(i ,mE) - wave(i,mE,3)
Tr = ql(i ,mT)
rhoer = rhoEr - 0.5d0*rhour**2/rhor
call SolveTrhok( Tr, rhoer, rhoWr, rkr, Nsp, ifail)
rhoCpr = avgtabip( Tr, rkr, cpk, Nsp )
gammar = RU / ( rhoCpr/rhoWr - RU ) + 1.d0
pr = rhoWr*RU*Tr
ar = dsqrt(gammar*pr/rhor)
if (rhor.le.0.d0.or.pr.le.0.d0) hllfix = .true.
c
if (hllfix) then
c if (roe) write (6,*) 'Switching to HLL in',i
c
rhol = 0.d0
rhoWl = 0.d0
do k = 1, Nsp
rkl(k) = qr(i-1,k)
rhol = rhol + qr(i-1,k)
rhoWl = rhoWl + qr(i-1,k)/Wk(k)
enddo
ul = qr(i-1,mu)/rhol
Tl = qr(i-1,mT)
pl = rhoWl*RU*Tl
rhoCpl = avgtabip( Tl, rkl, cpk, Nsp )
gammal = RU / ( rhoCpl/rhoWl - RU ) + 1.d0
al = dsqrt(gammal*pl/rhol)
c
rhor = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rkr(k) = ql(i ,k)
rhor = rhor + ql(i ,k)
rhoWr = rhoWr + ql(i ,k)/Wk(k)
enddo
ur = ql(i ,mu)/rhor
Tr = ql(i ,mT)
pr = rhoWr*RU*Tr
rhoCpr = avgtabip( Tr, rkr, cpk, Nsp )
gammar = RU / ( rhoCpr/rhoWr - RU ) + 1.d0
ar = dsqrt(gammar*pr/rhor)
c
sl = dmin1(ul-al,ur-ar)
sr = dmax1(ul+al,ur+ar)
c
do m=1,meqn
if (sl.ge.0.d0) amdq(i,m) = fr(i-1,m)
if (sr.le.0.d0) amdq(i,m) = fl(i,m)
if (sl.lt.0.d0.and.sr.gt.0.d0)
& amdq(i,m) = (sr*fr(i-1,m) - sl*fl(i,m) +
& sl*sr*(ql(i,m)-qr(i-1,m)))/ (sr-sl)
enddo
amdq(i,mT) = 0.d0
s(i,1) = sl
s(i,2) = 0.d0
s(i,3) = sr
endif
55 continue
endif
c
if (pfix) then
do 70 i=2-mbc,mx+mbc
amdr = 0.d0
rhol = 0.d0
rhor = 0.d0
do k = 1, Nsp
amdr = amdr + amdq(i,k)
rhol = rhol + qr(i-1,k)
rhor = rhor + ql(i ,k)
enddo
do k=1, Nsp
if (amdr.gt.0.d0) then
Z = qr(i-1,k)/rhol
else
Z = ql(i ,k)/rhor
endif
amdq(i,k) = Z*amdr
enddo
70 continue
endif
c
do 80 i = 2-mbc, mx+mbc
do 80 m=1,meqn
apdq(i,m) = -amdq(i,m)
80 continue
c
return
end
c
c
c ***********************************************************
c
double precision function Cpmix( Tl, Tr, hl, hr, Y )
implicit double precision(a-h,o-z)
include "ck.i"
c
dimension Y(*)
dimension hl(*), hr(*)
data Tol /1.d-6/
c
if( dabs(Tr-Tl).gt.Tol ) then
Cp = 0.d0
do k = 1, Nsp
Cp = Cp + (hr(k)-hl(k)) * Y(k)
enddo
Cp = Cp / (Tr-Tl)
else
T = 0.5d0*(Tr+Tl)
Cp = avgtabip( T, Y, cpk, Nsp )
endif
Cpmix = Cp
c
return
end