!-----------------------------------------------------------------------
! 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
!
! Copyright (C) 2002 Ralf Deiterding
! Brandenburgische Universitaet Cottbus
!
!-----------------------------------------------------------------------
subroutine physbd(u,mx,my,lb,ub,lbbnd,ubbnd,shapebnd,
& xc,dx,dir,bnd,mb,time,meqn)
implicit double precision (a-h,o-z)
include "cuser.i"
integer meqn, mx, my, mb, dir
integer lb(2), ub(2), lbbnd(2), ubbnd(2), shapebnd(2)
double precision u(meqn, mx, my), xc(2), dx(2), bnd(mb,2,2), time
! Local variables
integer i, j, imin, imax, jmin, jmax, jsym
integer stride, stridebnd, getindx
integer isx(4), isy(4)
c
data isx /1,-1,1,1/
data isy /1,1,-1,1/
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 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