c
c =====================================================
subroutine setaux(maxmx,maxmy,maxmz,meqn,mbc,ibx,iby,ibz,
& mx,my,mz,q,aux,maux,xc,yc,zc,dx,dy,dz,t,dt)
c =====================================================
c
c # Multi-dimensional computation of wave speeds, used in
c # multi-dimensional entropy correction in rpn2euefix.
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h, o-z)
c
include "ck.i"
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc,1-ibz*mbc:maxmz+ibz*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc,1-ibz*mbc:maxmz+ibz*mbc)
dimension rk(LeNsp), rkl(LeNsp), rkd(LeNsp), rkb(LeNsp)
c
mu = Nsp+1
mv = Nsp+2
mw = Nsp+3
mE = Nsp+4
mT = Nsp+5
c
do 5 k = 1-ibz*mbc, mz+ibz*mbc
do 5 j = 1-iby*mbc, my+iby*mbc
do 5 i = 1-ibx*mbc, mx+ibx*mbc
rho = 0.d0
rhoW = 0.d0
do m = 1, Nsp
rk(m) = q(m,i,j,k)
rho = rho + q(m,i,j,k)
rhoW = rhoW + q(m,i,j,k)/Wk(m)
enddo
rhoe = q(mE,i,j,k)-0.5d0*(q(mu,i,j,k)**2+q(mv,i,j,k)**2+
& q(mw,i,j,k)**2)/rho
call SolveTrhok(q(mT,i,j,k),rhoe,rhoW,rk,Nsp,ifail)
5 continue
c
do 10 k = 1-ibz*mbc, mz+ibz*mbc
do 10 j = 1-iby*mbc, my+iby*mbc
do 10 i = 1-ibx*mbc, mx+ibx*mbc
rho = 0.d0
rhoW = 0.d0
do m = 1, Nsp
rk(m) = q(m,i,j,k)
rho = rho + q(m,i,j,k)
rhoW = rhoW + q(m,i,j,k)/Wk(m)
enddo
u = q(mu,i,j,k)/rho
v = q(mv,i,j,k)/rho
w = q(mw,i,j,k)/rho
T0 = q(mT,i,j,k)
p = rhoW*RU*T0
rhoCp = avgtabip( T0, rk, cpk, Nsp )
gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0
a = dsqrt(gamma*p/rho)
c
if (i.gt.1-ibx*mbc) then
rhol = 0.d0
rhoWl = 0.d0
do m = 1, Nsp
rkl(m) = q(m,i-1,j,k)
rhol = rhol + q(m,i-1,j,k)
rhoWl = rhoWl + q(m,i-1,j,k)/Wk(m)
enddo
ul = q(mu,i-1,j,k)/rhol
Tl = q(mT,i-1,j,k)
pl = rhoWl*RU*Tl
rhoCpl = avgtabip( Tl, rkl, cpk, Nsp )
gammal = RU / ( rhoCpl/rhoWl - RU ) + 1.d0
al = dsqrt(gammal*pl/rhol)
aux(4,i,j,k) = dmax1(dabs(u-a-(ul-al)),dabs(u-ul),
& dabs(u+a-(ul+al)))
else
aux(4,i,j,k) = 0.d0
endif
c
if (j.gt.1-iby*mbc) then
rhod = 0.d0
rhoWd = 0.d0
do m = 1, Nsp
rkd(m) = q(m,i,j-1,k)
rhod = rhod + q(m,i,j-1,k)
rhoWd = rhoWd + q(m,i,j-1,k)/Wk(m)
enddo
vd = q(mv,i,j-1,k)/rhod
Td = q(mT,i,j-1,k)
pd = rhoWd*RU*Td
rhoCpd = avgtabip( Td, rkd, cpk, Nsp )
gammad = RU / ( rhoCpd/rhoWd - RU ) + 1.d0
ad = dsqrt(gammad*pd/rhod)
aux(5,i,j,k) = dmax1(dabs(v-a-(vd-ad)),dabs(v-vd),
& dabs(v+a-(vd+ad)))
else
aux(5,i,j,k) = 0.d0
endif
c
if (k.gt.1-ibz*mbc) then
rhob = 0.d0
rhoWb = 0.d0
do m = 1, Nsp
rkb(m) = q(m,i,j,k-1)
rhob = rhob + q(m,i,j,k-1)
rhoWb = rhoWb + q(m,i,j,k-1)/Wk(m)
enddo
wb = q(mw,i,j,k-1)/rhob
Tb = q(mT,i,j,k-1)
pb = rhoWb*RU*Tb
rhoCpb = avgtabip( Tb, rkb, cpk, Nsp )
gammab = RU / ( rhoCpb/rhoWb - RU ) + 1.d0
ab = dsqrt(gammab*pb/rhob)
aux(6,i,j,k) = dmax1(dabs(w-a-(wb-ab)),dabs(w-wb),
& dabs(w+a-(wb+ab)))
else
aux(6,i,j,k) = 0.d0
endif
10 continue
c
do 11 j = 1-iby*mbc, my+iby*mbc
do 11 k = 1-ibz*mbc, mz+ibz*mbc
aux(1,1-ibx*mbc,j,k) = 0.d0
aux(2,1-ibx*mbc,j,k) = 0.d0
aux(3,1-ibx*mbc,j,k) = 0.d0
11 continue
c
do 12 i = 1-ibx*mbc, mx+ibx*mbc
do 12 k = 1-ibz*mbc, mz+ibz*mbc
aux(1,i,1-iby*mbc,k) = 0.d0
aux(2,i,1-iby*mbc,k) = 0.d0
aux(3,i,1-iby*mbc,k) = 0.d0
12 continue
c
do 13 i = 1-ibx*mbc, mx+ibx*mbc
do 13 j = 1-iby*mbc, my+iby*mbc
aux(1,i,j,1-ibz*mbc) = 0.d0
aux(2,i,j,1-ibz*mbc) = 0.d0
aux(3,i,j,1-ibz*mbc) = 0.d0
13 continue
c
do 20 k = 2-ibz*mbc, mz+ibz*mbc
do 20 j = 2-iby*mbc, my+iby*mbc
do 20 i = 2-ibx*mbc, mx+ibx*mbc
aux(1,i,j,k) = dmax1(aux(4,i,j,k),
& aux(5,i,j,k),aux(5,i-1,j,k),
& aux(6,i,j,k),aux(6,i-1,j,k))
if (j.lt.my+iby*mbc)
& aux(1,i,j,k) = dmax1(aux(1,i,j,k),aux(5,i ,j+1,k),
& aux(5,i-1,j+1,k))
if (k.lt.mz+ibz*mbc)
& aux(1,i,j,k) = dmax1(aux(1,i,j,k),aux(6,i ,j,k+1),
& aux(6,i-1,j,k+1))
c
aux(2,i,j,k) = dmax1(aux(5,i,j,k),
& aux(4,i,j,k),aux(4,i,j-1,k),
& aux(6,i,j,k),aux(6,i,j-1,k))
if (i.lt.mx+ibx*mbc)
& aux(2,i,j,k) = dmax1(aux(2,i,j,k),aux(4,i+1,j ,k),
& aux(4,i+1,j-1,k))
if (k.lt.mz+ibz*mbc)
& aux(2,i,j,k) = dmax1(aux(2,i,j,k),aux(6,i,j ,k+1),
& aux(6,i,j-1,k+1))
c
aux(3,i,j,k) = dmax1(aux(6,i,j,k),
& aux(4,i,j,k),aux(4,i,j,k-1),
& aux(5,i,j,k),aux(5,i,j,k-1))
if (i.lt.mx+ibx*mbc)
& aux(3,i,j,k) = dmax1(aux(3,i,j,k),aux(4,i+1,j,k ),
& aux(4,i+1,j,k-1))
if (j.lt.my+iby*mbc)
& aux(3,i,j,k) = dmax1(aux(3,i,j,k),aux(5,i,j+1,k ),
& aux(5,i,j+1,k-1))
20 continue
c
return
end
c