c
c
c ==================================================================
subroutine rpn3eu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
& maux,auxl,auxr,wave,s,fl,fr)
c ==================================================================
c
c # solve Riemann problems for the 3D Euler equations using
c # an improved version of the Liou-Steffen Flux-Vector-Splitting
c
c # Yasuhiro Wada, Meng-Sing Liou "An accurate and robust flux
c # splitting scheme for shock and contact discontinuities",
c # SIAM J. Sci. Comput., Vol. 18, No.2, pp 633-657, May 1997.
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 ixyz=1
c # the y-direction if ixyz=2.
c # the z-direction if ixyz=3.
c
c # On output, wave contains the waves,
c # s the speeds, fl and fr the positive and negative 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 routines, 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 wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension s(1-mbc:maxm+mbc, mwaves)
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension fl(1-mbc:maxm+mbc, meqn)
dimension fr(1-mbc:maxm+mbc, meqn)
dimension auxl(1-mbc:maxm+mbc, maux, 3)
dimension auxr(1-mbc:maxm+mbc, maux, 3)
double precision l(5), r(5)
common /param/ gamma,gamma1
c
c # Method returns fluxes
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 1
c
c # set mu to point to the component of the system that corresponds
c # to momentum in the direction of this slice, mv and mw to the
c # orthogonal momentum:
c
if(ixyz .eq. 1)then
mu = 2
mv = 3
mw = 4
else if(ixyz .eq. 2)then
mu = 3
mv = 4
mw = 2
else
mu = 4
mv = 2
mw = 3
endif
c
c # AUSM Flux Vector Splitting
c
do 10 i=2-mbc,mx+mbc
rhol = qr(i-1,1)
rhor = ql(i ,1)
ul = qr(i-1,mu)/rhol
ur = ql(i ,mu)/rhor
vl = qr(i-1,mv)/rhol
vr = ql(i ,mv)/rhor
wl = qr(i-1,mw)/rhol
wr = ql(i ,mw)/rhor
pl = gamma1*(qr(i-1,5) - 0.5d0*(ul**2+vl**2+wl**2)*rhol)
pr = gamma1*(ql(i ,5) - 0.5d0*(ur**2+vr**2+wr**2)*rhor)
Hl = (qr(i-1,5)+pl)/rhol
Hr = (ql(i ,5)+pr)/rhor
al = dsqrt(gamma*pl/rhol)
ar = dsqrt(gamma*pr/rhor)
c
am = dmax1(al,ar)
alphal = 2.d0*(pl/rhol)/(pl/rhol+pr/rhor)
alphar = 2.d0*(pr/rhor)/(pl/rhol+pr/rhor)
c
ulp = 0.5d0*(ul+dabs(ul))
plp = pl*ulp/ul
if (dabs(ul).le.am) then
ulp = 0.25d0*alphal*(ul+al)**2/am + (1.d0-alphal)*ulp
plp = 0.25d0*pl*(ul+al)**2/am**2*(2.d0-ul/am)
endif
c
urm = 0.5d0*(ur-dabs(ur))
prm = pr*urm/ur
if (dabs(ur).le.am) then
urm = -0.25d0*alphar*(ur-ar)**2/am + (1.d0-alphar)*urm
prm = 0.25d0*pr*(ur-ar)**2/am**2*(2.d0+ur/am)
endif
c
c # Blending between AUSMV and AUSMD
c # sf=1.d0 gives AUSMV, sf=-1.d0 gives AUSMD
sf = dmin1(1.d0, 10.d0*dabs(pr-pl)/dmin1(pl,pr))
c
l(1) = 0.5d0*(ulp*rhol+dabs(ulp*rhol))
l(mu) = 0.5d0*((1.d0+sf)*ulp*rhol*ul + (1.d0-sf)*l(1)*ul) + plp
l(mv) = l(1)*vl
l(mw) = l(1)*wl
l(5) = l(1)*Hl
c
r(1) = 0.5d0*(urm*rhor-dabs(urm*rhor))
r(mu) = 0.5d0*((1.d0+sf)*urm*rhor*ur + (1.d0-sf)*r(1)*ur) + prm
r(mv) = r(1)*vr
r(mw) = r(1)*wr
r(5) = r(1)*Hr
c
do 20 m = 1,meqn
fl(i,m) = l(m) + r(m)
fr(i,m) = -fl(i,m)
20 continue
c
s(i,1) = 0.5d0*(ul-al + ur-ar)
s(i,2) = 0.5d0*(ul + ur)
s(i,3) = 0.5d0*(ul+al + ur+ar)
do 10 mws=1,mwaves
do 10 m=1,meqn
wave(i,m,mws) = 0.d0
10 continue
c
return
end
c
c