c c c ===================================================== subroutine rpn3eu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr, & maux,auxl,auxr,wave,s,fl,fr) c ===================================================== c c # Riemann solver for the 3D Euler equations c # The waves are computed using the Roe approximation. c c # This is quite a bit slower than the Roe solver, c # but may give more accurate solutions for some problems. 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 # 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 output, wave contains the waves, s the speeds c # and fl, fr the positive and negative Godunov 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, this routine is called with ql = qr c c Author: Randall J. LeVeque 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 auxl(1-mbc:maxm+mbc, maux, 3) dimension auxr(1-mbc:maxm+mbc, maux, 3) dimension fr(1-mbc:maxm+mbc, meqn) dimension fl(1-mbc:maxm+mbc, meqn) c c local arrays -- common block comroe is passed to rpt3eu c ------------ parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5 parameter (minmrp = -4) !# assumes at most mbc=5 dimension delta(5) dimension sl(2),sr(2) common /param/ gamma,gamma1 common /comroe/ u2v2w2(minmrp:maxmrp), & u(minmrp:maxmrp),v(minmrp:maxmrp),w(minmrp:maxmrp), & enth(minmrp:maxmrp),a(minmrp:maxmrp),g1a2(minmrp:maxmrp), & euv(minmrp:maxmrp) c c # Riemann solver returns flux differences c ------------ common /rpnflx/ mrpnflx mrpnflx = 1 c if (minmrp.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 momentum: c if(ixyz .eq. 1)then mu = 2 mv = 3 mw = 4 else if(ixyz .eq. 2)then mu = 3 mv = 4 mw = 2 else mu = 4 mv = 2 mw = 3 endif 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 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 splitting. c do 10 i = 2-mbc, mx+mbc rhsqrtl = dsqrt(qr(i-1,1)) rhsqrtr = dsqrt(ql(i,1)) pl = gamma1*(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)) pr = gamma1*(ql(i,5) - 0.5d0*(ql(i,mu)**2 + & ql(i,mv)**2 + ql(i,mw)**2)/ql(i,1)) rhsq2 = rhsqrtl + rhsqrtr 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 q 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) a4 = g1a2(i) * (euv(i)*delta(1) & + u(i)*delta(2) + v(i)*delta(3) + w(i)*delta(4) & - delta(5)) 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 c c # Compute the waves. c # Note that the 2-wave, 3-wave and 4-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)) 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) 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)) s(i,3) = u(i)+a(i) 20 continue c c c # compute Godunov flux f0 at each interface. c # Uses exact Riemann solver c c do 200 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 wl = qr(i-1,mw)/rhol wr = ql(i ,mw)/rhor pl = gamma1*(qr(i-1,5)-0.5*(qr(i-1,mu)*ul + qr(i-1,mv)*vl + & qr(i-1,mw)*wl )) pr = gamma1*(ql(i ,5)-0.5*(ql(i ,mu)*ur + ql(i ,mv)*vr + & ql(i ,mw)*wr)) c c # iterate to find pstar, ustar: c alpha = 1. pstar = 0.5*(pl+pr) wsr = dsqrt(pr*rhor) * phi(pstar/pr) wsl = dsqrt(pl*rhol) * phi(pstar/pl) c if (pl.eq.pr .and. rhol.eq.rhor) go to 60 c 40 do 50 iter=1,1000 p1 = (ul-ur+pr/wsr+pl/wsl) / (1./wsr + 1./wsl) pstar = dmax1(p1,1d-6)*alpha + (1.-alpha)*pstar wr1 = wsr wl1 = wsl wsr = dsqrt(pr*rhor) * phi(pstar/pr) wsl = dsqrt(pl*rhol) * phi(pstar/pl) if (dmax1(abs(wr1-wsr),dabs(wl1-wsl)) .lt. 1d-6) & go to 60 50 continue c c # nonconvergence: alpha = alpha/2. if (alpha .gt. 0.1) go to 40 write(6,*) 'no convergence',i,wr1,wsr,wl1,wsl wsr = .5*(wsr+wr1) wsl = .5*(wsl+wl1) c 60 continue ustar = (pl-pr+wsr*ur+wsl*ul) / (wsr+wsl) c c c # left wave: c ============ c if (pstar .gt. pl) then c c # shock: sl(1) = ul - wsl/rhol sr(1) = sl(1) rho1 = wsl/(ustar-sl(1)) c else c c # rarefaction: cl = dsqrt(gamma*pl/rhol) cstar = cl + 0.5*gamma1*(ul-ustar) sl(1) = ul-cl sr(1) = ustar-cstar rho1 = (pstar/pl)**(1./gamma) * rhol endif c c c c # right wave: c ============= c if (pstar .ge. pr) then c c # shock sl(2) = ur + wsr/rhor sr(2) = sl(2) rho2 = wsr/(sl(2)-ustar) c else c c # rarefaction: cr = dsqrt(gamma*pr/rhor) cstar = cr + 0.5*gamma1*(ustar-ur) sr(2) = ur+cr sl(2) = ustar+cstar rho2 = (pstar/pr)**(1./gamma)*rhor endif c c c # compute flux: c =============== c c # compute state (rhos,us,ps) at x/t = 0: c if (sl(1).gt.0) then rhos = rhol us = ul vs = vl ws = wl ps = pl else if (sr(1).le.0. .and. ustar.ge. 0.) then rhos = rho1 us = ustar vs = vl ws = wl ps = pstar else if (ustar.lt.0. .and. sl(2).ge. 0.) then rhos = rho2 us = ustar vs = vr ws = wr ps = pstar else if (sr(2).lt.0) then rhos = rhor us = ur vs = vr ws = wr ps = pr else if (sl(1).le.0. .and. sr(1).ge.0.) then c # transonic 1-rarefaction us = (gamma1*ul + 2.*cl)/(gamma+1.) e0 = pl/(rhol**gamma) rhos = (us**2/(gamma*e0))**(1./gamma1) ps = e0*rhos**gamma vs = vl ws = wl else if (sl(2).le.0. .and. sr(2).ge.0.) then c # transonic 3-rarefaction us = (gamma1*ur - 2.*cr)/(gamma+1.) e0 = pr/(rhor**gamma) rhos = (us**2/(gamma*e0))**(1./gamma1) ps = e0*rhos**gamma vs = vr ws = wr endif c fl(i,1) = rhos*us fl(i,mu) = rhos*us**2 + ps fl(i,mv) = rhos*us*vs fl(i,mw) = rhos*us*ws fl(i,5) = us*(gamma*ps/gamma1 + 0.5*rhos*(us**2+vs**2+ws**2)) 200 continue c do 220 m=1,6 do 220 i = 2-mbc, mx+mbc fr(i,m) = -fl(i,m) 220 continue c return end c c c double precision function phi(w) implicit double precision (a-h,o-z) common/param/ gamma,gamma1 c sqg = dsqrt(gamma) if (w .gt. 1.) then phi = dsqrt(w*(gamma+1.)/2. + gamma1/2.) else if (w .gt. 0.99999) then phi = sqg else if (w .gt. .999) then phi = sqg + (2*gamma**2 - 3.*gamma + 1) & *(w-1.) / (4.*sqg) else phi = gamma1*(1.-w) / (2.*sqg*(1.-w**(gamma1/(2.*gamma)))) endif return end