vtf-logo

src/3d/equations/euler/rpm/rpn3meuhllc.f

c
c
c     ==================================================================
      subroutine rpn3meu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
     &                   maux,auxl,auxr,wave,s,amdq,apdq)
c     ==================================================================
c
c     # Solve Riemann problems for the 3D two-component Euler equations 
c     # using HLLC. Use flux difference splitting formulation for full
c     # compatibility to Wave Propagation Method.
c     
c     # solve Riemann problems along one slice of data.
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 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     # amdq and apdq 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) 2003-2007 California Institute of Technology
c     # Ralf Deiterding, ralf@amroc.net
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 amdq(1-mbc:maxm+mbc, meqn)
      dimension apdq(1-mbc:maxm+mbc, meqn)
      dimension auxl(1-mbc:maxm+mbc, maux, 3)
      dimension auxr(1-mbc:maxm+mbc, maux, 3)
c
c     local arrays -- common block comroe is passed to rpt3eum
c                     
c     ------------
      parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5
      parameter (minmrp = -4)   !# assumes at most mbc=5
      dimension qls(5), qrs(5)
c
      common /comroe/ u2v2w2(minmrp:maxmrp),
     &     u(minmrp:maxmrp),v(minmrp:maxmrp),w(minmrp:maxmrp),
     &     enth(minmrp:maxmrp),a(minmrp:maxmrp),g1a2(minmrp:maxmrp),
     &     euv(minmrp:maxmrp),p(minmrp:maxmrp) 
c
c     # Riemann solver returns flux differences
c     ------------
      common /rpnflx/ mrpnflx
      mrpnflx = 0
c
      if (minmrp.gt.1-mbc .or. maxmrp .lt. maxm+mbc) then
         write(6,*) 'need to increase maxmrp in rpA'
         stop
      endif
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 = 2
         mv = 3
         mw = 4
      else if(ixyz .eq. 2)then
         mu = 3
         mv = 4
         mw = 2
      else
         mu = 4
         mv = 2
         mw = 3
      endif
c
c     # note that notation for u,v, and w reflects assumption that the 
c     # Riemann problems are in the x-direction with u in the normal
c     # direction and v and w in the orthogonal directions, but with the 
c     # above definitions of mu, mv, and mw the routine also works with 
c     # ixyz=2 and ixyz=3
c     # and returns, for example, f0 as the Godunov flux g0 for the
c     # Riemann problems u_t + g(u)_y = 0 in the y-direction.
c
c
c     # Compute the Roe-averaged variables needed in the Roe solver.
c     # These are stored in the common block comroe since they are
c     # later used in routine rpt3eu to do the transverse wave 
c     # splitting.
c
      do 10 i = 2-mbc, mx+mbc
         if (qr(i-1,1).le.0.d0.or.ql(i,1).le.0.d0) then 
            write (6,*) 'Unrecoverable error in density',i
            stop
         endif         
c
         rl = qr(i-1,1)
         ul = qr(i-1,mu)/rl
         vl = qr(i-1,mv)/rl
         wl = qr(i-1,mw)/rl
         El = qr(i-1,5)
         pl = (El-0.5d0*(ul**2+vl**2+wl**2)*rl-qr(i-1,7))/qr(i-1,6)
c     
         rr = ql(i  ,1)
         ur = ql(i  ,mu)/rr
         vr = ql(i  ,mv)/rr
         wr = ql(i  ,mw)/rr
         Er = ql(i  ,5)
         pr = (Er-0.5d0*(ur**2+vr**2+wr**2)*rr-ql(i  ,7))/ql(i  ,6)
c
         rhsqrtl = dsqrt(qr(i-1,1))
         rhsqrtr = dsqrt(ql(i,1))
         rhsq2 = rhsqrtl + rhsqrtr
         gamma1 = rhsq2 / ( qr(i-1,6)*rhsqrtl + ql(i,6)*rhsqrtr ) 
         xjota = ( pl*qr(i-1,6)*rhsqrtl + pr*ql(i,6)*rhsqrtr ) / rhsq2
         p(i) = xjota*gamma1
c
         u(i) = (qr(i-1,mu)/rhsqrtl + ql(i,mu)/rhsqrtr) / rhsq2
         v(i) = (qr(i-1,mv)/rhsqrtl + ql(i,mv)/rhsqrtr) / rhsq2
         w(i) = (qr(i-1,mw)/rhsqrtl + ql(i,mw)/rhsqrtr) / rhsq2
         enth(i) = (((qr(i-1,5)+pl)/rhsqrtl 
     &             + (ql(i,5)+pr)/rhsqrtr)) / rhsq2
         u2v2w2(i) = u(i)**2 + v(i)**2 + w(i)**2
c
         a2 = gamma1*(enth(i) - .5d0*u2v2w2(i))
         if (a2.le.0.d0) then 
            write (6,*) 'Unrecoverable error in speed of sound in',i
            stop
         endif         
         a(i) = dsqrt(a2)
         g1a2(i) = gamma1 / a2
         euv(i) = enth(i) - u2v2w2(i) 
c     
         sl = u(i)-a(i)
         sr = u(i)+a(i)
         ss = (pr-pl+rl*ul*(sl-ul)-rr*ur*(sr-ur))/
     &        (rl*(sl-ul)-rr*(sr-ur))
c
         qrs(1)  = rr*(sr-ur)/(sr-ss)
         qrs(mu) = qrs(1)*ss
         qrs(mv) = qrs(1)*vr
         qrs(mw) = qrs(1)*wr
         qrs(5)  = qrs(1)*(Er/rr+
     &        (ss-ur)*(ss+pr/(rr*(sr-ur))))
c
         qls(1)  = rl*(sl-ul)/(sl-ss)
         qls(mu) = qls(1)*ss
         qls(mv) = qls(1)*vl
         qls(mw) = qls(1)*wl
         qls(5)  = qls(1)*(El/rl+
     &        (ss-ul)*(ss+pl/(rl*(sl-ul))))
c
         do m=1,5
            wave(i,m,1) = qls(m) - qr(i-1,m)
            wave(i,m,2) = qrs(m) - qls(m)
            wave(i,m,3) = ql(i,m) - qrs(m)
         enddo
         do m=6,7
            wave(i,m,1) = 0.d0
            wave(i,m,2) = ql(i,m) - qr(i-1,m)
            wave(i,m,3) = 0.d0
         enddo
c
         s(i,1) = sl
         s(i,2) = ss
         s(i,3) = sr
c
         do m=1,meqn
            amdq(i,m) = 0.d0
            apdq(i,m) = 0.d0
            do mws=1,mwaves
               if (s(i,mws) .lt. 0.d0) then
                  amdq(i,m) = amdq(i,m) + s(i,mws)*wave(i,m,mws)
               else
                  apdq(i,m) = apdq(i,m) + s(i,mws)*wave(i,m,mws)
               endif
            enddo
         enddo
 10   continue
      return
      end

<