vtf-logo

src/3d/equations/euler/rpznd/ip3euzndrfl.f

c
c     Boundary conditions for ghost-fluid methods.
c 
c     Copyright (C) 2003-2007 California Institute of Technology
c     Ralf Deiterding, ralf@amroc.net
c
c     -----------------------------------------------------
c     Internal reflecting physical boundary conditions
c     for Euler equations
c     -----------------------------------------------------
c
c     Transformation of vector of conserved quantities
c     into primitives (rho1,rho2,u,v,w,p)
c
c     =====================================================
      subroutine it3euzndrfl(mx,my,mz,meqn,q,qt)
c     =====================================================
      implicit none
c     
      common /param/  gamma,gamma1,q0
      double precision    gamma,gamma1,q0
      integer   i, j, k, mx, my, mz, meqn
      double precision    q(meqn,mx,my,mz), qt(meqn,mx,my,mz)
c
      do 10 k = 1, mz
         do 10 j = 1, my
            do 10 i = 1, mx
               qt(1,i,j,k) = q(1,i,j,k)
               qt(2,i,j,k) = q(2,i,j,k)
               qt(3,i,j,k) = q(3,i,j,k)/(q(1,i,j,k)+q(2,i,j,k))
               qt(4,i,j,k) = q(4,i,j,k)/(q(1,i,j,k)+q(2,i,j,k))
               qt(5,i,j,k) = q(5,i,j,k)/(q(1,i,j,k)+q(2,i,j,k))
               qt(6,i,j,k) = gamma1*(q(6,i,j,k) - q(2,i,j,k)*q0 - 
     &           0.5d0*(q(3,i,j,k)**2+q(4,i,j,k)**2+q(5,i,j,k)**2)/
     &              (q(1,i,j,k)+q(2,i,j,k)))
 10   continue
c         
      return
      end
c
c     -----------------------------------------------------
c
c     Construction of reflective boundary conditions from
c     mirrored primitive values and application in
c     conservative form in local patch
c
c     =====================================================
      subroutine ip3euzndrfl(q,mx,my,mz,lb,ub,meqn,nc,idx,
     &     qex,xc,phi,vn,maux,auex,dx,time)
c     =====================================================

      implicit none

      common /param/  gamma,gamma1,q0
      double precision    gamma,gamma1,q0
      integer   mx, my, mz, meqn, maux, nc, idx(3,nc), lb(3), 
     &     ub(3)
      double precision    q(meqn,mx,my,mz), qex(meqn,nc), xc(3,nc), 
     &     phi(nc), vn(3,nc), auex(maux,nc), dx(3), time
c
c     Local variables
c
      integer   i, j, k, n, stride, getindx
      double precision    rho1, rho2, u, v, w, p, vl
c
      stride = (ub(1) - lb(1))/(mx-1)
c
      do 100 n = 1, nc

         i = getindx(idx(1,n), lb(1), stride)
         j = getindx(idx(2,n), lb(2), stride)
         k = getindx(idx(3,n), lb(3), stride)
c
         rho1 =  qex(1,n)
         rho2 =  qex(2,n)
         u    = -qex(3,n)       
         v    = -qex(4,n)
         w    = -qex(5,n)
         p    =  qex(6,n)
c
c        # Add boundary velocities if available
         if (maux.ge.3) then
            u = u + auex(1,n)
            v = v + auex(2,n)
            w = w + auex(3,n)
         endif
c
         vl = 2.d0*(u*vn(1,n)+v*vn(2,n)+w*vn(3,n))
         u = qex(3,n) + vl*vn(1,n) 
         v = qex(4,n) + vl*vn(2,n) 
         w = qex(5,n) + vl*vn(3,n) 
c
         q(1,i,j,k) = rho1
         q(2,i,j,k) = rho2
         q(3,i,j,k) = u*(rho1+rho2)
         q(4,i,j,k) = v*(rho1+rho2)
         q(5,i,j,k) = w*(rho1+rho2)
         q(6,i,j,k) = p/gamma1+0.5d0*(rho1+rho2)*(u**2+v**2+w**2)+
     &        rho2*q0
c
 100  continue
c
      return
      end
c
c     -----------------------------------------------------
c
c     Injection of conservative extrapolated values in local patch
c
c     =====================================================
      subroutine ip3euzndex(q,mx,my,mz,lb,ub,meqn,nc,idx,
     &     qex,xc,phi,vn,maux,auex,dx,time)
c     =====================================================
c
      implicit none
c
      common /param/  gamma,gamma1,q0
      double precision    gamma,gamma1,q0
      integer   mx, my, mz, meqn, maux, nc, idx(3,nc), lb(3), 
     &     ub(3)
      double precision    q(meqn,mx,my,mz), qex(meqn,nc), xc(3,nc), 
     &     phi(nc), vn(3,nc), auex(maux,nc), dx(3), time
c
c     Local variables
c
      integer   i, j, k, n, stride, getindx
      double precision    rho1, rho2, u, v, w, p, vl
c
      stride = (ub(1) - lb(1))/(mx-1)
c
      do 100 n = 1, nc

         i = getindx(idx(1,n), lb(1), stride)
         j = getindx(idx(2,n), lb(2), stride)
         k = getindx(idx(3,n), lb(3), stride)
c
         rho1 = qex(1,n)
         rho2 = qex(2,n)
         u    = qex(3,n)       
         v    = qex(4,n)
         w    = qex(5,n)
         p    = qex(6,n)
c
c        # Prescribe normal velocity vector
         vl = u*vn(1,n)+v*vn(2,n)+w*vn(3,n)
         u = vl*vn(1,n) 
         v = vl*vn(2,n) 
         w = vl*vn(3,n) 
c
         q(1,i,j,k) = rho1
         q(2,i,j,k) = rho2
         q(3,i,j,k) = u*(rho1+rho2)
         q(4,i,j,k) = v*(rho1+rho2)
         q(5,i,j,k) = w*(rho1+rho2)
         q(6,i,j,k) = p/gamma1+0.5d0*(rho1+rho2)*(u**2+v**2+w**2)+
     &        rho2*q0
c
 100  continue
c
      return
      end
c

<