c c c ========================================================= subroutine rp1eu(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux, & auxl,auxr,wave,s,fl,fr) c ========================================================= c c # solve Riemann problems for the 1D Euler equations using c # the Flux-Vector-Splitting of Vijayasundaram 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, s the speeds, c # 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 routine step1, 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 ql(1-mbc:maxmx+mbc, meqn) dimension qr(1-mbc:maxmx+mbc, meqn) dimension s(1-mbc:maxmx+mbc, mwaves) dimension wave(1-mbc:maxmx+mbc, meqn, mwaves) dimension fl(1-mbc:maxmx+mbc, meqn) dimension fr(1-mbc:maxmx+mbc, meqn) double precision el(3), er(3) common /param/ gamma,gamma1 c c # Method returns fluxes c ------------ common /rpnflx/ mrpnflx mrpnflx = 1 c do 10 i=2-mbc,mx+mbc rhol = qr(i-1,1) rhoul = qr(i-1,2) rhoEl = qr(i-1,3) rhor = ql(i ,1) rhour = ql(i ,2) rhoEr = ql(i ,3) pl = gamma1*(rhoEl - 0.5d0*rhoul**2/rhol) pr = gamma1*(rhoEr - 0.5d0*rhour**2/rhor) c rho = 0.5d0*(rhol + rhor ) rhou = 0.5d0*(rhoul + rhour) rhoE = 0.5d0*(rhoEl + rhoEr) u = rhou/rho p = gamma1*(rhoE - 0.5d0*rhou**2/rho) H = (rhoE+p)/rho if (p.le.0.d0.or.rho.le.0.d0.or.pl.le.0.d0.or.pr.le.0.d0) then write (6,*) 'Error in middle state in',i write (6,*) p,pl,pr,rho,rhol,rhor,rhoul,rhour,rhoEl,rhoEr endif a = dsqrt(gamma*p/rho) f = 0.5d0/a**2 c el(1) = 0.5d0*(u-a + dabs(u-a)) el(2) = 0.5d0*(u + dabs(u) ) el(3) = 0.5d0*(u+a + dabs(u+a)) er(1) = 0.5d0*(u-a - dabs(u-a)) er(2) = 0.5d0*(u - dabs(u) ) er(3) = 0.5d0*(u+a - dabs(u+a)) c zl = el(1)-el(3) zr = er(1)-er(3) ol = el(1)-2.d0*el(2)+el(3) or = er(1)-2.d0*er(2)+er(3) dul = a*(rhol*u-rhoul) dur = a*(rhor*u-rhour) dEl = gamma1*(rhoEl+0.5d0*rhol*u**2-rhoul*u) dEr = gamma1*(rhoEr+0.5d0*rhor*u**2-rhour*u) f1 = f*(zl*dul + ol*dEl + zr*dur + or*dEr) f2 = a*f*(ol*dul + zl*dEl + or*dur + zr*dEr) c fl(i,1) = rhol *el(2) + rhor *er(2) + f1 fl(i,2) = rhoul*el(2) + rhour*er(2) + u*f1 - f2 fl(i,3) = rhoEl*el(2) + rhoEr*er(2) + H*f1 - u*f2 c do 20 m = 1,meqn fr(i,m) = -fl(i,m) 20 continue c do 10 mw=1,mwaves s(i,mw) = dmax1(dabs(el(mw)),dabs(er(mw))) do 10 m=1,meqn wave(i,m,mw) = 0.d0 10 continue c return end c