c $Header: /home/proj/vtfng/vtf/amroc/rim/src/generic_rim2d.f,v 1.4 2009/02/05 02:52:29 cummings Exp $
subroutine funcaux(nx,ny,nbc)
implicit double precision (a-h,o-z)
include 'rim2d.i'
c
imin=1-nbc
imax=nx+nbc
jmin=1-nbc
jmax=ny+nbc
do j=jmin,jmax
do i=imin,imax
TT = p(i,j)/den(i,j)
ss(i,j) = dsqrt(gamma(i,j)*TT)
st5(i,j) =-rcon*z(i,j)*dexp(-eact/TT)
enddo
enddo
c
return
end
c
c------------------------------------------------------------------
c
subroutine traceback(i,j,ifc,pf,denf,uf,vf,zf,gammaf,dx,dy,dt)
implicit double precision (a-h,o-z)
include 'rim2d.i'
common/Bldir/dn1,dn2,tandn1,tandn2
common/Bltrackinfo/il,jl,ir,jr,dxl,dxr,
& dyl,dyr,rnx,rny,dxup,dyup,dxcfl
c
if(ifc.eq.1) then
il = i-1
jl = j
ir = i
jr = j
dxl = 0.5*dx ! (x(i,k)+x(i,k+1)) - xc(il,jl)
dyl = 0.d0 !0.5*(y(i,k)+y(i,k+1)) - yc(il,jl)
dxr =-0.5*dx ! (x(i,k)+x(i,k+1)) - xc(ir,jr)
dyr = 0.d0 !0.5*(y(i,k)+y(i,k+1)) - yc(ir,jr)
rnx = 1.d0 !rnx1(i,k)
rny = 0.d0 !rny1(i,k)
dxcfl= abs(dxl)+abs(dxr)
else
il = i
jl = j-1
ir = i
jr = j
dxl = 0.d0 !0.5*(x(i,k)+x(i+1,k)) - xc(il,jl)
dyl = 0.5*dy !(y(i,k)+y(i+1,k)) - yc(il,jl)
dxr = 0.d0 !0.5*(x(i,k)+x(i+1,k)) - xc(ir,jr)
dyr =-0.5*dy !(y(i,k)+y(i+1,k)) - yc(ir,jr)
rnx = 0.d0 !rnx2(i,k)
rny = 1.d0 !rny2(i,k)
dxcfl= abs(dyl)+abs(dyr)
endif
c
dxup = 0.75*0.5d0*dx !abs(0.5*(x(i,k)+x(i,k+1)) - xc(il,jl))
dyup = 0.75*0.5d0*dy !abs(0.5*(y(i,k)+y(i+1,k)) - yc(il,jl))
c write(*,10) '# ', dxup, dyup
c 10 format(a,f9.6,1x,f9.6)
c
ifind = 1
lside = 1
call rim(i,j,ifc,lside,pl,unl,denl,ifind,dt)
lside = -1
call rim(i,j,ifc,lside,pr,unr,denr,ifind,dt)
if(ifind.eq.0) call muscl(i,j,ifc,pl,denl,unl,pr,denr,unr,dt)
c
c - - - - - use Riemann solver for values at interfaces - - - - - -
c
approx = 1.0D-3
sm = 1.05d0
pratio = 1d0+2d0*gamma(i,j)/(gamma(i,j)+1d0)*(sm*sm-1d0)
Gammal = gamma(il,jl)
Gammar = gamma(ir,jr)
c
failriem = 0.0
call riemann(Gammar,Gammal,pr,pl,denr,denl,unr,unl,
& pf,unf,denfr,denfl,rmachr,rmachl,failriem,approx)
if(((pf/pr).gt.pratio).or.((pf/pl).gt.pratio)) then
approx=1.0D-12
failriem = 0.0
call riemann(Gammar,Gammal,pr,pl,denr,denl,unr,unl,
& pf,unf,denfr,denfl,rmachr,rmachl,failriem,approx)
endif
c
c call isentrope(i,j,pf,denf,ifind,dt)
c if(ifind.eq.0)
call denface(Gammar,Gammal,pr,pl,denr,
& denl,unr,unl,rmachr,rmachl,pf,unf,denfr,denfl,denf)
c
call tancom(i,j,utf,dt)
call isorate(i,j,zf,dt)
c
uf=unf*dn1+utf*tandn1
vf=unf*dn2+utf*tandn2
c
gammaf = Gammal
if(unf.lt.0d0) gammaf = Gammar
c
return
end
c
c-----------------------------------------------------------------
c
subroutine rim(j,k,ifc,lside,pup,unup,denup,ifind,dt)
implicit double precision (a-h,o-z)
include 'rim2d.i'
common/Bldir/dn1,dn2,tandn1,tandn2
common/Bltrackinfo/il,jl,ir,jr,dxl,dxr,
& dyl,dyr,rnx,rny,dxup,dyup,dxcfl
if(ifind.eq.0) go to 444
eps = 1d-18
c
If(lside.eq.1) Then
rimdn1=dn1
rimdn2=dn2
Else
rimdn1=-dn1
rimdn2=-dn2
Endif
c
divunx = dux(il,jl)*rimdn1+dvx(il,jl)*rimdn2
divuny = duy(il,jl)*rimdn1+dvy(il,jl)*rimdn2
gradRnx = dpx(il,jl)+den(il,jl)*ss(il,jl)*divunx
gradRny = dpy(il,jl)+den(il,jl)*ss(il,jl)*divuny
gradRn = dsqrt(gradRnx*gradRnx+gradRny*gradRny) + eps
c
divunxr = dux(ir,jr)*rimdn1+dvx(ir,jr)*rimdn2
divunyr = duy(ir,jr)*rimdn1+dvy(ir,jr)*rimdn2
gradRnxr = dpx(ir,jr)+den(ir,jr)*ss(ir,jr)*divunxr
gradRnyr = dpy(ir,jr)+den(ir,jr)*ss(ir,jr)*divunyr
gradRnr = dsqrt(gradRnxr*gradRnxr+gradRnyr*gradRnyr) + eps
c
RNnx = (gradRnx/gradRn)
RNny = (gradRny/gradRn)
pder = (dux(il,jl)+dvy(il,jl)-(divunx*rimdn1+divuny*rimdn2))
qr = q0*(gamma(il,jl)-1.0)*den(il,jl)*st5(il,jl)
Gn = gamma(il,jl)*p(il,jl)*pder + qr
ter = Gn/gradRn
c
RNnxr = gradRnxr/gradRnr
RNnyr = gradRnyr/gradRnr
pderr = (dux(ir,jr)+dvy(ir,jr)-(divunxr*rimdn1+divunyr*rimdn2))
qrr = q0*(gamma(ir,jr)-1.0)*den(ir,jr)*st5(ir,jr)
Gnr = gamma(ir,jr)*p(ir,jr)*pderr + qrr
terr = Gnr/gradRnr
c
uc = u(il,jl) + ss(il,jl)*rimdn1 + ter*RNnx
vc = v(il,jl) + ss(il,jl)*rimdn2 + ter*RNny
ucr = u(ir,jr) + ss(ir,jr)*rimdn1 + terr*RNnxr
vcr = v(ir,jr) + ss(ir,jr)*rimdn2 + terr*RNnyr
c
uchekn = uc*rnx+vc*rny
ucheknr = ucr*rnx+vcr*rny
c
If(lside.eq.1) Then
rck = 1.0
if(uchekn.lt.0d0.and.ucheknr.le.0D0) rck=-1.0
Else
rck =-1.0
if(ucheknr.gt.0d0.and.uchekn.ge.0D0) rck=1.0
Endif
c
IF(rck.ge.0.0) THEN
jj = il
kk = jl
dx = dxl
dy = dyl
ELSE
jj = ir
kk = jr
dx = dxr
dy = dyr
c
uc = ucr
vc = vcr
pder = pderr
RNnx = RNnxr
RNny = RNnyr
gradRn = gradRnr
ENDIF
c
dqrx = q0*(gamma(jj,kk)-1)*
& (ddenx(jj,kk)*st5(jj,kk) + den(jj,kk)*dst5x(jj,kk))
dqry = q0*(gamma(jj,kk)-1)*
& (ddeny(jj,kk)*st5(jj,kk) + den(jj,kk)*dst5y(jj,kk))
c
dGnx = dpx(jj,kk)*gamma(jj,kk)*pder + dqrx
dGny = dpy(jj,kk)*gamma(jj,kk)*pder + dqry
dterx = dGnx/gradRn
dtery = dGny/gradRn
c
durimx = dssx(jj,kk)*rimdn1 + dterx*RNnx
durimy = dssy(jj,kk)*rimdn1 + dtery*RNnx
dvrimx = dssx(jj,kk)*rimdn2 + dterx*RNny
dvrimy = dssy(jj,kk)*rimdn2 + dtery*RNny
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
a11 = dvy(jj,kk)+dvrimy+2.0/dt
a12 = -(duy(jj,kk)+durimy)
a21 = -(dvx(jj,kk)+dvrimx)
a22 = dux(jj,kk)+durimx+2.0/dt
b1 = 2.0*dx/dt-uc
b2 = 2.0*dy/dt-vc
det = a11*a22-a12*a21
c
xp = (a11*b1+a12*b2)/det
yp = (a21*b1+a22*b2)/det
If(abs(xp).gt.dxup.or.abs(yp).gt.dyup) go to 444
c
pup = p(jj,kk)+xp*dpx(jj,kk)+yp*dpy(jj,kk)
denup = den(jj,kk)+xp*ddenx(jj,kk)+yp*ddeny(jj,kk)
uup = u(jj,kk)+xp*dux(jj,kk)+yp*duy(jj,kk)
vup = v(jj,kk)+xp*dvx(jj,kk)+yp*dvy(jj,kk)
unup = uup*dn1+vup*dn2
c
uckp = dabs(u(jj,kk)+ss(jj,kk))
uckm = dabs(u(jj,kk)-ss(jj,kk))
vckp = dabs(v(jj,kk)+ss(jj,kk))
vckm = dabs(v(jj,kk)-ss(jj,kk))
uck=max(uckp,uckm)/dabs(dxcfl)
vck=max(vckp,vckm)/dabs(dxcfl)
If(ifc.eq.1.and.uck.gt.uchmax) uchmax=uck
If(ifc.eq.0.and.vck.gt.uchmax) uchmax=vck
c
Return
444 ifind = 0
Return
End
c
c------------------------------------------------------------------
c
subroutine isentrope(j,k,pf,denf,ifind,dt)
implicit double precision (a-h,o-z)
include 'rim2d.i'
common/Bltrackinfo/il,jl,ir,jr,dxl,dxr,
& dyl,dyr,rnx,rny,dxup,dyup,dxcfl
c
if(ifind.eq.0) go to 111
eps = 1d-18
c
gradR0x = dpx(il,jl)-ddenx(il,jl)*ss(il,jl)*ss(il,jl)
gradR0y = dpy(il,jl)-ddeny(il,jl)*ss(il,jl)*ss(il,jl)
gradR0 = dsqrt(gradR0x*gradR0x+gradR0y*gradR0y) + eps
G0 = q0*(gamma(il,jl)-1.0)*st5(il,jl)*den(il,jl)
c
gradR0xr = dpx(ir,jr)-ddenx(ir,jr)*ss(ir,jr)*ss(ir,jr)
gradR0yr = dpy(ir,jr)-ddeny(ir,jr)*ss(ir,jr)*ss(ir,jr)
gradR0r = dsqrt(gradR0xr*gradR0xr+gradR0yr*gradR0yr) + eps
G0r = q0*(gamma(ir,jr)-1.0)*st5(ir,jr)*den(ir,jr)
c
RNnx = gradR0x/gradR0
RNny = gradR0y/gradR0
RNnxr = gradR0xr/gradR0r
RNnyr = gradR0yr/gradR0r
c
u0 = G0*RNnx/gradR0
v0 = G0*RNny/gradR0
u0r = G0r*RNnxr/gradR0r
v0r = G0r*RNnyr/gradR0r
c
uc = u0 + u(il,jl)
vc = v0 + v(il,jl)
ucr = u0r + u(ir,jr)
vcr = v0r + v(ir,jr)
c
uchekn = uc*rnx+vc*rny
ucheknr = ucr*rnx+vcr*rny
c
If(uchekn.ge.0d0) Then
jj = il
kk = jl
dx = dxl
dy = dyl
Else
jj = ir
kk = jr
dx = dxr
dy = dyr
c
uc = ucr
vc = vcr
gradR0x = gradR0xr
gradR0y = gradR0yr
gradR0 = gradR0r
G0 = G0r
RNnx = RNnxr
RNny = RNnyr
ENDIF
c
dG0x = q0*(gamma(jj,kk)-1.0)*
& (st5(jj,kk)*ddenx(jj,kk)+dst5x(jj,kk)*den(jj,kk))
dG0y = q0*(gamma(jj,kk)-1.0)*
& (st5(jj,kk)*ddeny(jj,kk)+dst5y(jj,kk)*den(jj,kk))
c
dgradR0xdx = -2d0*ddenx(jj,kk)*ss(jj,kk)*dssx(jj,kk)
dgradR0xdy = -2d0*ddenx(jj,kk)*ss(jj,kk)*dssy(jj,kk)
dgradR0ydx = -2d0*ddeny(jj,kk)*ss(jj,kk)*dssx(jj,kk)
dgradR0ydy = -2d0*ddeny(jj,kk)*ss(jj,kk)*dssy(jj,kk)
c
ter3x=-G0*(gradR0x*dgradR0xdx+gradR0y*dgradR0ydx)/gradR0**2
ter3y=-G0*(gradR0x*dgradR0xdy+gradR0y*dgradR0ydy)/gradR0**2
c
du0x = (dG0x+ter3x)/gradR0*RNnx
du0y = (dG0y+ter3y)/gradR0*RNnx
dv0x = (dG0x+ter3x)/gradR0*RNny
dv0y = (dG0y+ter3y)/gradR0*RNny
c
a11 = dvy(jj,kk)+dv0y + 2.0/dt
a12 = -duy(jj,kk)+du0y
a21 = -dvx(jj,kk)+dv0x
a22 = dux(jj,kk) + du0x + 2.0/dt
b1 = 2.0*dx/dt-uc
b2 = 2.0*dy/dt-vc
det = a11*a22-a12*a21
c
xp = (a11*b1+a12*b2)/det
yp = (a21*b1+a22*b2)/det
If(abs(xp).gt.dxup.or.abs(yp).gt.dyup) go to 111
c
p0 = p(jj,kk) + xp*dpx(jj,kk) + yp*dpy(jj,kk)
den0 = den(jj,kk) + xp*ddenx(jj,kk) + yp*ddeny(jj,kk)
c
If(pf.le.p0) Then
denf = den0*(pf/p0)**(1.0/gamma(jj,kk))
Else
ter1 = (gamma(jj,kk)+1.0)+(gamma(jj,kk)-1.0)*pf/p0
ter2 = 0.5*p0/den0
denf = den0+(pf-p0)/(ter1*ter2)
Endif
Return
c
111 ifind =0
Return
End
c
c--------------------------------------------------------------
c
subroutine isorate(j,k,zf,dt)
implicit double precision (a-h,o-z)
include 'rim2d.i'
common/Bltrackinfo/il,jl,ir,jr,dxl,dxr,
& dyl,dyr,rnx,rny,dxup,dyup,dxcfl
c
eps = 1d-18
c
gradz = dzx(il,jl)*dzx(il,jl)+dzy(il,jl)*dzy(il,jl) + eps
gradzr = dzx(ir,jr)*dzx(ir,jr)+dzy(ir,jr)*dzy(ir,jr) + eps
c
uR =-st5(il,jl)*dzx(il,jl)/gradz
vR =-st5(il,jl)*dzy(il,jl)/gradz
uRr =-st5(ir,jr)*dzx(ir,jr)/gradzr
vRr =-st5(ir,jr)*dzy(ir,jr)/gradzr
c
uc = uR + u(il,jl)
vc = vR + v(il,jl)
ucr = uRr + u(ir,jr)
vcr = vRr + v(ir,jr)
c
uchekn = uc*rnx+vc*rny
ucheknr = ucr*rnx+vcr*rny
If(uchekn.ge.0.0.or.uchekn*ucheknr.le.0.0) Then
jj = il
kk = jl
dx = dxl
dy = dyl
ELSE
jj = ir
kk = jr
dx = dxr
dy = dyr
gradz = gradzr
uc = ucr
vc = vcr
ENDIF
c
duRx = -dst5x(jj,kk)*dzx(jj,kk)/gradz
duRy = -dst5y(jj,kk)*dzx(jj,kk)/gradz
dvRx = -dst5x(jj,kk)*dzy(jj,kk)/gradz
dvRy = -dst5y(jj,kk)*dzy(jj,kk)/gradz
c
a11 = dvy(jj,kk) + dvRy+2.0/dt
a12 = -(duy(jj,kk)+duRy)
a21 = -(dvx(jj,kk)+dvRx)
a22 = dux(jj,kk) + duRx+2.0/dt
b1 = 2.0*dx/dt-uc
b2 = 2.0*dy/dt-vc
det = a11*a22-a12*a21
c
xp = (a11*b1+a12*b2)/det
yp = (a21*b1+a22*b2)/det
zf = z(jj,kk)+xp*dzx(jj,kk)+yp*dzy(jj,kk)
c
If(dabs(xp).gt.dxup.or.dabs(yp).ge.dyup) Then
zf = z(jj,kk)+(dx*dzx(jj,kk)+dy*dzy(jj,kk)) +
& 0.5*dt*(st5(jj,kk)-
& u(jj,kk)*dzx(jj,kk)-v(jj,kk)*dzy(jj,kk))
Endif
c
Return
End
c
c--------------------------------------------------------------
c
subroutine tancom(j,k,utf,dt)
implicit double precision (a-h,o-z)
include 'rim2d.i'
common/Bldir/dn1,dn2,tandn1,tandn2
common/Bltrackinfo/il,jl,ir,jr,dxl,dxr,
& dyl,dyr,rnx,rny,dxup,dyup,dxcfl
c
eps = 1d-18
findchek = 0
c
gradRtx = (dux(il,jl)*tandn1+dvx(il,jl)*tandn2)
gradRty = (duy(il,jl)*tandn1+dvy(il,jl)*tandn2)
gradRt = dsqrt(gradRtx*gradRtx+gradRty*gradRty) + eps
Gt = (dpx(il,jl)*tandn1+dpy(il,jl)*tandn2)/den(il,jl)
c
gradRtxr = (dux(ir,jr)*tandn1+dvx(ir,jr)*tandn2)
gradRtyr = (duy(ir,jr)*tandn1+dvy(ir,jr)*tandn2)
gradRtr = dsqrt(gradRtxr*gradRtxr+gradRtyr*gradRtyr) + eps
Gtr = (dpx(ir,jr)*tandn1+dpy(ir,jr)*tandn2)/den(ir,jr)
c
RNvx = gradRtx/gradRt
RNvy = gradRty/gradRt
ut = Gt*RNvx/gradRt
vt = Gt*RNvy/gradRt
c
RNvxr = gradRtxr/gradRtr
RNvyr = gradRtyr/gradRtr
utr = Gtr*RNvxr/gradRtr
vtr = Gtr*RNvyr/gradRtr
c
222 findchek=findchek+1
c
uc = ut + u(il,jl)
vc = vt + v(il,jl)
ucr = utr + u(ir,jr)
vcr = vtr + v(ir,jr)
c
uchekn = uc*rnx+vc*rny
ucheknr = ucr*rnx+vcr*rny
If(uchekn.ge.0.0.or.uchekn*ucheknr.le.0.0) Then
jj = il
kk = jl
dx = dxl
dy = dyl
Else
jj = ir
kk = jr
dx = dxr
dy = dyr
ut = utr
vt = vtr
uc = ucr
vc = vcr
Endif
c
dutx = -ut*ddenx(jj,kk)/den(jj,kk)
duty = -ut*ddeny(jj,kk)/den(jj,kk)
dvtx = -vt*ddenx(jj,kk)/den(jj,kk)
dvty = -vt*ddeny(jj,kk)/den(jj,kk)
c
a11 = dvy(jj,kk)+dvty+2.0/dt
a12 = -(duy(jj,kk)+duty)
a21 = -(dvx(jj,kk)+dvtx)
a22 = dux(jj,kk)+dutx+2.0/dt
b1 = 2.0*dx/dt-uc
b2 = 2.0*dy/dt-vc
det = a11*a22-a12*a21
c
xp = (a11*b1+a12*b2)/det
yp = (a21*b1+a22*b2)/det
If((abs(xp).gt.dxup.or.abs(yp).gt.dyup).and.
& findchek.eq.1) Then
ut = 0.0
vt = 0.0
utr = 0.0
vtr = 0.0
go to 222
Endif
c
utan = u(jj,kk)+(xp*dux(jj,kk)+yp*duy(jj,kk))
vtan = v(jj,kk)+(xp*dvx(jj,kk)+yp*dvy(jj,kk))
utf = utan*tandn1+vtan*tandn2
c
Return
End
c
c--------------------------------------------------------------------
c
subroutine muscl(j,k,ifc,pl,denl,unl,pr,denr,unr,dt)
implicit double precision (a-h,o-z)
include 'rim2d.i'
common/Bldir/dn1,dn2,tandn1,tandn2
common/Bltrackinfo/il,jl,ir,jr,dxl,dxr,
& dyl,dyr,rnx,rny,dxup,dyup,dxcfl
u4l = u(il,jl) + ss(il,jl)*rnx
v4l = v(il,jl) + ss(il,jl)*rny
u4r = u(ir,jr) + ss(ir,jr)*rnx
v4r = v(ir,jr) + ss(ir,jr)*rny
c
dudx4l = dux(il,jl) + dssx(il,jl)*rnx
dudy4l = duy(il,jl) + dssy(il,jl)*rnx
dvdx4l = dvx(il,jl) + dssx(il,jl)*rny
dvdy4l = dvy(il,jl) + dssy(il,jl)*rny
c
dudx4r = dux(ir,jr) + dssx(ir,jr)*rnx
dudy4r = duy(ir,jr) + dssy(ir,jr)*rnx
dvdx4r = dvx(ir,jr) + dssx(ir,jr)*rny
dvdy4r = dvy(ir,jr) + dssy(ir,jr)*rny
c
slchekl = u4l*rnx + v4l*rny
slchekr = u4r*rnx + v4r*rny
If(slchekl.ge.0d0.or.slchekr.gt.0d0) Then
jj = il
kk = jl
u4 = u4l
v4 = v4l
dudx4 = dudx4l
dvdx4 = dvdx4l
dudy4 = dudy4l
dvdy4 = dvdy4l
dx = dxl
dy = dyl
Else
jj = ir
kk = jr
u4 = u4r
v4 = v4r
dudx4 = dudx4r
dvdx4 = dvdx4r
dudy4 = dudy4r
dvdy4 = dvdy4r
dx = dxr
dy = dyr
Endif
c
a11 = dvdy4+2.0/dt
a12 = -dudy4
a21 = -dvdx4
a22 = dudx4+2.0/dt
b1 = 2.0*dx/dt-u4
b2 = 2.0*dy/dt-v4
det = a11*a22-a12*a21
c
xs4 = (a11*b1+a12*b2)/det
ys4 = (a21*b1+a22*b2)/det
pl = p(jj,kk)+xs4*dpx(jj,kk)+ys4*dpy(jj,kk)
denl = den(jj,kk)+xs4*ddenx(jj,kk)+ys4*ddeny(jj,kk)
ul = u(jj,kk)+xs4*dux(jj,kk)+ys4*duy(jj,kk)
vl = v(jj,kk)+xs4*dvx(jj,kk)+ys4*dvy(jj,kk)
unl = ul*rnx+vl*rny
c
c
u1l = u(il,jl) - ss(il,jl)*rnx
v1l = v(il,jl) - ss(il,jl)*rny
u1r = u(ir,jr) - ss(ir,jr)*rnx
v1r = v(ir,jr) - ss(ir,jr)*rny
c
dudx1l = dux(il,jl) - dssx(il,jl)*rnx
dudy1l = duy(il,jl) - dssy(il,jl)*rnx
dvdx1l = dvx(il,jl) - dssx(il,jl)*rny
dvdy1l = dvy(il,jl) - dssy(il,jl)*rny
c
dudx1r = dux(ir,jr) - dssx(ir,jr)*rnx
dudy1r = duy(ir,jr) - dssy(ir,jr)*rnx
dvdx1r = dvx(ir,jr) - dssx(ir,jr)*rny
dvdy1r = dvy(ir,jr) - dssy(ir,jr)*rny
c
slchekl = u1l*rnx+v1l*rny
slchekr = u1r*rnx+v1r*rny
If(slchekr.le.0.0.or.slchekl.lt.0.0) Then
jj = ir
kk = jr
u1 = u1r
v1 = v1r
dudx1 = dudx1r
dvdx1 = dvdx1r
dudy1 = dudy1r
dvdy1 = dvdy1r
dx = dxr
dy = dyr
Else
jj = il
kk = jl
u1 = u1l
v1 = v1l
dudx1 = dudx1l
dvdx1 = dvdx1l
dudy1 = dudy1l
dvdy1 = dvdy1l
dx = dxl
dy = dyl
Endif
a11 = dvdy1+2.0/dt
a12 = -dudy1
a21 = -dvdx1
a22 = dudx1+2.0/dt
b1 = 2.0*dx/dt-u1
b2 = 2.0*dy/dt-v1
det = a11*a22-a12*a21
c
xs1 = (a11*b1+a12*b2)/det
ys1 = (a21*b1+a22*b2)/det
pr = p(jj,kk)+xs1*dpx(jj,kk)+ys1*dpy(jj,kk)
denr = den(jj,kk)+xs1*ddenx(jj,kk)+ys1*ddeny(jj,kk)
ur = u(jj,kk)+xs1*dux(jj,kk)+ys1*duy(jj,kk)
vr = v(jj,kk)+xs1*dvx(jj,kk)+ys1*dvy(jj,kk)
unr = ur*rnx+vr*rny
c
If(ifc.eq.1) Then
umax = max(abs(u4),abs(u1))/dxcfl
Else
umax = max(abs(v4),abs(v1))/dxcfl
endif
If(umax.gt.uchmax) uchmax = umax
c
Return
End
c
c--------------------------------------------------------------------
c
subroutine riemann(Gamma1,Gamma4,p1,p4,den1,den4,u1,u4,
& pf,uf,denfr,denfl,rmachr,rmachl,failriem,approx)
implicit double precision (a-h,o-z)
common/RFUNC/pp1,pp4,ssnd1,ssnd4,du
common/RGAM/G1,G4
external FF,GG,FFR,FFD,GGD,FFRD,SS,SSD
c
eps = 1d-10
G1 = Gamma1
G4 = Gamma4
snd1 = dsqrt(Gamma1*p1/den1)
snd4 = dsqrt(Gamma4*p4/den4)
pp1 = p1
pp4 = p4
ssnd1 = snd1
ssnd4 = snd4
du = u4-u1
c
c if(dabs(u1).lt.eps) u1 = 0d0
c if(dabs(u4).lt.eps) u4 = 0d0
c
c
c***********************************************************
c Right-facing shock and left-facing shock
c***********************************************************
c
If(du.lt.-eps) go to 20
c
r0=0.5d0*ABS(du)/snd1
Call NEWTON(GG,GGD,r0,r,fail,approx)
c
If(fail.eq.1.0) go to 20
If(r.lt.-eps.or.r.gt.du/snd1+eps) go to 20
c
rmachr=0.25*(Gamma1+1.)*r+dsqrt((0.25*(Gamma1+1.)*r)**2.+1.0)
pf=p1*(1.0+2.*Gamma1/(Gamma1+1.0)*(rmachr*rmachr-1.0))
denrat = (Gamma1+1.)*rmachr*rmachr /
& (2.0+(Gamma1-1.)*rmachr*rmachr)
denfr=den1*denrat
uf=r*snd1+u1
der=(u4-uf)/snd4
rmachl = 0.25*(Gamma4+1.)*der +
& dsqrt((0.25*(Gamma4+1.)*der)**2.+1.0)
denratl = (Gamma4+1.)*rmachl*rmachl /
& (2.0+(Gamma4-1.)*rmachl*rmachl)
denfl=den4*denratl
c
Return
c
20 If(p4-p1.lt.-eps) go to 30
c***********************************************************
c Right-facing shock and left-facing expansion
c***********************************************************
c
If(-du.ge.2.0*snd4/(Gamma4-1.0)+eps) go to 40
c
r0=0.5*ABS(du)/snd1
Call NEWTON(FF,FFD,r0,r,fail,approx)
If(fail.eq.1.0) then
If(p4.eq.p1) then
go to 30
else
go to 40
end if
end if
comp=du/snd1+2.0/(Gamma4-1.0)*snd4/snd1
If(r.lt.-eps.or.r.ge.comp+eps.or.r.lt.du/snd1-eps) then
If(p4.eq.p1) then
go to 30
else
go to 40
end if
end if
c
rmachr=0.25*(Gamma1+1.)*r+dsqrt((0.25*(Gamma1+1.)*r)**2.+1.)
pf=p1*( 1.0+2.*Gamma1/(Gamma1+1.0)*(rmachr*rmachr-1.) )
denrat=(Gamma1+1.)*rmachr*rmachr /
& (2.+(Gamma1-1.)*rmachr*rmachr)
uf=r*snd1+u1
denfr=den1*denrat
denfl=(pf/p4)**(1./Gamma4)*den4
rmachl=1.0
c
Return
c
c***********************************************************
c Right-facing expansion and left-facing shock
c***********************************************************
c
30 If(-du.ge.2.0*snd1/(Gamma1-1.0)+eps) go to 40
c
rs0=0.5*ABS(du)/snd4
Call NEWTON(FFR,FFRD,rs0,rs,fail,approx)
If(fail.eq.1.0) go to 40
c
comp=du/snd4+2.0/(Gamma1-1.0)*snd1/snd4
If(rs.lt.-eps.or.rs.ge.comp+eps.or.rs.lt.du/snd4-eps ) go to 40
c
rmachl=0.25*(Gamma4+1.)*rs+dsqrt( (0.25*(Gamma4+1.)*rs)**2.+1. )
pf=p4*( 1.0+2.*Gamma4/(Gamma4+1.0)*(rmachl*rmachl-1.) )
denratl=(Gamma4+1.)*rmachl*rmachl/( 2.+(Gamma4-1.)*rmachl*rmachl )
uf=-rs*snd4+u4
denfr=(pf/p1)**(1./Gamma1)*den1
denfl=den4*denratl
rmachr=1.0
c
Return
c
c
c***********************************************************
c Right-facing expansion and left-facing expansion
c***********************************************************
c
40 comp1=2.0*(snd1/(Gamma1-1.0)+snd4/(Gamma4-1.0))
If(-du.ge.comp1+eps.or.du.gt.eps) then
failriem=1.0
Return
end if
c
r0=0.5d0*dabs(du)/snd4
Call NEWTON(SS,SSD,r0,r,fail,approx)
if(abs(r).lt.eps) r=0d0
If(fail.eq.1.0) then
failriem=1.0
Return
end if
c
comp=du/snd1+2.0/(Gamma4-1.0)*snd4/snd1
If(r.gt.eps.or.r.le.-2.0/(Gamma1-1.)-eps
& .or.r.ge.comp+eps.or.r.lt.du/snd1-eps) then
failriem=1.0
Return
end if
c
pf = (1.0+0.5*(Gamma1-1.0)*r )**(2.*Gamma1/(Gamma1-1.))*p1
uf = r*snd1+u1
denfr = (pf/p1)**(1.0/Gamma1)*den1
denfl = (pf/p4)**(1.0/Gamma4)*den4
rmachr = 1d0
rmachl = 1d0
c
Return
End
c
c----------------------------------------------------------
c
Double Precision Function FF(r)
c
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum1=0.25*(Gamma1+1.0)
dum2=1.0+0.5*(Gamma4-1.0)*du/snd4
FF=1.0+Gamma1*dum1*r*r+Gamma1*r*dsqrt(dum1*dum1*r*r+1.)-p4/p1*
& ( dum2-.5*(Gamma4-1.)*snd1/snd4*r )**(2.*Gamma4/(Gamma4-1.))
c
Return
End
c
c----------------------------------------------------------
Double Precision Function FFD(r)
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum1=0.25*(Gamma1+1.0)
dum2=1.0+0.5*(Gamma4-1.0)*du/snd4
dum3=dsqrt(dum1*dum1*r*r+1.0)
FFD=2.*Gamma1*dum1*r+Gamma1*dum3+Gamma1*dum1*dum1/dum3*r*r
& +Gamma4*p4/p1*snd1/snd4*( dum2-.5*(Gamma4-1.)*snd1/snd4*r )**
& ( (Gamma4+1.0)/(Gamma4-1.0) )
c
Return
End
c
c----------------------------------------------------------
Double Precision Function GG(r)
c
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum1 = 0.25*(Gamma1+1.0)
dum2 = du/snd4-snd1/snd4*r
dum3 = dsqrt(dum1*dum1*r*r+1.0)
dum4 = 0.25*(Gamma4+1.0)
GG=1.0-p4/p1+Gamma1*dum1*r*r+Gamma1*r*dum3
& -p4/p1*Gamma4*dum4*dum2*dum2-p4/p1*Gamma4*dum2*
& dsqrt( dum4*dum4*dum2*dum2+1.0 )
c
Return
End
c
c----------------------------------------------------------
Double Precision Function GGD(r)
c
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum1 = 0.25*(Gamma1+1.0)
dum2 = du/snd4-snd1/snd4*r
dum3 = dsqrt(dum1*dum1*r*r+1.0)
dum4 = 0.25*(Gamma4+1.0)
GGD=2.*Gamma1*dum1*r+Gamma1*dum3+Gamma1*dum1*dum1/dum3*r*r
& + p4/p1*snd1/snd4*Gamma4*dsqrt(dum4*dum4*dum2*dum2+1.0)
& + p4/p1*snd1/snd4*Gamma4*2.*dum4*dum2
& + p4/p1*snd1/snd4*Gamma4*dum4*dum4*dum2*dum2
& /dsqrt( dum4*dum4*dum2*dum2+1.0 )
c
Return
End
c
c----------------------------------------------------------
Double Precision Function FFR(r)
c
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum1=0.25*(Gamma4+1.0)
dum2=1.0+0.5*(Gamma1-1.0)*du/snd1
FFR=1.0+Gamma4*dum1*r*r+Gamma4*r*dsqrt(dum1*dum1*r*r+1.)-p1/p4*
& (dum2-.5*(Gamma1-1.0)*snd4/snd1*r)**(2.0*Gamma1/(Gamma1-1.0))
c
Return
End
c
c----------------------------------------------------------
Double Precision Function FFRD(r)
c
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum1=0.25*(Gamma4+1.0)
dum2=1.0+0.5*(Gamma1-1.0)*du/snd1
dum3=dsqrt(dum1*dum1*r*r+1.0)
FFRD=2.*Gamma4*dum1*r+Gamma4*dum3+Gamma4*dum1*dum1/dum3*r*r
& +Gamma1*p1/p4*snd4/snd1*( dum2-.5*(Gamma1-1.)*snd4/snd1*r )
& **( (Gamma1+1.0)/(Gamma1-1.0) )
c
Return
End
c
c----------------------------------------------------------
Double Precision Function SS(r)
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
c
dum = 1d0+0.5d0*(Gamma4-1d0)*du/snd4
SS = (1d0+0.5d0*(Gamma1-1d0)*r )**(2d0*Gamma1/(Gamma1-1d0))-
& p4/p1*( dum-0.5d0*(Gamma4-1d0)*snd1/snd4*r )**
& (2d0*Gamma4/(Gamma4-1d0))
c
Return
End
c
c----------------------------------------------------------
Double Precision Function SSD(r)
c
Implicit Double Precision (a-h,o-z)
Common/RFUNC/p1,p4,snd1,snd4,du
Common/RGAM/Gamma1,Gamma4
dum=1.0+0.5*(Gamma4-1.0)*du/snd4
SSD=Gamma1*( 1.0+0.5*(Gamma1-1.0)*r )
& **((Gamma1+1.0)/(Gamma1-1.0))+
& Gamma4*p4/p1*snd1/snd4*( dum-0.5*(Gamma4-1.0)*snd1/snd4*r )
& **((Gamma4+1.0)/(Gamma4-1.0))
c
Return
End
c
c----------------------------------------------------------
c
Subroutine Newton(F,FD,x0,x,fail,approx)
c
Implicit Double Precision (a-h,o-z)
External F,FD
c
fail=0.0
x=x0
Do i=1,20
comp = x
x = x-F(x)/FD(x)
If(dabs(x-comp).lt.approx) Return
End do
fail=1.0
c
Return
End
c
c----------------------------------------------------------
c
Subroutine Denface(Gammar,Gammal,pr,pl,denr,denl,ur,ul,rmachr,
& rmachl,pf,uf,denfr,denfl,denf)
Implicit Double Precision (a-h,o-z)
c
sndr=DSQRT(Gammar*pr/denr)
sndl=DSQRT(Gammal*pl/denl)
Vr=ur+rmachr*sndr
Vl=ul-rmachl*sndl
c
If(pf.ge.pl) then
Vlstar=Vl
else
Vlstar=uf-DSQRT(Gammal*pf/denfl)
end if
If(pf.ge.pr) then
Vrstar=Vr
else
Vrstar=uf+DSQRT(Gammar*pf/denfr)
end if
c
denf=denfl
If(Vlstar.le.0.0.and.uf.gt.0.0) RETURN
denf=denfr
If(Vrstar.ge.0.0.and.uf.le.0.0) RETURN
c
uf=( (Gammar-1.0)*ur-2.*sndr )/(Gammar+1.0)
dum=1.0+0.5*(Gammar-1.0)*(uf-ur)/sndr
If(dum.le.0.0) GO TO 20
pf=pr*dum**(2.0*Gammar/(Gammar-1.0))
c
denf=denr*(pf/pr)**(1.0/Gammar)
If(Vrstar.lt.0.0.and.Vr.gt.0.0) RETURN
c
20 uf=( (Gammal-1.0)*ul+2.*sndl )/(Gammal+1.0)
dum=1.0+0.5*(Gammal-1.0)*(ul-uf)/sndl
If(dum.le.0.0) GO TO 30
pf=pl*dum**(2.0*Gammal/(Gammal-1.0))
denf=denl*(pf/pl)**(1.0/Gammal)
If(Vlstar.gt.0.0.and.Vl.lt.0.0) RETURN
c
30 uf=ur
pf=pr
denf=denr
If(Vr.le.0.0) RETURN
c
uf=ul
pf=pl
denf=denl
c
RETURN
END
c----------------------------------------------------------------