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----------------------------------------------------------------