applications/euler/2d/Ramp/src/physbd2.f

!-----------------------------------------------------------------------
! Physical boundary conditions
! Interface:
!   mx,my    := shape of grid function
!
!   u(,) := grid function
!
!   lb(2) := lower bound for grid
!   ub(2) := upper bound for grid
!   lbbnd(2) := lower bound for boundary region
!   ubnd(2) := upper bound for boundary region
!   shapebnd(2) := shape of boundary region 
!   xc(2) := lower left corner of grid
!   dx(2) := grid spacing
!   dir := at which side of the grid is the boundary?
!   bnd(,2,2) := lower left and upper right corner of global grid and 
!      of mb-1 internal boundary regions 
!-----------------------------------------------------------------------

      subroutine physbd(u,mx,my,lb,ub,lbbnd,ubbnd,shapebnd,
     &     xc,dx,dir,bnd,mb,time,args,argc)

      implicit double precision (a-h,o-z)
      include  "cuser.i"

      integer   argc, args, mx, my, mb, dir
      integer   lb(2), ub(2), lbbnd(2), ubbnd(2), shapebnd(2)
      real*8    u(args, mx, my), xc(2), dx(2), bnd(mb,2,2), time

!      Local variables
      integer   i, j, imin, imax, jmin, jmax, jsym
      integer   stride, stridebnd, meqn, getindx
      integer   isx(4), isy(4)
c
      data isx /1,-1,1,1/
      data isy /1,1,-1,1/

      meqn = args

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!      See definition of member-function extents() in BBox.h 
!      for calculation of stride
         
      stride = (ub(1) - lb(1))/(mx-1)
      
      if (meqn .ne. 4) then
         write(0,*)'ERROR in physbd: '
      else

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Find working domain
         
         imin = getindx(max(lbbnd(1), lb(1)), lb(1), stride)
         imax = getindx(min(ubbnd(1), ub(1)), lb(1), stride)
         
         jmin = getindx(max(lbbnd(2), lb(2)), lb(2), stride)
         jmax = getindx(min(ubbnd(2), ub(2)), lb(2), stride)
         
         if(imax .gt. mx .or. jmax .gt. my .or. 
     &        imin .lt. 1 .or. jmin .lt. 1) then
            write(0,*)'INDEX ERROR in physbd'
         end if



!        Left Side --- Inflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         if (dir.eq.0) then
            do 200 i = imax, imin, -1
               do 200 j = jmin, jmax
                  u(1,i,j) = rhoshk
                  u(2,i,j) = rhoshk * ushk
                  u(3,i,j) = rhoshk * vshk
                  u(4,i,j) = rhoshk *(eshk +.5d0*(ushk*ushk+vshk*vshk))
 200        continue
            return
         end if

!        Right Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         if (dir.eq.1) then
            do 210 i = imin, imax
               do 210 j = jmin, jmax
                  do 210 m = 1, meqn
                     u(m,i,j) = u(m,i-1,j)
 210        continue
            return
         end if
         
!        Bottom Side --- reflecting (for ramp part else inflow)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         if (dir.eq.2) then
            do 220 j = jmax, jmin, -1
               jsym = 2*jmax+1-j
               do 220 i = imin, imax
                  xcen = xc(1) + (i-.5d0)*dx(1)
                  if (xcen .le. disp) then
                     u(1,i,j) = rhoshk
                     u(2,i,j) = rhoshk*ushk
                     u(3,i,j) = rhoshk*vshk
                     u(4,i,j) = rhoshk*(eshk+.5d0*(ushk*ushk+vshk*vshk))
                  else
                     do 221 m=1,meqn
                        u(m,i,j)=u(m,i,jsym)*isy(m)
 221                 continue
                  end if
 220        continue
            return
         end if

!        Top Side --- inflow (left of actual shock position) 
!                     outflow (right of actual shock position) 
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         if (dir.eq.3) then
            x0      =  disp + 10.d0*time/cos(pi/6.d0)
            y0      =  bnd(1,1,2)
            do 230 j = jmin, jmax
               do 230 i = imin, imax
                  xcen = xc(1) + (i-.5d0)*dx(1)
                  ycen = xc(2) + (j-.5d0)*dx(2)
                  call cellave(xcen-dx(1)/2.d0,ycen-dx(2)/2.d0,
     $                 dx(1),dx(2),wl)
                  rho = (1.d0-wl)*rhoamb + wl*rhoshk
                  uu = (1.d0-wl)*uamb + wl*ushk
                  v = (1.d0-wl)*vamb + wl*vshk
                  p = (1.d0-wl)*pamb + wl*pshk
                  u(1,i,j) = rho
                  u(2,i,j) = rho*uu
                  u(3,i,j) = rho*v
                  u(4,i,j) = p/gamma1 + .5d0*rho*(uu*uu+v*v)
 230        continue
            return
         end if
         
      end if
      
      return
      end