vtf-logo

src/2d/equations/euler/rpznd/rpn2euzndefixg.f

c
c     =========================================================
      subroutine rpn2euznd(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
     &                     auxl,auxr,wave,s,amdq,apdq)
c     =========================================================
c
c     # solve Riemann problems for the 2D ZND-Euler equations using Roe's 
c     # approximate Riemann solver.  
c     # Scheme is blended with HLL for robustness and uses a multi-dimensional 
c     # entropy correction to prevent the carbuncle phenomenon.
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, 
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 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 auxl(1-mbc:maxm+mbc, maux)
      dimension auxr(1-mbc:maxm+mbc, maux)
      dimension apdq(1-mbc:maxm+mbc, meqn)
      dimension amdq(1-mbc:maxm+mbc, meqn)
c
c     local arrays -- common block comroe is passed to rpt2eu
c     ------------
      parameter (maxm2 = 10005)  !# assumes at most 10000x10000 grid with mbc=5
      parameter (minm2 = -4)     !# assumes at most mbc=5
      dimension delta(5), fl(minm2:maxm2,5), fr(minm2:maxm2,5)
      logical efix, pfix, hll, roe, hllfix
      common /param/  gamma,gamma1,q0
      common /comroe/ u2v2(minm2:maxm2),u(minm2:maxm2),v(minm2:maxm2),
     &       enth(minm2:maxm2),a(minm2:maxm2),Y(2,minm2:maxm2)
c
      data efix /.true./   !# use entropy fix
      data pfix /.true./   !# use Larrouturou's positivity fix for species
      data hll  /.true./   !# use HLL instead of Roe solver, if unphysical values occur
      data roe  /.true./   !# turn off Roe solver when debugging HLL
c
c     # Riemann solver returns fluxes
c     ------------
      common /rpnflx/ mrpnflx
      mrpnflx = 1
c
      if (minm2.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then
         write(6,*) 'need to increase maxm2 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 to the orthogonal
c     # momentum:
c
      if (ixy.eq.1) then
         mu = 3
         mv = 4
      else
         mu = 4
         mv = 3
      endif
c
c     # note that notation for u and v reflects assumption that the 
c     # Riemann problems are in the x-direction with u in the normal
c     # direciton and v in the orthogonal direcion, but with the above
c     # definitions of mu and mv the routine also works with ixy=2
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 rpt2eu to do the transverse wave splitting.
c
      do 10 i=2-mbc,mx+mbc
c
         pl = gamma1*(qr(i-1,5) - qr(i-1,2)*q0 - 
     &        0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2)/(qr(i-1,1)+qr(i-1,2)))
         pr = gamma1*(ql(i,  5) - ql(i,  2)*q0 - 
     &        0.5d0*(ql(i,  mu)**2+ql(i,  mv)**2)/(ql(i,  1)+ql(i,  2)))
         rhsqrtl = dsqrt(qr(i-1,1) + qr(i-1,2))  
         rhsqrtr = dsqrt(ql(i,  1) + ql(i,  2))
         rhsq2 = rhsqrtl + rhsqrtr
         u(i) = (qr(i-1,mu)/rhsqrtl + ql(i,mu)/rhsqrtr) / rhsq2
         v(i) = (qr(i-1,mv)/rhsqrtl + ql(i,mv)/rhsqrtr) / rhsq2
         u2v2(i) = u(i)**2 + v(i)**2
         enth(i) = (((qr(i-1,5)+pl)/rhsqrtl
     &             + (ql(i  ,5)+pr)/rhsqrtr)) / rhsq2
         Y(1,i) = (qr(i-1,1)/rhsqrtl + ql(i,1)/rhsqrtr) / rhsq2
         Y(2,i) = (qr(i-1,2)/rhsqrtl + ql(i,2)/rhsqrtr) / rhsq2
c        # speed of sound
         a2 = gamma1*(enth(i) - 0.5d0*u2v2(i) - Y(2,i)*q0)
         a(i) = dsqrt(a2) 
c
   10 continue
c
      do 30 i=2-mbc,mx+mbc
c
c        # find a1 thru a3, the coefficients of the 4 eigenvectors:
c
         do k = 1, 5
            delta(k) = ql(i,k) - qr(i-1,k)
         enddo
         drho = delta(1) + delta(2)
c
         a2 = gamma1/a(i)**2 * (drho*0.5d0*u2v2(i) - delta(2)*q0 
     &      - (u(i)*delta(mu)+v(i)*delta(mv)) + delta(5))
         a3 = delta(mv) - v(i)*drho
         a4 = 0.5d0*( a2 - ( u(i)*drho - delta(mu) )/a(i) )
         a1 = a2 - a4 
c
c        # Compute the waves.
c
c        # 1-wave
         wave(i,1,1)  = a1*Y(1,i)
         wave(i,2,1)  = a1*Y(2,i)
         wave(i,mu,1) = a1*(u(i) - a(i))
         wave(i,mv,1) = a1*v(i)
         wave(i,5,1)  = a1*(enth(i) - u(i)*a(i))
         s(i,1) = u(i)-a(i)
c
c        # 2-wave
         wave(i,1,2)  = delta(1) - Y(1,i)*a2
         wave(i,2,2)  = delta(2) - Y(2,i)*a2         
         wave(i,mu,2) = (drho - a2)*u(i)
         wave(i,mv,2) = (drho - a2)*v(i) + a3
         wave(i,5,2)  = (drho - a2)*0.5d0*u2v2(i) + 
     &                  q0*(delta(2) - Y(2,i)*a2)  + a3*v(i)
         s(i,2) = u(i)
c
c        # 3-wave
         wave(i,1,3)  = a4*Y(1,i)
         wave(i,2,3)  = a4*Y(2,i)
         wave(i,mu,3) = a4*(u(i) + a(i))
         wave(i,mv,3) = a4*v(i)
         wave(i,5,3)  = a4*(enth(i) + u(i)*a(i))
         s(i,3) = u(i)+a(i)
c                  
   30 continue
c
c     # compute fluxes as
c     # F(Ur,Ul) = 0.5*( f(Ur)+f(Ul) - |A|(Ur-Ul) )
c     --------------------------
c
      call flx2(ixy,maxm,meqn,mbc,mx,qr,maux,auxr,apdq)
      call flx2(ixy,maxm,meqn,mbc,mx,ql,maux,auxl,amdq)
c
      do 35 i = 1-mbc, mx+mbc
         do 35 m=1,meqn
            fl(i,m) = amdq(i,m)
            fr(i,m) = apdq(i,m)
 35   continue      
c
      if (roe) then
         do 40 i = 2-mbc, mx+mbc
            do 40 m=1,meqn
               amdq(i,m) = 0.5d0*(fl(i,m)+fr(i-1,m))
 40      continue
c     
         do 50 i = 2-mbc, mx+mbc
            do 50 m=1,meqn
               sw = 0.d0
               do 60 mw=1,mwaves
                  sl = dabs(s(i,mw))
                  if (efix.and.dabs(s(i,mw)).lt.auxl(i,ixy))
     &               sl = s(i,mw)**2/(2.d0*auxl(i,ixy))+
     &                    0.5d0*auxl(i,ixy)
                  sw = sw + sl*wave(i,m,mw)
 60            continue
               amdq(i,m) = amdq(i,m) - 0.5d0*sw
 50      continue
      endif
c
      if (pfix) then
         do 70 i=2-mbc,mx+mbc
            amdr = amdq(i,1)+amdq(i,2)
            rhol = qr(i-1,1)+qr(i-1,2)
            rhor = ql(i  ,1)+ql(i  ,2)
            do 70 m=1,2
               if (amdr.gt.0.d0) then
                  Z = qr(i-1,m)/rhol
               else
                  Z = ql(i  ,m)/rhor
               endif
               amdq(i,m) = Z*amdr
 70      continue    
      endif
c
      if (hll) then
         do 55 i = 2-mbc, mx+mbc
            hllfix = .false.
            if (.not.roe) hllfix = .true.
c     
            rho1l = qr(i-1,1)  + wave(i,1,1)
            rho2l = qr(i-1,2)  + wave(i,2,1)
            rhoul = qr(i-1,mu) + wave(i,mu,1)
            rhovl = qr(i-1,mv) + wave(i,mv,1)
            El    = qr(i-1,5)  + wave(i,5,1)
            pl = gamma1*(El - rho2l*q0 - 
     &           0.5d0*(rhoul**2 + rhovl**2)/(rho1l+rho2l))
            if (rho1l+rho2l.le.0.d0.or.pl.le.0.d0) 
     &         hllfix = .true.
c     
            rho1r = ql(i,1)  - wave(i,1,3)
            rho2r = ql(i,2)  - wave(i,2,3)
            rhour = ql(i,mu) - wave(i,mu,3)
            rhovr = ql(i,mv) - wave(i,mv,3)
            Er    = ql(i,5 ) - wave(i,5,3)
            pr = gamma1*(Er - rho2r*q0 - 
     &           0.5d0*(rhour**2 + rhovr**2)/(rho1r+rho2r))
            if (rho1r+rho2r.le.0.d0.or.pr.le.0.d0) 
     &         hllfix = .true.
c     
            if (hllfix) then
c               if (roe) write (6,*) 'Switching to HLL in',i
c     
               rl = qr(i-1,1) + qr(i-1,2)
               ul = qr(i-1,mu)/rl
               pl = gamma1*(qr(i-1,5) - qr(i-1,2)*q0 - 
     &              0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2)/rl)
               al = dsqrt(gamma*pl/rl)
c     
               rr = ql(i  ,1) + ql(i  ,2)
               ur = ql(i  ,mu)/rr
               pr = gamma1*(ql(i  ,5) - ql(i  ,2)*q0 - 
     &              0.5d0*(ql(i  ,mu)**2+ql(i  ,mv)**2)/rr)
               ar = dsqrt(gamma*pr/rr)
c     
               sl = dmin1(ul-al,ur-ar)
               sr = dmax1(ul+al,ur+ar)
c
               do m=1,meqn
                  if (sl.ge.0.d0) amdq(i,m) = fr(i-1,m)
                  if (sr.le.0.d0) amdq(i,m) = fl(i,m)
                  if (sl.lt.0.d0.and.sr.gt.0.d0) 
     &               amdq(i,m) = (sr*fr(i-1,m) - sl*fl(i,m) + 
     &                            sl*sr*(ql(i,m)-qr(i-1,m)))/ (sr-sl)
               enddo
               s(i,1) = sl
               s(i,2) = 0.d0
               s(i,3) = sr
            endif     
 55      continue
      endif
c
      do 80 i = 2-mbc, mx+mbc
         do 80 m=1,meqn
            apdq(i,m) = -amdq(i,m)
 80   continue
c
      return
      end
c

<