• VTF
  • FSI
  • AMROC
  • SFC
  • Motion
  • STLIB
  • Main Page
  • src/2d/equations/euler/rp/saux2efix.f

    c
    c     =====================================================
          subroutine setaux(maxmx,maxmy,meqn,mbc,ibx,iby,mx,my,q,
         &     aux,maux,xc,yc,dx,dy,t,dt)
    c     =====================================================
    c
    c     # Multi-dimensional computation of wave speeds, used in 
    c     # multi-dimensional entropy correction in rpn2euefix.
    c
    c     # Copyright (C) 2002 Ralf Deiterding
    c     # Brandenburgische Universitaet Cottbus
    c
          implicit double precision (a-h, o-z)
    c      
          dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc, 
         &     1-iby*mbc:maxmy+iby*mbc)
          dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc, 
         &     1-iby*mbc:maxmy+iby*mbc)
          common /param/  gamma,gamma1
    c
          do 10 j = 1-iby*mbc, my+iby*mbc
             do 10 i = 1-ibx*mbc, mx+ibx*mbc
                p = gamma1*(q(4,i,j) - 0.5d0*(q(2,i,j)**2 + 
         &                  q(3,i,j)**2)/q(1,i,j))
                a = dsqrt(gamma*p/q(1,i,j))
                u = q(2,i,j)/q(1,i,j)
                v = q(3,i,j)/q(1,i,j)
    c     
                if (i.gt.1-ibx*mbc) then
                   pl = gamma1*(q(4,i-1,j) - 0.5d0*(q(2,i-1,j)**2 + 
         &                      q(3,i-1,j)**2)/q(1,i-1,j))
                   al = dsqrt(gamma*pl/q(1,i-1,j))
                   ul = q(2,i-1,j)/q(1,i-1,j)
                   aux(3,i,j) = dmax1(dmax1(dabs(u-a-(ul-al)),dabs(u-ul)),
         &                                  dabs(u+a-(ul+al)))
                else
                   aux(3,i,j) = 0.d0
                endif
    c
                if (j.gt.1-iby*mbc) then
                   pd = gamma1*(q(4,i,j-1) - 0.5d0*(q(2,i,j-1)**2 + 
         &                      q(3,i,j-1)**2)/q(1,i,j-1))
                   ad = dsqrt(gamma*pd/q(1,i,j-1))
                   vd = q(3,i,j-1)/q(1,i,j-1)
                   aux(4,i,j) = dmax1(dmax1(dabs(v-a-(vd-ad)),dabs(v-vd)),
         &                                  dabs(v+a-(vd+ad)))
                else
                   aux(4,i,j) = 0.d0
                endif
     10   continue
    c
          do 11 j = 1-iby*mbc, my+iby*mbc
             aux(1,1-ibx*mbc,j) = 0.d0
             aux(2,1-ibx*mbc,j) = 0.d0
     11   continue
    c
          do 12 i = 1-ibx*mbc, mx+ibx*mbc
             aux(1,i,1-iby*mbc) = 0.d0
             aux(2,i,1-iby*mbc) = 0.d0
     12   continue
    c
          do 20 j = 2-iby*mbc, my+iby*mbc
             do 20 i = 2-ibx*mbc, mx+ibx*mbc
                aux(1,i,j) = dmax1(dmax1(aux(3,i,j),aux(4,i  ,j  )),
         &                                          aux(4,i-1,j  ))
                if (j.lt.my+iby*mbc) 
         &           aux(1,i,j) = dmax1(dmax1(aux(1,i,j),aux(4,i  ,j+1)),
         &                                               aux(4,i-1,j+1))
    c
                aux(2,i,j) = dmax1(dmax1(aux(4,i,j),aux(3,i  ,j  )),
         &                                          aux(3,i  ,j-1))
                if (i.lt.mx+ibx*mbc) 
         &           aux(2,i,j) = dmax1(dmax1(aux(2,i,j),aux(3,i+1,j  )),
         &                                               aux(3,i+1,j-1))
     20   continue
    c
          return
          end
    c
    

<