c
c =========================================================
subroutine rpn3euznd(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
& auxl,auxr,wave,s,dfl,dfr)
c =========================================================
c
c # Riemann solver for the 3D ZND-Euler equations.
c # The waves are computed using the Roe approximation.
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 # dfl the left-going flux difference A^- \Delta q
c # dfr 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 # Author: Ralf Deiterding (based on rpn3euexact.f)
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 dfl(1-mbc:maxm+mbc, meqn)
dimension dfr(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)
dimension f0(minmrp:maxmrp,6), fl(minmrp:maxmrp,6), fr(minmrp:maxmrp,6)
dimension sl(2),sr(2)
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
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 Godunov flux f0:
c --------------------------
c
c # compute Godunov flux f0 at each interface.
c # Uses exact Riemann solver
c
do 200 i = 2-mbc, mx+mbc
c
rhol = qr(i-1,1) + qr(i-1,2)
rhor = ql(i ,1) + qr(i ,2)
Y2l = qr(i-1,2)/rhol
Y2r = ql(i ,2)/rhor
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,6) - qr(i-1,2)*q0 -
& 0.5d0*(ul**2+vl**2+wl**2)*rhol)
pr = gamma1*(ql(i, 6) - ql(i, 2)*q0 -
& 0.5d0*(ur**2+vr**2+wr**2)*rhor)
c
c # iterate to find pstar, ustar:
c
alpha = 1.
pstar = 0.5*(pl+pr)
wsr = dsqrt(pr*rhor) * phi(pstar/pr)
wsl = dsqrt(pl*rhol) * phi(pstar/pl)
c if (pl.eq.pr .and. rhol.eq.rhor) go to 60
c
40 do 50 iter=1,1000
p1 = (ul-ur+pr/wsr+pl/wsl) / (1./wsr + 1./wsl)
pstar = dmax1(p1,1d-6)*alpha + (1.-alpha)*pstar
wr1 = wsr
wl1 = wsl
wsr = dsqrt(pr*rhor) * phi(pstar/pr)
wsl = dsqrt(pl*rhol) * phi(pstar/pl)
if (dmax1(abs(wr1-wsr),dabs(wl1-wsl)) .lt. 1d-6)
& go to 60
50 continue
c
c # nonconvergence:
alpha = alpha/2.
if (alpha .gt. 0.1) go to 40
write(6,*) 'no convergence',i,wr1,wsr,wl1,wsl
wsr = .5*(wsr+wr1)
wsl = .5*(wsl+wl1)
c
60 continue
ustar = (pl-pr+wsr*ur+wsl*ul) / (wsr+wsl)
c
c # left wave:
c ============
c
if (pstar .gt. pl) then
c
c # shock:
sl(1) = ul - wsl/rhol
sr(1) = sl(1)
rho1 = wsl/(ustar-sl(1))
c
else
c
c # rarefaction:
cl = dsqrt(gamma*pl/rhol)
cstar = cl + 0.5*gamma1*(ul-ustar)
sl(1) = ul-cl
sr(1) = ustar-cstar
rho1 = (pstar/pl)**(1./gamma) * rhol
endif
c
c
c
c # right wave:
c =============
c
if (pstar .ge. pr) then
c
c # shock
sl(2) = ur + wsr/rhor
sr(2) = sl(2)
rho2 = wsr/(sl(2)-ustar)
c
else
c
c # rarefaction:
cr = dsqrt(gamma*pr/rhor)
cstar = cr + 0.5*gamma1*(ustar-ur)
sr(2) = ur+cr
sl(2) = ustar+cstar
rho2 = (pstar/pr)**(1./gamma)*rhor
endif
c
c
c # compute flux:
c ===============
c
c # compute state (rhos,us,ps) at x/t = 0:
c
if (sl(1).gt.0) then
rhos = rhol
us = ul
ps = pl
vs = vl
ws = wl
Y2s = Y2l
else if (sr(1).le.0. .and. ustar.ge. 0.) then
rhos = rho1
us = ustar
ps = pstar
vs = vl
ws = wl
Y2s = Y2l
else if (ustar.lt.0. .and. sl(2).ge. 0.) then
rhos = rho2
us = ustar
ps = pstar
vs = vr
ws = wr
Y2s = Y2r
else if (sr(2).lt.0) then
rhos = rhor
us = ur
ps = pr
vs = vr
ws = wr
Y2s = Y2r
else if (sl(1).le.0. .and. sr(1).ge.0.) then
c # transonic 1-rarefaction
us = (gamma1*ul + 2.*cl)/(gamma+1.)
e0 = pl/(rhol**gamma)
rhos = (us**2/(gamma*e0))**(1./gamma1)
ps = e0*rhos**gamma
vs = vl
ws = wl
Y2s = Y2l
else if (sl(2).le.0. .and. sr(2).ge.0.) then
c # transonic 3-rarefaction
us = (gamma1*ur - 2.*cr)/(gamma+1.)
e0 = pr/(rhor**gamma)
rhos = (us**2/(gamma*e0))**(1./gamma1)
ps = e0*rhos**gamma
vs = vr
ws = wr
Y2s = Y2r
endif
c
f0(i,1) = (1.d0-Y2s)*rhos*us
f0(i,2) = Y2s*rhos*us
f0(i,mu) = rhos*us**2 + ps
f0(i,mv) = rhos*us*vs
f0(i,mw) = rhos*us*ws
f0(i,6) = us*(gamma*ps/gamma1 + Y2s*rhos*q0 +
& 0.5*rhos*(us**2+vs**2+ws**2))
200 continue
c
c # compute fluxes in each cell:
c
call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,auxr,dfr)
call flx3(ixyz,maxm,meqn,mbc,mx,ql,maux,auxl,dfl)
c
do 210 m=1,6
do 210 i = 1-mbc, mx+mbc
fr(i,m) = dfr(i,m)
fl(i,m) = dfl(i,m)
210 continue
c
c # compute the leftgoing and rightgoing flux differences:
do 220 m=1,6
do 220 i = 2-mbc, mx+mbc
dfl(i,m) = f0(i,m) - fr(i-1,m)
dfr(i,m) = fl(i,m) - f0(i,m)
220 continue
c
return
end
c
c
double precision function phi(w)
implicit double precision (a-h,o-z)
common /param/ gamma,gamma1,q0
c
sqg = dsqrt(gamma)
if (w .gt. 1.) then
phi = dsqrt(w*(gamma+1.)/2. + gamma1/2.)
else if (w .gt. 0.99999) then
phi = sqg
else if (w .gt. .999) then
phi = sqg + (2*gamma**2 - 3.*gamma + 1)
& *(w-1.) / (4.*sqg)
else
phi = gamma1*(1.-w) / (2.*sqg*(1.-w**(gamma1/(2.*gamma))))
endif
return
end