c
c =========================================================
subroutine rpn3euznd(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,amdq,apdq)
c =========================================================
c
c # solve Riemann problems for the 3D ZND-Euler equations using Roe's
c # approximate Riemann solver.
c # Scheme is blended with HLL for robustness and uses a multi-dimensional
c # entropy correction to prevent the carbuncle phenomenon.
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,
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 routines, this routine is called with ql = qr
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h,o-z)
c
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 apdq(1-mbc:maxm+mbc, meqn)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension auxl(1-mbc:maxm+mbc, maux, 3)
dimension auxr(1-mbc:maxm+mbc, maux, 3)
c
c local arrays -- common block comroe is passed to rpt3euznd
c ------------
parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5
parameter (minmrp = -4) !# assumes at most mbc=5
dimension delta(6), fl(minmrp:maxmrp,6), fr(minmrp:maxmrp,6)
logical efix, pfix, hll, roe, hllfix
common /param/ gamma,gamma1,q0
common /comroe/ u2v2w2(minmrp:maxmrp),
& u(minmrp:maxmrp),v(minmrp:maxmrp),w(minmrp:maxmrp),
& enth(minmrp:maxmrp),a(minmrp:maxmrp),Y(2,minmrp:maxmrp)
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 flux differences
c ------------
common /rpnflx/ mrpnflx
mrpnflx = 0
c
if (minmrp.gt.1-mbc .or. maxmrp .lt. maxm+mbc) then
write(6,*) 'need to increase maxmrp 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 and mw to the
c # orthogonal momentums:
c
if(ixyz .eq. 1)then
mu = 3
mv = 4
mw = 5
else if(ixyz .eq. 2)then
mu = 4
mv = 5
mw = 3
else
mu = 5
mv = 3
mw = 4
endif
c
c # note that notation for u,v, and w reflects assumption that the
c # Riemann problems are in the x-direction with u in the normal
c # direction and v and w in the orthogonal directions, but with the
c # above definitions of mu, mv, and mw the routine also works with
c # ixyz=2 and ixyz = 3
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 rpt3euznd to do the transverse wave splitting.
c
do 10 i=2-mbc,mx+mbc
c
pl = gamma1*(qr(i-1,6) - qr(i-1,2)*q0 -
& 0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2+qr(i-1,mw)**2)/
& (qr(i-1,1)+qr(i-1,2)))
pr = gamma1*(ql(i, 6) - ql(i, 2)*q0 -
& 0.5d0*(ql(i, mu)**2+ql(i, mv)**2+ql(i, mw)**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,mu)/rhsqrtl + ql(i,mu)/rhsqrtr) / rhsq2
v(i) = (qr(i-1,mv)/rhsqrtl + ql(i,mv)/rhsqrtr) / rhsq2
w(i) = (qr(i-1,mw)/rhsqrtl + ql(i,mw)/rhsqrtr) / rhsq2
u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2
enth(i) = (((qr(i-1,6)+pl)/rhsqrtl
& + (ql(i ,6)+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*u2v2w2(i) - Y(2,i)*q0)
a(i) = dsqrt(a2)
c
10 continue
c
do 30 i=2-mbc,mx+mbc
c
c # find a1 thru a5, the coefficients of the 5 eigenvectors:
c
do k = 1, 6
delta(k) = ql(i,k) - qr(i-1,k)
enddo
drho = delta(1) + delta(2)
c
a2 = gamma1/a(i)**2 * (drho*0.5d0*u2v2w2(i) - delta(2)*q0
& - (u(i)*delta(mu)+v(i)*delta(mv)+w(i)*delta(mw))
& + delta(6))
a3 = delta(mv) - v(i)*drho
a4 = delta(mw) - w(i)*drho
a5 = 0.5d0*( a2 - ( u(i)*drho - delta(mu) )/a(i) )
a1 = a2 - a5
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,mu,1) = a1*(u(i) - a(i))
wave(i,mv,1) = a1*v(i)
wave(i,mw,1) = a1*w(i)
wave(i,6,1) = a1*(enth(i) - u(i)*a(i))
s(i,1) = u(i)-a(i)
c
c # 2-wave
wave(i,1,2) = delta(1) - Y(1,i)*a2
wave(i,2,2) = delta(2) - Y(2,i)*a2
wave(i,mu,2) = (drho - a2)*u(i)
wave(i,mv,2) = (drho - a2)*v(i) + a3
wave(i,mw,2) = (drho - a2)*w(i) + a4
wave(i,6,2) = (drho - a2)*0.5d0*u2v2w2(i) +
& q0*(delta(2) - Y(2,i)*a2) + a3*v(i) + a4*w(i)
s(i,2) = u(i)
c
c # 3-wave
wave(i,1,3) = a5*Y(1,i)
wave(i,2,3) = a5*Y(2,i)
wave(i,mu,3) = a5*(u(i) + a(i))
wave(i,mv,3) = a5*v(i)
wave(i,mw,3) = a5*w(i)
wave(i,6,3) = a5*(enth(i) + u(i)*a(i))
s(i,3) = u(i)+a(i)
c
30 continue
c
c # compute flux differences as
c # (+/-)
c # A (Ur-Ul) = 0.5*( f(Ur)-f(Ul) +/- |A|(Ur-Ul) )
c --------------------------
call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,auxr,apdq)
call flx3(ixyz,maxm,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 mws=1,mwaves
sl = dabs(s(i,mws))
if (efix.and.dabs(s(i,mws)).lt.auxl(i,ixyz,2))
& sl = s(i,mws)**2/(2.d0*auxl(i,ixyz,2))+
& 0.5d0*auxl(i,ixyz,2)
sw = sw + sl*wave(i,m,mws)
60 continue
amdq(i,m) = amdq(i,m) - 0.5d0*sw
apdq(i,m) = amdq(i,m) + sw
50 continue
endif
c
if (hll) then
do 55 i = 2-mbc, mx+mbc
hllfix = .false.
c
rho1l = qr(i-1,1) + wave(i,1,1)
rho2l = qr(i-1,2) + wave(i,2,1)
rhoul = qr(i-1,mu) + wave(i,mu,1)
rhovl = qr(i-1,mv) + wave(i,mv,1)
rhowl = qr(i-1,mw) + wave(i,mw,1)
El = qr(i-1,6) + wave(i,6,1)
pl = gamma1*(El - rho2l*q0 -
& 0.5d0*(rhoul**2 + rhovl**2 + rhowl**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,mu) - wave(i,mu,3)
rhovr = ql(i,mv) - wave(i,mv,3)
rhowr = ql(i,mw) - wave(i,mw,3)
Er = ql(i,6 ) - wave(i,6,3)
pr = gamma1*(Er - rho2r*q0 -
& 0.5d0*(rhour**2 + rhovr**2 + rhowr**2)/(rho1r+rho2r))
if (rho1r+rho2r.le.0.d0.or.pr.le.0.d0)
& hllfix = .true.
c
if (hllfix) then
c write (6,*) 'Switching to HLL in',i
c
rl = qr(i-1,1) + qr(i-1,2)
ul = qr(i-1,mu)/rl
pl = gamma1*(qr(i-1,6) - qr(i-1,2)*q0 - 0.5d0*
& (qr(i-1,mu)**2+qr(i-1,mv)**2+qr(i-1,mw)**2)/rl)
al = dsqrt(gamma*pl/rl)
c
rr = ql(i ,1) + ql(i ,2)
ur = ql(i ,mu)/rr
pr = gamma1*(ql(i ,6) - ql(i ,2)*q0 - 0.5d0*
& (ql(i ,mu)**2+ql(i ,mv)**2+ql(i ,mw)**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)
apdq(i,m) = -(fg-fl(i ,m))
enddo
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 = amdq(i,1)+amdq(i,2)
apdr = apdq(i,1)+apdq(i,2)
rhol = qr(i-1,1)+qr(i-1,2)
rhor = ql(i ,1)+ql(i ,2)
do 70 m=1,2
if (qr(i-1,mu)+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,mu)
apdq(i,m) = Z*apdr - (Z-ql(i ,m)/rhor)*ql(i ,mu)
70 continue
endif
c
return
end
c