!----------------------------------------------------------------------- ! 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 ! ! Copyright (C) 2003-2007 California Institute of Technology ! Ralf Deiterding, ralf@cacr.caltech.edu ! !----------------------------------------------------------------------- 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) !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! 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 go to (100,200,300,400) dir+1 ! Left Side --- Inflow/Outflow !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 100 continue do 110 i = imax, imin, -1 if (xc(1)+(i-0.5d0)*dx(1).lt.bnd(1,1,1)) then do 120 j = jmin, jmax do 120 m = 1, meqn u(m,i,j) = u(m,i+1,j) 120 continue endif 110 continue return ! Right Side --- Inflow/Outflow or Reflection !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200 continue do 210 i = imin, imax if (xc(1)+(i-0.5d0)*dx(1).gt.bnd(1,2,1)) then if (ibnd.eq.1) then do 220 j = jmin, jmax do 220 m = 1, meqn u(m,i,j) = u(m,i-1,j) 220 continue else isym = 2*imin-1-i do 230 j = jmin, jmax do 230 m = 1, meqn u(m,i,j) = u(m,isym,j)*isx(m) 230 continue endif endif 210 continue return ! Bottom Side --- Reflection !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 300 continue do 310 j = jmax, jmin, -1 if (xc(2)+(j-0.5d0)*dx(2).lt.bnd(1,1,2)) then jsym = 2*jmax+1-j do 320 i = imin, imax do 320 m = 1, meqn u(m,i,j) = u(m,i,jsym)*isy(m) 320 continue endif 310 continue return ! Top Side --- Reflection !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 400 continue do 410 j = jmin, jmax if (xc(2)+(j-0.5d0)*dx(2).gt.bnd(1,2,2)) then jsym = 2*jmin-1-j do 420 i = imin, imax do 420 m = 1, meqn u(m,i,j) = u(m,i,jsym)*isy(m) 420 continue endif 410 continue return return end