!----------------------------------------------------------------------- ! Physical boundary conditions for 1d Euler equations. ! Zero order outflow at both sides. ! Interface: ! mx,my := shape of grid function ! ! u(,) := grid function ! ! lb(1) := lower bound for grid ! ub(1) := upper bound for grid ! lbbnd(1) := lower bound for boundary region ! ubnd(1) := upper bound for boundary region ! shapebnd(1) := shape of boundary region ! xc(1) := lower left corner of grid ! dx(1) := grid spacing ! dir := at which side of the grid is the boundary? ! bnd(,2,1) := lower left and upper right corner of global grid and ! of mb-1 internal boundary regions !----------------------------------------------------------------------- subroutine physbd(u,mx,lb,ub,lbbnd,ubbnd,shapebnd, & xc,dx,dir,bnd,mb,time,meqn) implicit none include 'heairacesf6.i' include 'cuser.i' integer meqn, mx, mb, dir integer lb(1), ub(1), lbbnd(1), ubbnd(1), shapebnd(1) double precision u(meqn, mx), xc(1), dx(1), bnd(mb,2,1), time ! Local variables integer i, m, imin, imax, isym integer stride, stridebnd, getindx integer isx(10) double precision gamma1 c data isx /1,-1,1,1,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) if(imax .gt. mx .or. imin .lt. 1) then write(0,*)'INDEX ERROR in physbd' end if if(trunc.eq.1) then c # set up for A. BCs do i=n_bc, 1, -1 if ( time .gt. time_bc(i) ) then u_Sair = u_bc(i) p_Sair = p_bc(i) rho_Sair = rho_bc(i) go to 100 endif enddo 100 continue end if ! Left Side --- Reflecting !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if (dir.eq.0) then if (trunc.eq.1) then c # load for the A BCs gamma1 = gamma_air-1.0d0 do 110 i = imax, imin, -1 u(1,i) = rho_Sair u(2,i) = u(1,i) * u_Sair u(3,i) = 0.0d0 u(4,i) = 0.0d0 u(5,i) = p_Sair/gamma1+rho_Sair*(u_Sair**2)/2.0d0 u(6,i) = 1.d0*u(1,i) u(7,i) = 0.d0 u(8,i) = p_Sair/(rho_Sair*rgas_air) 110 continue return else do 200 i = imax, imin, -1 isym = 2*imax+1-i do 200 m = 1, meqn u(m,i) = isx(m)*u(m,isym) 200 continue return endif end if ! Right Side --- Reflecting !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if (dir.eq.1) then do 210 i = imin, imax isym = 2*imin-1-i do 210 m = 1, meqn u(m,i) = isx(m)*u(m,isym) 210 continue return end if c return end