c c ========================================================= subroutine rpn3euznd(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,maux, & auxl,auxr,wave,s,amdq,apdq) c ========================================================= c c # Hybrid Roe-solver for the 2D ZND-Euler equations c # solve Riemann problems along one slice of data. c # Scheme is blended with HLL for robustness. 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, c # s the speeds, c # amdq the left-going flux difference A^- \Delta q c # apdq the right-going flux difference A^+ \Delta q 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: Ralf Deiterding c implicit double precision (a-h,o-z) c 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 apdq(1-mbc:maxm+mbc, meqn) dimension amdq(1-mbc:maxm+mbc, meqn) dimension auxl(1-mbc:maxm+mbc, maux, 3) dimension auxr(1-mbc:maxm+mbc, maux, 3) c c local arrays -- common block comroe is passed to rpt3euznd c ------------ parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5 parameter (minmrp = -4) !# assumes at most mbc=5 dimension delta(6), fl(minmrp:maxmrp,6), fr(minmrp:maxmrp,6) logical efix, pfix, hll, roe, hllfix common /param/ gamma,gamma1,q0 common /comroe/ u2v2w2(minmrp:maxmrp), & u(minmrp:maxmrp),v(minmrp:maxmrp),w(minmrp:maxmrp), & enth(minmrp:maxmrp),a(minmrp:maxmrp),Y(2,minmrp:maxmrp) c data efix /.true./ !# use entropy fix for transonic rarefactions data pfix /.true./ !# use Larrouturou's positivity fix for species data hll /.true./ !# use HLL instead of Roe solver, if unphysical values occur data roe /.true./ !# turn off Roe solver when debugging HLL 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 momentums: c if(ixyz .eq. 1)then mu = 3 mv = 4 mw = 5 else if(ixyz .eq. 2)then mu = 4 mv = 5 mw = 3 else mu = 5 mv = 3 mw = 4 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 rpt3euznd to do the transverse wave splitting. c do 10 i=2-mbc,mx+mbc c pl = gamma1*(qr(i-1,6) - qr(i-1,2)*q0 - & 0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2+qr(i-1,mw)**2)/ & (qr(i-1,1)+qr(i-1,2))) pr = gamma1*(ql(i, 6) - ql(i, 2)*q0 - & 0.5d0*(ql(i, mu)**2+ql(i, mv)**2+ql(i, mw)**2)/ & (ql(i, 1)+ql(i, 2))) rhsqrtl = dsqrt(qr(i-1,1) + qr(i-1,2)) rhsqrtr = dsqrt(ql(i, 1) + ql(i, 2)) 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 u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2 enth(i) = (((qr(i-1,6)+pl)/rhsqrtl & + (ql(i ,6)+pr)/rhsqrtr)) / rhsq2 Y(1,i) = (qr(i-1,1)/rhsqrtl + ql(i,1)/rhsqrtr) / rhsq2 Y(2,i) = (qr(i-1,2)/rhsqrtl + ql(i,2)/rhsqrtr) / rhsq2 c # speed of sound a2 = gamma1*(enth(i) - 0.5d0*u2v2w2(i) - Y(2,i)*q0) a(i) = dsqrt(a2) c 10 continue c do 30 i=2-mbc,mx+mbc c c # find a1 thru a5, the coefficients of the 5 eigenvectors: c do k = 1, 6 delta(k) = ql(i,k) - qr(i-1,k) enddo drho = delta(1) + delta(2) c a2 = gamma1/a(i)**2 * (drho*0.5d0*u2v2w2(i) - delta(2)*q0 & - (u(i)*delta(mu)+v(i)*delta(mv)+w(i)*delta(mw)) + & delta(6)) a3 = delta(mv) - v(i)*drho a4 = delta(mw) - w(i)*drho a5 = 0.5d0*( a2 - ( u(i)*drho - delta(mu) )/a(i) ) a1 = a2 - a5 c c # Compute the waves. c c # 1-wave wave(i,1,1) = a1*Y(1,i) wave(i,2,1) = a1*Y(2,i) 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,6,1) = a1*(enth(i) - u(i)*a(i)) s(i,1) = u(i)-a(i) c c # 2-wave wave(i,1,2) = delta(1) - Y(1,i)*a2 wave(i,2,2) = delta(2) - Y(2,i)*a2 wave(i,mu,2) = (drho - a2)*u(i) wave(i,mv,2) = (drho - a2)*v(i) + a3 wave(i,mw,2) = (drho - a2)*w(i) + a4 wave(i,6,2) = (drho - a2)*0.5d0*u2v2w2(i) + & q0*(delta(2) - Y(2,i)*a2) + a3*v(i) + a4*w(i) s(i,2) = u(i) c c # 3-wave wave(i,1,3) = a5*Y(1,i) wave(i,2,3) = a5*Y(2,i) 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,6,3) = a5*(enth(i) + u(i)*a(i)) s(i,3) = u(i)+a(i) c 30 continue c call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,auxr,apdq) call flx3(ixyz,maxm,meqn,mbc,mx,ql,maux,auxl,amdq) c do 35 i = 1-mbc, mx+mbc do 35 m=1,meqn fl(i,m) = amdq(i,m) fr(i,m) = apdq(i,m) 35 continue c if (.not.roe) go to 900 c c # compute flux differences amdq and apdq. c --------------------------------------- c if (efix) go to 110 c c # no entropy fix c ---------------- c c # amdq = SUM s*wave over left-going waves c # apdq = SUM s*wave over right-going waves c do 100 m=1,meqn do 100 i=2-mbc, mx+mbc amdq(i,m) = 0.d0 apdq(i,m) = 0.d0 do 90 mws=1,mwaves if (s(i,mws) .lt. 0.d0) then amdq(i,m) = amdq(i,m) + s(i,mws)*wave(i,m,mws) endif 90 continue 100 continue go to 900 110 continue c c # With entropy fix c ------------------ c c # compute flux differences amdq and apdq. c # First compute amdq as sum of s*wave for left going waves. c # Incorporate entropy fix by adding a modified fraction of wave c # if s should change sign. c do 200 i=2-mbc,mx+mbc c c # check 1-wave: c --------------- c rk1 = qr(i-1,1) rk2 = qr(i-1,2) rhou = qr(i-1,mu) rhov = qr(i-1,mv) rhow = qr(i-1,mw) rhoE = qr(i-1,6) rho = rk1 + rk2 p = gamma1*(rhoE - rk2*q0 - 0.5d0*(rhou**2+rhov**2+rhow**2)/rho) if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200 c = dsqrt(gamma*p/rho) s0 = rhou/rho - c !# u-c in left state (cell i-1) * write(6,*) 'left state 0', a(i), c, T c c # check for fully supersonic case: if (s0.ge.0.d0 .and. s(i,1).gt.0.d0) then c # everything is right-going do 60 m=1,meqn amdq(i,m) = 0.d0 60 continue go to 200 endif c rk1 = rk1 + wave(i,1,1) rk2 = rk2 + wave(i,2,1) rhou = rhou + wave(i,mu,1) rhov = rhov + wave(i,mv,1) rhow = rhow + wave(i,mw,1) rhoE = rhoE + wave(i,6,1) rho = rk1 + rk2 p = gamma1*(rhoE - rk2*q0 - 0.5d0*(rhou**2+rhov**2+rhow**2)/rho) if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200 c = dsqrt(gamma*p/rho) s1 = rhou/rho - c !# u-c to right of 1-wave * write(6,*) 'left state 1', a(i), c, T c if (s0.lt.0.d0 .and. s1.gt.0.d0) then c # transonic rarefaction in the 1-wave sfract = s0 * (s1-s(i,1)) / (s1-s0) else if (s(i,1) .lt. 0.d0) then c # 1-wave is leftgoing sfract = s(i,1) else c # 1-wave is rightgoing sfract = 0.d0 !# this shouldn't happen since s0 < 0 endif do 120 m=1,meqn amdq(i,m) = sfract*wave(i,m,1) 120 continue c c # check 2-wave: c --------------- c if (s(i,2) .ge. 0.d0) go to 200 !# 2-wave is rightgoing do 140 m=1,meqn amdq(i,m) = amdq(i,m) + s(i,2)*wave(i,m,2) 140 continue c c # check 3-wave: c --------------- c rk1 = ql(i,1) rk2 = ql(i,2) rhou = ql(i,mu) rhov = ql(i,mv) rhow = ql(i,mw) rhoE = ql(i,6) rho = rk1 + rk2 p = gamma1*(rhoE - rk2*q0 - 0.5d0*(rhou**2+rhov**2+rhow**2)/rho) if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200 c = dsqrt(gamma*p/rho) s3 = rhou/rho + c !# u+c in right state (cell i) * write(6,*) 'right state 1', a(i), c, T c rk1 = rk1 - wave(i,1,3) rk2 = rk2 - wave(i,2,3) rhou = rhou - wave(i,mu,3) rhov = rhov - wave(i,mv,3) rhow = rhow - wave(i,mw,3) rhoE = rhoE - wave(i,6,3) rho = rk1 + rk2 p = gamma1*(rhoE - rk2*q0 - 0.5d0*(rhou**2+rhov**2+rhow**2)/rho) if ((rho.le.0.d0.or.p.le.0.d0).and.hll) go to 200 c = dsqrt(gamma*p/rho) s2 = rhou/rho + c !# u+c to left of 3-wave * write(6,*) 'right state 0', a(i), c, T c if (s2 .lt. 0.d0 .and. s3.gt.0.d0) then c # transonic rarefaction in the 3-wave sfract = s2 * (s3-s(i,3)) / (s3-s2) else if (s(i,3) .lt. 0.d0) then c # 3-wave is leftgoing sfract = s(i,3) else c # 3-wave is rightgoing go to 200 endif c do 160 m=1,meqn amdq(i,m) = amdq(i,m) + sfract*wave(i,m,3) 160 continue 200 continue c 900 continue c if (hll) then do 350 i = 2-mbc, mx+mbc hllfix = .false. if (.not.roe) hllfix = .true. c rho1l = qr(i-1,1) + wave(i,1,1) rho2l = qr(i-1,2) + wave(i,2,1) rhoul = qr(i-1,mu) + wave(i,mu,1) rhovl = qr(i-1,mv) + wave(i,mv,1) rhowl = qr(i-1,mw) + wave(i,mw,1) El = qr(i-1,6) + wave(i,6,1) pl = gamma1*(El - rho2l*q0 - & 0.5d0*(rhoul**2 + rhovl**2 + rhowl**2)/(rho1l+rho2l)) if (rho1l+rho2l.le.0.d0.or.pl.le.0.d0) & hllfix = .true. c rho1r = ql(i,1) - wave(i,1,3) rho2r = ql(i,2) - wave(i,2,3) rhour = ql(i,mu) - wave(i,mu,3) rhovr = ql(i,mv) - wave(i,mv,3) rhowr = ql(i,mw) - wave(i,mw,3) Er = ql(i,5 ) - wave(i,5,3) pr = gamma1*(Er - rho2r*q0 - & 0.5d0*(rhour**2 + rhovr**2 + rhowr**2)/(rho1r+rho2r)) if (rho1r+rho2r.le.0.d0.or.pr.le.0.d0) & hllfix = .true. c if (hllfix) then * if (roe) write (6,*) 'Switching to HLL in',i c rl = qr(i-1,1) + qr(i-1,2) ul = qr(i-1,mu)/rl pl = gamma1*(qr(i-1,6) - qr(i-1,2)*q0 - 0.5d0* & (qr(i-1,mu)**2+qr(i-1,mv)**2+qr(i-1,mw)**2)/rl) al = dsqrt(gamma*pl/rl) c rr = ql(i ,1) + ql(i ,2) ur = ql(i ,mu)/rr pr = gamma1*(ql(i ,6) - ql(i ,2)*q0 - 0.5d0* & (ql(i ,mu)**2+ql(i ,mv)**2+ql(i ,mw)**2)/rr) ar = dsqrt(gamma*pr/rr) c sl = dmin1(ul-al,ur-ar) sr = dmax1(ul+al,ur+ar) c do m=1,meqn if (sl.ge.0.d0) fg = fr(i-1,m) if (sr.le.0.d0) fg = fl(i,m) if (sl.lt.0.d0.and.sr.gt.0.d0) & fg = (sr*fr(i-1,m) - sl*fl(i,m) + & sl*sr*(ql(i,m)-qr(i-1,m)))/ (sr-sl) amdq(i,m) = fg-fr(i-1,m) enddo s(i,1) = sl s(i,2) = 0.d0 s(i,3) = sr endif 350 continue endif c if (pfix) then do 300 i=2-mbc,mx+mbc amdr = amdq(i,1)+amdq(i,2) rhol = qr(i-1,1)+qr(i-1,2) rhor = ql(i ,1)+ql(i ,2) do 300 m=1,2 if (qr(i-1,mu)+amdr.gt.0.d0) then Z = qr(i-1,m)/rhol else Z = ql(i ,m)/rhor endif amdq(i,m) = Z*amdr + (Z-qr(i-1,m)/rhol)*qr(i-1,mu) 300 continue endif c do 400 i = 2-mbc, mx+mbc do 400 m=1,meqn amdq(i,m) = fr(i-1,m) + amdq(i,m) apdq(i,m) = -amdq(i,m) 400 continue c return end c