c
c
c
c ==================================================================
subroutine rpn3meuqs(ixyz,maxm,meqn,mwaves,mbc,mx,qlo,qro,
& maux,auxl,auxr,wave,s,amdq,apdq)
c ==================================================================
c
c # Two-component Roe-solver for the Euler equations with quasi-stationary
c # source term balancing.
c
c # solve Riemann problems along one slice of data.
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 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
c # auxl(i,ma,2) contains auxiliary data for cells along this slice,
c # where ma=1,maux in the case where maux=method(7) > 0.
c # auxl(i,ma,1) and auxl(i,ma,3) contain auxiliary data along
c # neighboring slices that generally aren't needed in the rpn3 routine.
c
c
c # On output, wave contains the waves, s the speeds,
c # and amdq, apdq the decomposition of the flux difference
c # f(qr(i-1)) - f(ql(i))
c # into leftgoing and rightgoing parts respectively.
c # With the Roe solver we have
c # amdq = A^- \Delta q and apdq = A^+ \Delta q
c # where A is the Roe matrix. An entropy fix can also be incorporated
c # into the flux differences.
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 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 qlo(1-mbc:maxm+mbc, meqn)
dimension qro(1-mbc:maxm+mbc, meqn)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(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 rpt3eu
c
c ------------
parameter (maxmrp = 1002) !# assumes atmost max(mx,my,mz) = 1000 with mbc=2
dimension delta(7), delh(7)
logical efix
common /comroe/ u2v2w2(-1:maxmrp),
& u(-1:maxmrp),v(-1:maxmrp),w(-1:maxmrp),enth(-1:maxmrp),
& a(-1:maxmrp),g1a2(-1:maxmrp),euv(-1:maxmrp),p(-1:maxmrp)
dimension ql(-1:maxmrp, 7), qr(-1:maxmrp, 7)
c
data efix /.true./ !# use entropy fix for transonic rarefactions
c
if (-1.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 = 2
mv = 3
mw = 4
ma = 1
else if(ixyz .eq. 2)then
mu = 3
mv = 4
mw = 2
ma = 2
else
mu = 4
mv = 2
mw = 3
ma = 3
endif
c
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 # Compute delh for left states
do 11 i=2-mbc,mx+mbc
gamma1 = 1.d0/qlo(i,6)
gamma = gamma1 + 1.d0
pinf = qlo(i,7) * gamma1/gamma
G = (gamma**2-gamma+2.)/gamma1
r = qlo(i,1)
ru = qlo(i,mu)
rv = qlo(i,mv)
rw = qlo(i,mw)
e = qlo(i,5)
delh(1) = auxl(i,ma,2)/(gamma1*(2.*(e-pinf)*
& gamma/r-ru**2*G/r**2))
delh(5) = auxl(i,ma,2)*r/(2.*e*gamma1)*(1.+ru**2*(3.-gamma)/
& (gamma1*(2.*(e-pinf)*gamma*r-G*ru**2)))
ql(i,1) = r*(1.-delh(1))
ql(i,mu) = ru
ql(i,mv) = rv
ql(i,mw) = rw
ql(i,5) = e*(1.-delh(5))
ql(i,6) = qlo(i,6)
ql(i,7) = qlo(i,7)
11 continue
c # Compute delh for right states
do 12 i=1-mbc,mx+mbc-1
gamma1 = 1.d0/qro(i,6)
gamma = gamma1 + 1.d0
pinf = qro(i,7) * gamma1/gamma
G = (gamma**2-gamma+2.)/gamma1
r = qro(i,1)
ru = qro(i,mu)
rv = qro(i,mv)
rw = qro(i,mw)
e = qro(i,5)
G = (gamma**2-gamma+2.)/gamma1
delh(1) = auxr(i,ma,2)/(gamma1*(2.*(e-pinf)*
* gamma/r-ru**2*G/r**2))
delh(5) = auxr(i,ma,2)*r/(2.*e*gamma1)*(1.+ru**2*(3.-gamma)/
& (gamma1*(2.*(e-pinf)*gamma*r-G*ru**2)))
qr(i,1) = r*(1.+delh(1))
qr(i,mu) = ru
qr(i,mv) = rv
qr(i,mw) = rw
qr(i,5) = e*(1.+delh(5))
qr(i,6) = qro(i,6)
qr(i,7) = qro(i,7)
12 continue
c do 13 i=2-mbc,mx+mbc-1
c pr = (qr(i,5) - 0.5d0*(qr(i,2)**2 + qr(i,3)**2 +
c & qr(i,4)**2)/qr(i,1) - qr(i,7) ) / qr(i,6)
c pl = (ql(i,5) - 0.5d0*(ql(i,2)**2 + ql(i,3)**2 +
c & ql(i,4)**2)/ql(i,1) - ql(i,7) ) / ql(i,6)
c cm = qr(i,mu)**2/qr(i,1)+pr-ql(i,mu)**2/ql(i,1)-pl
c ce = (qr(i,5)+pr)*qr(i,mu)/qr(i,1)-
c & (ql(i,5)+pl)*ql(i,mu)/ql(i,1)
c write(6,*) 'cm ',cm,qlo(i,1)*auxl(i,ma,2),cm-qlo(i,1)*auxl(i,ma,2)
c write(6,*) 'ce ',ce,qlo(i,mu)*auxl(i,ma,2),ce-qlo(i,mu)*auxl(i,ma,2)
c
c write(6,*) 'new ',ql(i,1),qlo(i,1),qr(i,1),ql(i,5),qlo(i,5),qr(i,5)
c 13 continue
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 rpt3eu to do the transverse wave
c # splitting.
c
do 10 i = 2-mbc, mx+mbc
rhsqrtl = dsqrt(qr(i-1,1))
rhsqrtr = dsqrt(ql(i,1))
pl = (qr(i-1,5) - 0.5d0*(qr(i-1,mu)**2 + qr(i-1,mv)**2 +
& qr(i-1,mw)**2)/qr(i-1,1) - qr(i-1,7) ) / qr(i-1,6)
pr = (ql(i,5) - 0.5d0*(ql(i,mu)**2 + ql(i,mv)**2 +
& ql(i,mw)**2)/ql(i,1) - ql(i,7) ) / ql(i,6)
rhsq2 = rhsqrtl + rhsqrtr
gamma1 = rhsq2 / ( qr(i-1,6)*rhsqrtl + ql(i,6)*rhsqrtr )
xjota = ( pl*qr(i-1,6)*rhsqrtl + pr*ql(i,6)*rhsqrtr ) / rhsq2
p(i) = xjota*gamma1
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
enth(i) = (((qr(i-1,5)+pl)/rhsqrtl
& + (ql(i,5)+pr)/rhsqrtr)) / rhsq2
u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2
a2 = gamma1*(enth(i) - .5d0*u2v2w2(i))
a(i) = dsqrt(a2)
g1a2(i) = gamma1 / a2
euv(i) = enth(i) - u2v2w2(i)
10 continue
c
c
c # now split the jump in q1d at each interface into waves
c
c # find a1 thru a5, the coefficients of the 5 eigenvectors:
do 20 i = 2-mbc, mx+mbc
delta(1) = ql(i,1) - qr(i-1,1)
delta(2) = ql(i,mu) - qr(i-1,mu)
delta(3) = ql(i,mv) - qr(i-1,mv)
delta(4) = ql(i,mw) - qr(i-1,mw)
delta(5) = ql(i,5) - qr(i-1,5)
delta(6) = ql(i,6) - qr(i-1,6)
delta(7) = ql(i,7) - qr(i-1,7)
a4 = g1a2(i) * (euv(i)*delta(1)
& + u(i)*delta(2) + v(i)*delta(3) + w(i)*delta(4)
& - delta(5) + p(i)*delta(6) + delta(7))
a2 = delta(3) - v(i)*delta(1)
a3 = delta(4) - w(i)*delta(1)
a5 = (delta(2) + (a(i)-u(i))*delta(1) - a(i)*a4) / (2.d0*a(i))
a1 = delta(1) - a4 - a5
a6 = delta(6)
a7 = delta(7)
c
c # Compute the waves.
c # Note that the 2,3,6,7-wave travel at the same speed
c # and are lumped together in wave(.,.,2). The 5-wave is then stored
c # in wave(.,.,3).
c
wave(i,1,1) = a1
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,5,1) = a1*(enth(i) - u(i)*a(i))
wave(i,6,1) = 0.d0
wave(i,7,1) = 0.d0
s(i,1) = u(i)-a(i)
c
wave(i,1,2) = a4
wave(i,mu,2) = a4*u(i)
wave(i,mv,2) = a4*v(i) + a2
wave(i,mw,2) = a4*w(i) + a3
wave(i,5,2) = a4*0.5d0*u2v2w2(i)+a2*v(i)+a3*w(i)+a6*p(i)+a7
wave(i,6,2) = a6
wave(i,7,2) = a7
s(i,2) = u(i)
c
wave(i,1,3) = a5
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,5,3) = a5*(enth(i)+u(i)*a(i))
wave(i,6,3) = 0.d0
wave(i,7,3) = 0.d0
s(i,3) = u(i)+a(i)
20 continue
c
c
c # compute flux differences amdq and apdq.
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 mws=1,mwaves
if (s(i,mws) .lt. 0.d0) then
amdq(i,m) = amdq(i,m) + s(i,mws)*wave(i,m,mws)
else
apdq(i,m) = apdq(i,m) + s(i,mws)*wave(i,m,mws)
endif
90 continue
100 continue
go to 900
c
c-----------------------------------------------------
c
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
rhoim1 = qr(i-1,1)
pim1 = (qr(i-1,5) - 0.5d0*(qr(i-1,mu)**2 + qr(i-1,mv)**2 +
& qr(i-1,mw)**2)/qr(i-1,1) - qr(i-1,7) ) / qr(i-1,6)
gamma1 = 1.d0/qr(i-1,6)
gamma = gamma1 + 1.d0
pinf = qr(i-1,7)*gamma1/gamma
cim1 = dsqrt(gamma*(pim1+pinf)/rhoim1)
s0 = qr(i-1,mu)/rhoim1 - cim1 !# u-c in left state (cell i-1)
c
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
rho1 = qr(i-1,1) + wave(i,1,1)
rhou1 = qr(i-1,mu) + wave(i,mu,1)
rhov1 = qr(i-1,mv) + wave(i,mv,1)
rhow1 = qr(i-1,mw) + wave(i,mw,1)
en1 = qr(i-1,5) + wave(i,5,1)
p1 = (en1 - 0.5d0*(rhou1**2 + rhov1**2 + rhow1**2)/rho1
& - qr(i-1,7) ) / qr(i-1,6)
c1 = dsqrt(gamma*(p1+pinf)/rho1)
s1 = rhou1/rho1 - c1 !# u-c to right of 1-wave
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-,3- and 4- waves are 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
rhoi = ql(i,1)
pi = (ql(i,5) - 0.5d0*(ql(i,mu)**2 + ql(i,mv)**2 +
& ql(i,mw)**2)/ql(i,1) - ql(i,7) ) / ql(i,6)
gamma1 = 1.d0/ql(i,6)
gamma = gamma1 + 1.d0
pinf = ql(i,7)*gamma1/gamma
ci = dsqrt(gamma*(pi+pinf)/rhoi)
s3 = ql(i,mu)/rhoi + ci !# u+c in right state (cell i)
c
rho2 = ql(i,1) - wave(i,1,3)
rhou2 = ql(i,mu) - wave(i,mu,3)
rhov2 = ql(i,mv) - wave(i,mv,3)
rhow2 = ql(i,mw) - wave(i,mw,3)
en2 = ql(i,5) - wave(i,5,3)
p2 = (en2 - 0.5d0*(rhou2**2+rhov2**2+rhow2**2)/rho2
& - ql(i,7)) / ql(i,6)
c2 = dsqrt(gamma*(p2+pinf)/rho2)
s2 = rhou2/rho2 + c2 !# u+c to left of 3-wave
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,5
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 mws=1,mwaves
df = df + s(i,mws)*wave(i,m,mws)
210 continue
apdq(i,m) = df - amdq(i,m)
220 continue
c
900 continue
return
end