c
c =========================================================
subroutine rpn2eurhok(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the thermally perfect 2D multi-component
c # Euler equations using Roe's 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 # This data is along a slice in the x-direction if ixy=1
c # or the y-direction if ixy=2.
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, Georg Bader
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h,o-z)
c
include "ck.i"
c
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension s(1-mbc:maxm+mbc, mwaves)
dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(1-mbc:maxm+mbc, meqn)
c
c # local storage
c ---------------
parameter (maxm2 = 10005) !# assumes at most 10000x10000 grid with mbc=5
parameter (minm2 = -4) !# assumes at most mbc=5
common /comroe/ u(minm2:maxm2), v(minm2:maxm2), u2v2(minm2:maxm2),
& enth(minm2:maxm2), a(minm2:maxm2), g1a2(minm2:maxm2),
& dpY(minm2:maxm2), Y(LeNsp,minm2:maxm2), pk(LeNsp,minm2:maxm2)
logical efix, pfix
c
c define local arrays
c
dimension delta(LeNsp+3)
dimension rkl(LeNsp), rkr(LeNsp)
dimension hkl(LeNsp), hkr(LeNsp)
c
data efix /.true./ !# use entropy fix for transonic rarefactions
data pfix /.true./ !# use Larrouturou's positivity fix for species
c
c # Riemann solver returns flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 0
c
if (minm2.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then
write(6,*) 'need to increase maxm2 in rpA'
stop
endif
c
c # set mu to point to the component of the system that corresponds
c # to momentum in the direction of this slice, mv to the orthogonal
c # momentum:
c
if (ixy.eq.1) then
mu = Nsp+1
mv = Nsp+2
else
mu = Nsp+2
mv = Nsp+1
endif
mE = Nsp+3
mT = Nsp+4
c
c # note that notation for u and v reflects assumption that the
c # Riemann problems are in the x-direction with u in the normal
c # direciton and v in the orthogonal direcion, but with the above
c # definitions of mu and mv the routine also works with ixy=2
c # and returns, for example, f0 as the Godunov flux g0 for the
c # Riemann problems u_t + g(u)_y = 0 in the y-direction.
c
c
c # compute the Roe-averaged variables needed in the Roe solver.
c # These are stored in the common block comroe since they are
c # later used in routine rpt2eurhok to do the transverse wave splitting.
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 # left/right Temperatures already calculated
c
rhoel = qr(i-1,mE)-0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2)/rhol
call SolveTrhok(qr(i-1,mT),rhoel,rhoWl,rkl,Nsp,ifail)
rhoer = ql(i ,mE)-0.5d0*(ql(i ,mu)**2+ql(i ,mv)**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
v(i) = (qr(i-1,mv)/rhsqrtl + ql(i,mv)/rhsqrtr) / rhsq2
u2v2(i) = u(i)**2 + v(i)**2
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
c
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
pk(k,i) = 0.5d0*u2v2(i) - hk + tmp / Wk(k)
enddo
c
c # compute speed of sound
c
dpY(i) = 0.d0
do k = 1, Nsp
dpY(i) = dpY(i) + pk(k,i) * Y(k,i)
enddo
a2 = dpY(i) + enth(i)-u2v2(i)
g1a2(i) = 1.d0 / a2
a(i) = dsqrt(gamma1*a2)
c
20 continue
c
c
do 30 i=2-mbc,mx+mbc
c
c # find a1 thru a4, the coefficients of the Nsp+3 eigenvectors:
c
dpdr = 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)
enddo
delta(mu) = ql(i,mu) - qr(i-1,mu)
delta(mv) = ql(i,mv) - qr(i-1,mv)
delta(mE) = ql(i,mE) - qr(i-1,mE)
c
a2 = g1a2(i)*(dpdr - ( u(i)*delta(mu) + v(i)*delta(mv) )
& + delta(mE) )
a3 = delta(mv) - v(i)*drho
a4 = 0.5d0*( a2 - ( u(i)*drho - delta(mu) )/a(i) )
a1 = a2 - a4
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) = a4*Y(k,i)
enddo
c # 1-wave
wave(i,mu,1) = a1*(u(i) - a(i))
wave(i,mv,1) = a1*v(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) = (drho - a2)*u(i)
wave(i,mv,2) = (drho - a2)*v(i) + a3
wave(i,mE,2) = (drho - a2)*u2v2(i)
& - dpdr + dpY(i)*a2 + a3*v(i)
wave(i,mT,2) = 0.d0
s(i,2) = u(i)
c
c # 3-wave
wave(i,mu,3) = a4*(u(i) + a(i))
wave(i,mv,3) = a4*v(i)
wave(i,mE,3) = a4*(enth(i) + u(i)*a(i))
wave(i,mT,3) = 0.d0
s(i,3) = u(i)+a(i)
c
30 continue
c
c # compute Godunov flux f0:
c --------------------------
c
c
if (efix) go to 110
c
c # no entropy fix
c ----------------
c
c # amdq = SUM s*wave over left-going waves
c # apdq = SUM s*wave over right-going waves
c
do 100 m=1,meqn
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
110 continue
c
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
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rkl(k) = qr(i-1,k)
rho = rho + rkl(k)
rhoW = rhoW + rkl(k)/Wk(k)
enddo
rhou = qr(i-1,mu)
rhov = qr(i-1,mv)
rhoE = qr(i-1,mE)
T = qr(i-1,mT)
rhoCp = avgtabip( T, rkl, cpk, Nsp )
gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0
p = rhoW*RU*T
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
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rkl(k) = rkl(k) + wave(i,k,1)
rho = rho + rkl(k)
rhoW = rhoW + rkl(k)/Wk(k)
enddo
rhou = rhou + wave(i,mu,1)
rhov = rhov + wave(i,mv,1)
rhoE = rhoE + wave(i,mE,1)
rhoe = rhoE - 0.5d0*(rhou**2+rhov**2)/rho
if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then
write(6,*) 'Temperature out of range', T
write(6,*) 'state vector 1 before'
write(6,*) (rkl(k),k=1,Nsp)
endif
call SolveTrhok( T, rhoe, rhoW, rkl, Nsp, ifail)
rhoCp = avgtabip( T, rkl, cpk, Nsp )
if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then
write(6,*) 'Temperature out of range', T
write(6,*) 'state vector 1 after'
write(6,*) (rkl(k),k=1,Nsp)
endif
gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0
p = rhoW*RU*T
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
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rkr(k) = ql(i,k)
rho = rho + rkr(k)
rhoW = rhoW + rkr(k)/Wk(k)
enddo
rhou = ql(i,mu)
rhov = ql(i,mv)
rhoE = ql(i,mE)
T = ql(i,mT)
rhoCp = avgtabip( T, rkr, cpk, Nsp )
gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0
p = rhoW*RU*T
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
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rkr(k) = rkr(k) - wave(i,k,3)
rho = rho + rkr(k)
rhoW = rhoW + rkr(k)/Wk(k)
enddo
rhou = rhou - wave(i,mu,3)
rhov = rhov - wave(i,mv,3)
rhoE = rhoE - wave(i,mE,3)
rhoe = rhoE - 0.5d0*(rhou**2+rhov**2)/rho
if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then
write(6,*) 'Temperature out of range', T
write(6,*) 'state vector 1 before'
write(6,*) (rkr(k),k=1,Nsp)
endif
call SolveTrhok( T, rhoe, rhoW, rkr, Nsp, ifail)
rhoCp = avgtabip( T, rkr, cpk, Nsp )
if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then
write(6,*) 'Temperature out of range', T
write(6,*) 'state vector 1 after'
write(6,*) (rkr(k),k=1,Nsp)
endif
gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0
p = rhoW*RU*T
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
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,meqn
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
c
900 continue
c
if (pfix) then
do 70 i=2-mbc,mx+mbc
amdr = 0.d0
apdr = 0.d0
rhol = 0.d0
rhor = 0.d0
do k = 1, Nsp
amdr = amdr + amdq(i,k)
apdr = apdr + apdq(i,k)
rhol = rhol + qr(i-1,k)
rhor = rhor + ql(i ,k)
enddo
do 70 k=1,Nsp
if (qr(i-1,mu)+amdr.gt.0.d0) then
Z = qr(i-1,k)/rhol
else
Z = ql(i ,k)/rhor
endif
amdq(i,k) = Z*amdr + (Z-qr(i-1,k)/rhol)*qr(i-1,mu)
apdq(i,k) = Z*apdr - (Z-ql(i ,k)/rhor)*ql(i ,mu)
70 continue
endif
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