c c ========================================================= subroutine rp1euznd(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux, & auxl,auxr,wave,s,amdq,apdq) c ========================================================= c c # solve Riemann problems for the 1D ZND-Euler equations using Roe's c # approximate Riemann solver. 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 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 # Copyright (C) 2002 Ralf Deiterding c # Brandenburgische Universitaet Cottbus c implicit double precision (a-h,o-z) c dimension wave(1-mbc:maxmx+mbc, meqn, mwaves) dimension s(1-mbc:maxmx+mbc, mwaves) dimension ql(1-mbc:maxmx+mbc, meqn) dimension qr(1-mbc:maxmx+mbc, meqn) dimension auxl(1-mbc:maxmx+mbc, maux) dimension auxr(1-mbc:maxmx+mbc, maux) dimension apdq(1-mbc:maxmx+mbc, meqn) dimension amdq(1-mbc:maxmx+mbc, meqn) c c # local storage c --------------- parameter (max2 = 100002) !# assumes at most 100000 grid points with mbc=2 dimension u(-1:max2), enth(-1:max2), a(-1:max2), smax(-1:max2) dimension delta(4), Y(2,-1:max2), fl(-1:max2,4), fr(-1:max2,4) logical efix, pfix, hll, roe, hllfix common /param/ gamma,gamma1,q0 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 solver if unphysical values occur data roe /.true./ !# use Roe solver c c # Riemann solver returns flux differences c ------------ common /rpnflx/ mrpnflx mrpnflx = 1 c if (-1.gt.1-mbc .or. max2 .lt. maxmx+mbc) then write(6,*) 'need to increase max2 in rp' stop endif c c # Compute Roe-averaged quantities: c do 10 i=2-mbc,mx+mbc c rhol = qr(i-1,1)+qr(i-1,2) rhor = ql(i ,1)+ql(i ,2) ul = qr(i-1,3)/rhol ur = ql(i ,3)/rhor pl = gamma1*(qr(i-1,4) - qr(i-1,2)*q0 - 0.5d0*ul**2*rhol) pr = gamma1*(ql(i, 4) - ql(i, 2)*q0 - 0.5d0*ur**2*rhor) al = dsqrt(gamma*pl/rhol) ar = dsqrt(gamma*pr/rhor) rhsqrtl = dsqrt(rhol) rhsqrtr = dsqrt(rhor) rhsq2 = rhsqrtl + rhsqrtr u(i) = (qr(i-1,3)/rhsqrtl + ql(i,3)/rhsqrtr) / rhsq2 enth(i) = (((qr(i-1,4)+pl)/rhsqrtl & + (ql(i ,4)+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*u(i)**2 - Y(2,i)*q0) a(i) = dsqrt(a2) smax(i) = dmax1(dmax1(dabs(ur-ar-(ul-al)),dabs(ur-ul)), & dabs(ur+ar-(ul+al))) c 10 continue c do 30 i=2-mbc,mx+mbc c c # find a1 thru a3, the coefficients of the 4 eigenvectors: c do k = 1, 4 delta(k) = ql(i,k) - qr(i-1,k) enddo drho = delta(1) + delta(2) c a2 = gamma1/a(i)**2 * (drho*0.5d0*u(i)**2 - delta(2)*q0 & - u(i)*delta(3) + delta(4)) a3 = 0.5d0*( a2 - ( u(i)*drho - delta(3) )/a(i) ) a1 = a2 - a3 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,3,1) = a1*(u(i) - a(i)) wave(i,4,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,3,2) = (drho - a2)*u(i) wave(i,4,2) = (drho - a2)*0.5d0*u(i)**2 + & q0*(delta(2) - Y(2,i)*a2) s(i,2) = u(i) c c # 3-wave wave(i,1,3) = a3*Y(1,i) wave(i,2,3) = a3*Y(2,i) wave(i,3,3) = a3*(u(i) + a(i)) wave(i,4,3) = a3*(enth(i) + u(i)*a(i)) s(i,3) = u(i)+a(i) c c write (6,*) i,s(i,1),s(i,2),s(i,3) c 30 continue c c # compute fluxes as c # F(Ur,Ul) = 0.5*( f(Ur)+f(Ul) - |A|(Ur-Ul) ) c -------------------------- c call flx1(maxmx,meqn,mbc,mx,qr,maux,auxr,apdq) call flx1(maxmx,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 (roe) then do 40 i = 2-mbc, mx+mbc do 40 m=1,meqn amdq(i,m) = 0.5d0*(fl(i,m)+fr(i-1,m)) 40 continue c c tau = 0.1d0 do 50 i = 2-mbc, mx+mbc do 50 m=1,meqn sw = 0.d0 do 60 mw=1,mwaves sl = dabs(s(i,mw)) c # Alternative (worse results for 2nd order) c # for f=1.8, 10 pts/L1/2, 0.25 stable, 0.125 instable c if (efix) sl = sl + 0.5d0*0.d0*smax(i) c c # Also possible: Skip the linearly degenerate field if (efix.and.dabs(s(i,mw)).lt.smax(i).and.mw.ne.2) & sl = s(i,mw)**2/(2.d0*smax(i))+ & 0.5d0*smax(i) c if (efix.and.dabs(s(i,mw)).lt.smax(i)) c & sl = tau*(s(i,mw)**2/(2.d0*smax(i))+ c & 0.5d0*smax(i))+(1.d0-tau)*sl sw = sw + sl*wave(i,m,mw) 60 continue amdq(i,m) = amdq(i,m) - 0.5d0*sw 50 continue endif c if (hll) then do 55 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,3) + wave(i,3,1) El = qr(i-1,4) + wave(i,4,1) pl = gamma1*(El - rho2l*q0 - 0.5d0*rhoul**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,3) - wave(i,3,3) Er = ql(i,4) - wave(i,4,3) pr = gamma1*(Er - rho2r*q0 - 0.5d0*rhour**2/(rho1r+rho2r)) if (rho1r+rho2r.le.0.d0.or.pr.le.0.d0) & hllfix = .true. c if (hllfix) then c if (roe) write (6,*) 'Switching to HLL in',i c rl = qr(i-1,1) + qr(i-1,2) ul = qr(i-1,3)/rl pl = gamma1*(qr(i-1,4) - qr(i-1,2)*q0 - & 0.5d0*qr(i-1,3)**2/rl) al = dsqrt(gamma*pl/rl) c rr = ql(i ,1) + ql(i ,2) ur = ql(i ,3)/rr pr = gamma1*(ql(i ,4) - ql(i ,2)*q0 - & 0.5d0*ql(i ,3)**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) amdq(i,m) = fr(i-1,m) if (sr.le.0.d0) amdq(i,m) = fl(i,m) if (sl.lt.0.d0.and.sr.gt.0.d0) & amdq(i,m) = (sr*fr(i-1,m) - sl*fl(i,m) + & sl*sr*(ql(i,m)-qr(i-1,m)))/ (sr-sl) enddo s(i,1) = sl s(i,2) = 0.d0 s(i,3) = sr endif 55 continue endif c if (pfix) then do 70 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 70 m=1,2 if (amdr.gt.0.d0) then Z = qr(i-1,m)/rhol else Z = ql(i ,m)/rhor endif amdq(i,m) = Z*amdr 70 continue endif c do 80 i = 2-mbc, mx+mbc do 80 m=1,meqn apdq(i,m) = -amdq(i,m) 80 continue c return end c