c c ========================================================= subroutine rp1euznd(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux, & auxl,auxr,wave,s,dfl,dfr) c ========================================================= c c # Riemann solver for the 1D ZND-Euler equations. c # The waves are computed using the Roe approximation. 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 clawpack 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 dfr(1-mbc:maxmx+mbc, meqn) dimension dfl(1-mbc:maxmx+mbc, meqn) c 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) common /param/ gamma,gamma1,q0 c c define local arrays c dimension delta(4), Y(2,-1:max2) dimension f0(-1:max2,4), fl(-1:max2,4), fr(-1:max2,4) dimension sl(2), sr(2) c c # Riemann solver returns flux differences c ------------ common /rpnflx/ mrpnflx mrpnflx = 0 c if (-1.gt.1-mbc .or. max2 .lt. maxmx+mbc) then write(6,*) 'need to increase max2 in rp' stop endif c c c # Compute Roe-averaged quantities: c do 10 i=2-mbc,mx+mbc c pl = gamma1*(qr(i-1,4) - qr(i-1,2)*q0 - & 0.5d0*qr(i-1,3)**2/(qr(i-1,1)+qr(i-1,2))) pr = gamma1*(ql(i, 4) - ql(i, 2)*q0 - & 0.5d0*ql(i, 3)**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,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) 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 30 continue c c # compute Godunov flux f0: c -------------------------- c c # compute Godunov flux f0 at each interface. c # Uses exact Riemann solver c do 200 i = 2-mbc, mx+mbc c rhol = qr(i-1,1) + qr(i-1,2) rhor = ql(i ,1) + qr(i ,2) Y2l = qr(i-1,2)/rhol Y2r = ql(i ,2)/rhor 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) c c # iterate to find pstar, ustar: c alpha = 1. pstar = 0.5*(pl+pr) wr = dsqrt(pr*rhor) * phi(pstar/pr) wl = dsqrt(pl*rhol) * phi(pstar/pl) c if (pl.eq.pr .and. rhol.eq.rhor) go to 60 c 40 do 50 iter=1,20 p1 = (ul-ur+pr/wr+pl/wl) / (1./wr + 1./wl) pstar = dmax1(p1,1d-6)*alpha + (1.-alpha)*pstar wr1 = wr wl1 = wl wr = dsqrt(pr*rhor) * phi(pstar/pr) wl = dsqrt(pl*rhol) * phi(pstar/pl) if (dmax1(abs(wr1-wr),dabs(wl1-wl)) .lt. 1d-6) & go to 60 50 continue c c # nonconvergence: alpha = alpha/2. if (alpha .gt. 0.001) go to 40 write(6,*) 'no convergence',wr1,wr,wl1,wl wr = .5*(wr+wr1) wl = .5*(wl+wl1) c 60 continue ustar = (pl-pr+wr*ur+wl*ul) / (wr+wl) c c # left wave: c ============ c if (pstar .gt. pl) then c c # shock: sl(1) = ul - wl/rhol sr(1) = sl(1) rho1 = wl/(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 # right wave: c ============= c if (pstar .ge. pr) then c c # shock sl(2) = ur + wr/rhor sr(2) = sl(2) rho2 = wr/(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 # 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 ps = pl Y2s = Y2l else if (sr(1).le.0. .and. ustar.ge. 0.) then rhos = rho1 us = ustar ps = pstar Y2s = Y2l else if (ustar.lt.0. .and. sl(2).ge. 0.) then rhos = rho2 us = ustar ps = pstar Y2s = Y2r else if (sr(2).lt.0) then rhos = rhor us = ur ps = pr Y2s = Y2r 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 Y2s = Y2l 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 Y2s = Y2r endif c f0(i,1) = (1.d0-Y2s)*rhos*us f0(i,2) = Y2s*rhos*us f0(i,3) = rhos*us**2 + ps f0(i,4) = us*(gamma*ps/gamma1 + Y2s*rhos*q0 + 0.5*rhos*us**2) 200 continue c c c # compute fluxes in each cell: c call flx1(maxmx,meqn,mbc,mx,qr,maux,auxr,dfr) call flx1(maxmx,meqn,mbc,mx,ql,maux,auxl,dfl) c do 210 m=1,meqn do 210 i = 1-mbc, mx+mbc fr(i,m) = dfr(i,m) fl(i,m) = dfl(i,m) 210 continue c c # compute the leftgoing and rightgoing flux differences: do 220 m=1,meqn do 220 i = 2-mbc, mx+mbc dfl(i,m) = f0(i,m) - fr(i-1,m) dfr(i,m) = fl(i,m) - f0(i,m) 220 continue c return end c c double precision function phi(w) implicit double precision (a-h,o-z) common /param/ gamma,gamma1,q0 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