c c c ===================================================== subroutine rpn3euznd(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr, & maux,auxl,auxr,wave,s,fl,fr) c ===================================================== c c # solve Riemann problems for the 3D ZND-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 # This data is along a slice in the x-direction if ixyz=1 c # the y-direction if ixyz=2. c # the z-direction if ixyz=3. c 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 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: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) dimension auxl(1-mbc:maxm+mbc, maux, 3) dimension auxr(1-mbc:maxm+mbc, maux, 3) double precision el(3), er(3) common /param/ gamma,gamma1,q0 c c # Riemann solver returns flux differences 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 and mw to the c # orthogonal momentums: c if(ixyz .eq. 1)then mu = 3 mv = 4 mw = 5 else if(ixyz .eq. 2)then mu = 4 mv = 5 mw = 3 else mu = 5 mv = 3 mw = 4 endif c do 10 i=2-mbc,mx+mbc rho1l = qr(i-1,1) rho1r = ql(i ,1) rho2l = qr(i-1,2) rho2r = ql(i ,2) rhoul = qr(i-1,mu) rhour = ql(i ,mu) rhovl = qr(i-1,mv) rhovr = ql(i ,mv) rhowl = qr(i-1,mw) rhowr = ql(i ,mw) rhoEl = qr(i-1,6) rhoEr = ql(i ,6) rhol = rho1l+rho2l rhor = rho1r+rho2r c rho = 0.5d0*(rhol + rhor ) rho1 = 0.5d0*(rho1l + rho1r) rho2 = 0.5d0*(rho2l + rho2r) rhou = 0.5d0*(rhoul + rhour) rhov = 0.5d0*(rhovl + rhovr) rhoE = 0.5d0*(rhoEl + rhoEr) c Y1 = rho1/rho Y2 = rho2/rho u = rhou/rho v = rhov/rho w = rhov/rho p = gamma1*(rhoE - rho2*q0 - 0.5d0*rho*(u**2+v**2+w**2)) H = (rhoE+p)/rho if (p.le.0.d0.or.rho.le.0.d0) & write (6,*) 'Error in middle state in',i,p,pl,pr, & rho,rhol,rhor,a,al,ar a = dsqrt(gamma*p/rho) f = 0.5d0/a**2 c el1 = 0.5d0*(u-a + dabs(u-a)) el2 = 0.5d0*(u + dabs(u) ) el3 = 0.5d0*(u+a + dabs(u+a)) er1 = 0.5d0*(u-a - dabs(u-a)) er2 = 0.5d0*(u - dabs(u) ) er3 = 0.5d0*(u+a - dabs(u+a)) c zl = el1-el3 ol = el1-2.d0*el2+el3 zr = er1-er3 or = er1-2.d0*er2+er3 dul = a*(rhol*u-rhoul) dur = a*(rhor*u-rhour) dEl = gamma1*(rhoEl-rho2l*q0+0.5d0*rhol*(u**2+v**2+w**2)- & rhoul*u-rhovl*v-rhowl*w) dEr = gamma1*(rhoEr-rho2r*q0+0.5d0*rhor*(u**2+v**2+w**2)- & rhour*u-rhovr*v-rhowr*w) 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) = rho1l*el2 + rho1r*er2 + Y1*f1 fl(i,2) = rho2l*el2 + rho2r*er2 + Y2*f1 fl(i,mu) = rhoul*el2 + rhour*er2 + u*f1 - f2 fl(i,mv) = rhovl*el2 + rhovr*er2 + v*f1 fl(i,mw) = rhowl*el2 + rhowr*er2 + w*f1 fl(i,6) = rhoEl*el2 + rhoEr*er2 + H*f1 - u*f2 c do 20 m = 1,meqn fr(i,m) = -fl(i,m) 20 continue c s(i,1) = u-a s(i,2) = u s(i,3) = u+a do 10 mws=1,mwaves do 10 m=1,meqn wave(i,m,mws) = 0.d0 10 continue c return end c