c c c ===================================================== subroutine rpn2eu(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,maux, & auxl,auxr,wave,s,fl,fr) c ===================================================== c c # solve Riemann problems for the 2D Euler equations using c # an improved version of the Liou-Steffen Flux-Vector-Splitting c c # Yasuhiro Wada, Meng-Sing Liou "An accurate and robust flux c # splitting scheme for shock and contact discontinuities", c # SIAM J. Sci. Comput., Vol. 18, No.2, pp 633-657, May 1997. 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 # On output, wave contains the waves, c # s the speeds, fl and fr the positive and negative flux. 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, rp is called with ql = qr = q. c c # Copyright (C) 2002 Ralf Deiterding c # Brandenburgische Universitaet Cottbus c implicit double precision (a-h,o-z) 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 fl(1-mbc:maxm+mbc, meqn) dimension fr(1-mbc:maxm+mbc, meqn) double precision l(4), r(4) common /param/ gamma,gamma1 c c # Method returns fluxes c ------------ common /rpnflx/ mrpnflx mrpnflx = 1 c c # set mu to point to the component of the system that corresponds c # to momentum in the direction of this slice, mv to the orthogonal c # momentum: c if (ixy.eq.1) then mu = 2 mv = 3 else mu = 3 mv = 2 endif c c # AUSM Flux Vector Splitting c do 10 i=2-mbc,mx+mbc rhol = qr(i-1,1) rhor = ql(i ,1) ul = qr(i-1,mu)/rhol ur = ql(i ,mu)/rhor vl = qr(i-1,mv)/rhol vr = ql(i ,mv)/rhor pl = gamma1*(qr(i-1,4) - 0.5d0*(qr(i-1,mu)**2+ & qr(i-1,mv)**2)/rhol) pr = gamma1*(ql(i ,4) - 0.5d0*(ql(i ,mu)**2+ & ql(i ,mv)**2)/rhor) Hl = (qr(i-1,4)+pl)/rhol Hr = (ql(i ,4)+pr)/rhor al = dsqrt(gamma*pl/rhol) ar = dsqrt(gamma*pr/rhor) c am = dmax1(al,ar) alphal = 2.d0*(pl/rhol)/(pl/rhol+pr/rhor) alphar = 2.d0*(pr/rhor)/(pl/rhol+pr/rhor) c ulp = 0.5d0*(ul+dabs(ul)) plp = pl*ulp/ul if (dabs(ul).le.am) then ulp = 0.25d0*alphal*(ul+al)**2/am + (1.d0-alphal)*ulp plp = 0.25d0*pl*(ul+al)**2/am**2*(2.d0-ul/am) endif c urm = 0.5d0*(ur-dabs(ur)) prm = pr*urm/ur if (dabs(ur).le.am) then urm = -0.25d0*alphar*(ur-ar)**2/am + (1.d0-alphar)*urm prm = 0.25d0*pr*(ur-ar)**2/am**2*(2.d0+ur/am) endif c c # Blending between AUSMV and AUSMD c # sf=1.d0 gives AUSMV, sf=-1.d0 gives AUSMD sf = dmin1(1.d0, 10.d0*dabs(pr-pl)/dmin1(pl,pr)) c l(1) = 0.5d0*(ulp*rhol+dabs(ulp*rhol)) l(mu) = 0.5d0*((1.d0+sf)*ulp*rhol*ul + (1.d0-sf)*l(1)*ul) + plp l(mv) = l(1)*vl l(4) = l(1)*Hl c r(1) = 0.5d0*(urm*rhor-dabs(urm*rhor)) r(mu) = 0.5d0*((1.d0+sf)*urm*rhor*ur + (1.d0-sf)*r(1)*ur) + prm r(mv) = r(1)*vr r(4) = r(1)*Hr c do 20 m = 1,meqn fl(i,m) = l(m) + r(m) fr(i,m) = -fl(i,m) 20 continue c s(i,1) = 0.5d0*(ul-al + ur-ar) s(i,2) = 0.5d0*(ul + ur) s(i,3) = 0.5d0*(ul+al + ur+ar) do 10 mw=1,mwaves do 10 m=1,meqn wave(i,m,mw) = 0.d0 10 continue c return end c c