c c ========================================================= subroutine rp1eurhok(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux, & auxl,auxr,wave,s,amdq,apdq) c ========================================================= c c # solve Riemann problems for the thermally perfect 1D multi-component c # Euler equations using Roe's 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 # 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 routine step1, rp is called with ql = qr = q. c c # Copyright (C) 2002 Ralf Deiterding, Georg Bader 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 amdq(1-mbc:maxmx+mbc, meqn) dimension apdq(1-mbc:maxmx+mbc, meqn) c c # local storage c --------------- parameter (max2 = 10002) !# assumes at most 10000 grid points with mbc=2 dimension u(-1:max2),enth(-1:max2),a(-1:max2) dimension g1a2(-1:max2) logical efix, pfix c c define local arrays c include "ck.i" c dimension delta(LeNsp+2) dimension rkl(LeNsp), rkr(LeNsp) dimension hkl(LeNsp), hkr(LeNsp) dimension Y(LeNsp,-1:max2), pk(LeNsp,-1:max2) c data efix /.true./ !# use entropy fix for transonic rarefactions data pfix /.true./ !# use Larrouturou's positivity fix for species c c # Riemann solver returns fluxes c ------------ common /rpnflx/ mrpnflx mrpnflx = 0 c mu = Nsp+1 mE = Nsp+2 mT = Nsp+3 c c # Compute Roe-averaged quantities: c do 20 i=2-mbc,mx+mbc rhol = 0.d0 rhor = 0.d0 do k = 1, Nsp rkl(k) = qr(i-1,k) rkr(k) = ql(i ,k) rhol = rhol + rkl(k) rhor = rhor + rkr(k) enddo if( rhol.le.1.d-10 ) then write(6,*) 'negative total density, left', rhol stop endif if( rhor.le.1.d-10 ) then write(6,*) 'negative total density, right', rhor stop endif c c # compute left/right rho/W and rho*Cp c rhoWl = 0.d0 rhoWr = 0.d0 do k = 1, Nsp rhoWl = rhoWl + rkl(k)/Wk(k) rhoWr = rhoWr + rkr(k)/Wk(k) enddo c c # left/right Temperatures already calculated c Tl = qr(i-1,mT) Tr = ql(i ,mT) pl = rhoWl*RU*Tl pr = rhoWr*RU*Tr c c # compute quantities for rho-average c rhsqrtl = dsqrt(rhol) rhsqrtr = dsqrt(rhor) rhsq2 = rhsqrtl + rhsqrtr c c # find rho-averaged specific velocity and enthalpy c u(i) = (qr(i-1,mu)/rhsqrtl + & ql(i ,mu)/rhsqrtr) / rhsq2 enth(i) = (((qr(i-1,mE)+pl)/rhsqrtl & + (ql(i ,mE)+pr)/rhsqrtr)) / rhsq2 c c # compute rho-averages for T, cp, and W c T = (Tl * rhsqrtl + Tr * rhsqrtr) / rhsq2 W = rhsq2 / (rhoWl/rhsqrtl + rhoWr/rhsqrtr) c c # evaluate left/right entropies and mean cp c call tabintp( Tl, hkl, hms, Nsp ) call tabintp( Tr, hkr, hms, Nsp ) do k = 1, Nsp Y(k,i) = (rkl(k)/rhsqrtl + rkr(k)/rhsqrtr) / rhsq2 enddo Cp = Cpmix( Tl, Tr, hkl, hkr, Y(1,i) ) gamma1 = RU / ( W*Cp - RU ) gamma = gamma1 + 1.d0 c c # find rho-averaged specific enthalpies, c # compute rho-averaged mass fractions and c # compute partial pressure derivatives c tmp = gamma * RU * T / gamma1 * ht = 0.d0 do k = 1, Nsp hk = (hkl(k)*rhsqrtl + hkr(k)*rhsqrtr) / rhsq2 * ht = ht + Y(k,i)*(hkl(k)*rhsqrtl + hkr(k)*rhsqrtr) / rhsq2 pk(k,i) = 0.5d0*u(i)**2 - hk + tmp / Wk(k) enddo c * write (6,4) qr(i-1,mE)+pl, ql(i,mE)+pr, * & ht+0.5d0*u(i)**2, enth(i), ht+0.5d0*u(i)**2-enth(i) * 4 format(e16.8,e16.8,e16.8,e16.8,e24.16) c c # compute speed of sound c a2 = enth(i)-u(i)**2 do k = 1, Nsp a2 = a2 + Y(k,i) * pk(k,i) enddo g1a2(i) = 1.d0 / a2 a(i) = dsqrt(gamma1*a2) c 20 continue c c do 30 i=2-mbc,mx+mbc c c # find a1 thru a3, the coefficients of the mE eigenvectors: c dpdr = 0.d0 dpY = 0.d0 drho = 0.d0 do k = 1, Nsp delta(k) = ql(i,k) - qr(i-1,k) drho = drho + delta(k) dpdr = dpdr + pk(k,i) * delta(k) dpY = dpY + pk(k,i) * Y(k,i) enddo delta(mu) = ql(i,mu) - qr(i-1,mu) delta(mE) = ql(i,mE) - qr(i-1,mE) c a2 = g1a2(i)*(dpdr - u(i)*delta(mu) + delta(mE)) a3 = 0.5d0*( a2 - ( u(i)*drho - delta(mu) )/a(i) ) a1 = a2 - a3 c c c # Compute the waves. c # Note that the 1+k-waves, for 1 .le. k .le. Nsp travel at c # the same speed and are lumped together in wave(.,.,2). c # The 3-wave is then stored in wave(.,.,3). c do k = 1, Nsp c # 1-wave wave(i,k,1) = a1*Y(k,i) c # 2-wave wave(i,k,2) = delta(k) - Y(k,i)*a2 c # 3-wave wave(i,k,3) = a3*Y(k,i) enddo c # 1-wave wave(i,mu,1) = a1*(u(i)-a(i)) wave(i,mE,1) = a1*(enth(i) - u(i)*a(i)) wave(i,mT,1) = 0.d0 s(i,1) = u(i)-a(i) c c # 2-wave wave(i,mu,2) = u(i)*(drho - a2) wave(i,mE,2) = u(i)**2*(drho - a2) - dpdr + dpY*a2 wave(i,mT,2) = 0.d0 s(i,2) = u(i) c c # 3-wave wave(i,mu,3) = a3*(u(i)+a(i)) wave(i,mE,3) = a3*(enth(i)+u(i)*a(i)) wave(i,mT,3) = 0.d0 s(i,3) = u(i)+a(i) c 30 continue c c # compute Godunov flux f0: 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 mw=1,mwaves if (s(i,mw) .lt. 0.d0) then amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw) else apdq(i,m) = apdq(i,m) + s(i,mw)*wave(i,m,mw) 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 rho = 0.d0 rhoW = 0.d0 do k = 1, Nsp rkl(k) = qr(i-1,k) rho = rho + rkl(k) rhoW = rhoW + rkl(k)/Wk(k) enddo rhou = qr(i-1,mu) rhoE = qr(i-1,mE) T = qr(i-1,mT) rhoCp = avgtabip( T, rkl, cpk, Nsp ) gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0 p = rhoW*RU*T 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 rho = 0.d0 rhoW = 0.d0 do k = 1, Nsp rkl(k) = rkl(k) + wave(i,k,1) rho = rho + rkl(k) rhoW = rhoW + rkl(k)/Wk(k) enddo rhou = rhou + wave(i,mu,1) rhoE = rhoE + wave(i,mE,1) rhoe = rhoE - 0.5d0*rhou**2/rho if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then write(6,*) 'Temperature out of range', T write(6,*) 'state vector 1 before' write(6,*) (rkl(k),k=1,Nsp) endif call SolveTrhok( T, rhoe, rhoW, rkl, Nsp, ifail) rhoCp = avgtabip( T, rkl, cpk, Nsp ) if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then write(6,*) 'Temperature out of range', T write(6,*) 'state vector 1 after' write(6,*) (rkl(k),k=1,Nsp) endif gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0 p = rhoW*RU*T 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 rho = 0.d0 rhoW = 0.d0 do k = 1, Nsp rkr(k) = ql(i,k) rho = rho + rkr(k) rhoW = rhoW + rkr(k)/Wk(k) enddo rhou = ql(i,mu) rhoE = ql(i,mE) T = ql(i,mT) rhoCp = avgtabip( T, rkr, cpk, Nsp ) gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0 p = rhoW*RU*T 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 rho = 0.d0 rhoW = 0.d0 do k = 1, Nsp rkr(k) = rkr(k) - wave(i,k,3) rho = rho + rkr(k) rhoW = rhoW + rkr(k)/Wk(k) enddo rhou = rhou - wave(i,mu,3) rhoE = rhoE - wave(i,mE,3) rhoe = rhoE - 0.5d0*rhou**2/rho if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then write(6,*) 'Temperature out of range', T write(6,*) 'state vector 1 before' write(6,*) (rkr(k),k=1,Nsp) endif call SolveTrhok( T, rhoe, rhoW, rkr, Nsp, ifail) rhoCp = avgtabip( T, rkr, cpk, Nsp ) if ( TabS.gt.T*TABFAC .or. T*TABFAC.gt.TabE ) then write(6,*) 'Temperature out of range', T write(6,*) 'state vector 1 after' write(6,*) (rkr(k),k=1,Nsp) endif gamma = RU / ( rhoCp/rhoW - RU ) + 1.d0 p = rhoW*RU*T 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 c # compute the rightgoing flux differences: c # df = SUM s*wave is the total flux difference and apdq = df - amdq c do 220 m=1,meqn do 220 i = 2-mbc, mx+mbc df = 0.d0 do 210 mw=1,mwaves df = df + s(i,mw)*wave(i,m,mw) 210 continue apdq(i,m) = df - amdq(i,m) 220 continue c 900 continue c if (pfix) then do 70 i=2-mbc,mx+mbc amdr = 0.d0 apdr = 0.d0 rhol = 0.d0 rhor = 0.d0 do k = 1, Nsp amdr = amdr + amdq(i,k) apdr = apdr + apdq(i,k) rhol = rhol + qr(i-1,k) rhor = rhor + ql(i ,k) enddo do 70 k=1,Nsp if (qr(i-1,mu)+amdr.gt.0.d0) then Z = qr(i-1,k)/rhol else Z = ql(i ,k)/rhor endif amdq(i,k) = Z*amdr + (Z-qr(i-1,k)/rhol)*qr(i-1,mu) apdq(i,k) = Z*apdr - (Z-ql(i ,k)/rhor)*ql(i ,mu) 70 continue endif c return end c c c *********************************************************** c double precision function Cpmix( Tl, Tr, hl, hr, Y ) implicit double precision(a-h,o-z) include "ck.i" c dimension Y(*) dimension hl(*), hr(*) data Tol /1.d-6/ c if( dabs(Tr-Tl).gt.Tol ) then Cp = 0.d0 do k = 1, Nsp Cp = Cp + (hr(k)-hl(k)) * Y(k) enddo Cp = Cp / (Tr-Tl) else T = 0.5d0*(Tr+Tl) Cp = avgtabip( T, Y, cpk, Nsp ) endif Cpmix = Cp c return end