!-----------------------------------------------------------------------
! Physical boundary conditions for 3d Euler equations.
! Outflow at all sides.
! Physical boundary conditions
! Interface:
! mx,my,mz := shape of grid
!
! u(,) := grid function
!
! lb(3) := lower bound for grid
! ub(3) := upper bound for grid
! lbbnd(3) := lower bound for boundary region
! ubnd(3) := upper bound for boundary region
! shapebnd(3) := shape of boundary region
! xc(3) := lower left corner of grid
! dx(3) := grid spacing
! dir := at which side of the grid is the boundary?
! bnd(,2,3) := lower left and upper right corner of global grid and
! of mb-1 internal boundary regions
!
! Copyright (C) 2002 Ralf Deiterding
! Brandenburgische Universitaet Cottbus
!
! Copyright (C) 2003-2007 California Institute of Technology
! Ralf Deiterding, ralf@cacr.caltech.edu
!
!-----------------------------------------------------------------------
subroutine physbd(u,mx,my,mz,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, mz, mb, dir
integer lb(3), ub(3), lbbnd(3), ubbnd(3), shapebnd(3)
double precision u(meqn,mx,my,mz), xc(3), dx(3), bnd(mb,2,3), time
! Local variables
integer i, j, k, imin, imax, jmin, jmax, kmin, kmax, m
integer stride, getindx, isym, jsym, ksym
integer isx(6), isz(6)
data isx /1,1,-1,1,1,1/
data isz /1,1,1,1,-1,1/
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! See definition of member-function extents() in BBox.h
! for calculation of stride
stride = (ub(1) - lb(1))/(mx-1)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 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)
kmin = getindx(max(lbbnd(3), lb(3)), lb(3), stride)
kmax = getindx(min(ubbnd(3), ub(3)), lb(3), stride)
if (imax .gt. mx .or. jmax .gt. my .or. kmax .gt. mz .or.
& imin .lt. 1 .or. jmin .lt. 1 .or. kmin .lt. 1) then
write(0,*)'INDEX ERROR in physbd'
end if
go to (100,200,300,400,500,600) dir+1
! Left Side --- Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100 continue
do 110 i = imax, imin, -1
if (xc(1)+(i-0.5d0)*dx(1).lt.bnd(1,1,1)) then
isym = 2*imax+1-i
do 120 k = kmin, kmax
do 120 j = jmin, jmax
do 120 m = 1, meqn
u(m,i,j,k) = u(m,isym,j,k)*isx(m)
120 continue
endif
110 continue
return
! Right Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200 continue
do 210 i = imin, imax
if (xc(1)+(i-0.5d0)*dx(1).gt.bnd(1,2,1)) then
do 220 k = kmin, kmax
do 220 j = jmin, jmax
do 220 m = 1, meqn
u(m,i,j,k) = u(m,i-1,j,k)
220 continue
endif
210 continue
return
! Bottom Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300 continue
do 310 j = jmax, jmin, -1
if (xc(2)+(j-0.5d0)*dx(2).lt.bnd(1,1,2)) then
do 320 k = kmin, kmax
do 320 i = imin, imax
do 320 m = 1, meqn
u(m,i,j,k) = u(m,i,j+1,k)
320 continue
endif
310 continue
return
! Top Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400 continue
do 410 j = jmin, jmax
if (xc(2)+(j-0.5d0)*dx(2).gt.bnd(1,2,2)) then
do 420 k = kmin, kmax
do 420 i = imin, imax
do 420 m = 1, meqn
u(m,i,j,k) = u(m,i,j-1,k)
420 continue
endif
410 continue
return
! Front Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
500 continue
do 510 k = kmax, kmin, -1
if (xc(3)+(k-0.5d0)*dx(3).lt.bnd(1,1,3)) then
do 520 j = jmin, jmax
do 520 i = imin, imax
do 520 m = 1, meqn
u(m,i,j,k) = u(m,i,j,k+1)
520 continue
endif
510 continue
return
! Back Side --- Reflection
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
600 continue
if (Nback.eq.1) then
do 610 k = kmin, kmax
if (xc(3)+(k-0.5d0)*dx(3).gt.bnd(1,2,3)) then
ksym = 2*kmin-1-k
do 620 j = jmin, jmax
do 620 i = imin, imax
xcen = xc(1) + (i-.5d0)*dx(1)
ycen = xc(2) + (j-.5d0)*dx(2)
r = dsqrt(xcen*xcen+ycen*ycen)
if (r.lt.rfback) then
do m = 1, meqn
u(m,i,j,k) = u(m,i,j,ksym)*isz(m)
enddo
else
do m = 1, meqn
u(m,i,j,k) = u(m,i,j,k-1)
enddo
endif
620 continue
endif
610 continue
else
do 630 k = kmin, kmax
if (xc(3)+(k-0.5d0)*dx(3).gt.bnd(1,2,3)) then
do 640 j = jmin, jmax
do 640 i = imin, imax
do 640 m = 1, meqn
u(m,i,j,k) = u(m,i,j,k-1)
640 continue
endif
630 continue
endif
return
end