c c c ===================================================== subroutine rpn2euznd(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr, & maux,auxl,auxr,wave,s,fl,fr) c ===================================================== c c # solve Riemann problems for the 2D ZND-Euler equations using c # Steger & Warming - Flux Vector Splitting 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 # This data is along a slice in the x-direction if ixy=1 c # or the y-direction if ixy=2. 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 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: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 fl(1-mbc:maxm+mbc, meqn) dimension fr(1-mbc:maxm+mbc, meqn) double precision el(3), er(3) common /param/ gamma,gamma1,q0 c c # Method returns fluxes c ------------ common /rpnflx/ mrpnflx mrpnflx = 1 c c # set mu to point to the component of the system that corresponds c # to momentum in the direction of this slice, mv to the orthogonal c # momentum: c if (ixy.eq.1) then mu = 3 mv = 4 else mu = 4 mv = 3 endif c c # Steger & Warming - Flux Vector Splitting c do 10 i=2-mbc,mx+mbc rhol = qr(i-1,1)+qr(i-1,2) rhor = ql(i ,1)+ql(i ,2) Y1l = qr(i-1,1)/rhol Y2l = qr(i-1,2)/rhol Y1r = ql(i ,1)/rhor Y2r = ql(i ,2)/rhor ul = qr(i-1,mu)/rhol ur = ql(i ,mu)/rhor vl = qr(i-1,mv)/rhol vr = ql(i ,mv)/rhor pl = gamma1*(qr(i-1,5) - qr(i-1,2)*q0 - & 0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2)/rhol) pr = gamma1*(ql(i ,5) - ql(i ,2)*q0 - & 0.5d0*(ql(i ,mu)**2+ql(i ,mv)**2)/rhor) Hl = (qr(i-1,5)+pl)/rhol Hr = (ql(i ,5)+pr)/rhor c al2 = gamma*pl/rhol al = dsqrt(al2) ar2 = gamma*pr/rhor ar = dsqrt(ar2) c el(1) = 0.5d0*(ul-al + dabs(ul-al)) el(2) = 0.5d0*(ul + dabs(ul) ) el(3) = 0.5d0*(ul+al + dabs(ul+al)) er(1) = 0.5d0*(ur-ar - dabs(ur-ar)) er(2) = 0.5d0*(ur - dabs(ur) ) er(3) = 0.5d0*(ur+ar - dabs(ur+ar)) c facl = 0.5d0*rhol/gamma facr = 0.5d0*rhor/gamma c taul = facl*(el(1) + 2.d0*gamma1*el(2) + el(3)) taur = facr*(er(1) + 2.d0*gamma1*er(2) + er(3)) zetal = al*facl*(el(1)-el(3)) zetar = ar*facr*(er(1)-er(3)) c fl(i,1) = Y1l*taul + Y1r*taur fl(i,2) = Y2l*taul + Y2r*taur fl(i,mu) = ul*taul - zetal + ur*taur - zetar fl(i,mv) = vl*taul + vr*taur fl(i,5) = Hl*taul - ul*zetal - 2.d0*el(2)*facl*al2 + & Hr*taur - ur*zetar - 2.d0*er(2)*facr*ar2 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