c
c Boundary conditions for ghost-fluid methods.
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@cacr.caltech.edu
c
c -----------------------------------------------------
c Internal reflecting physical boundary conditions
c for Euler equations
c -----------------------------------------------------
c
c Transformation of vector of conserved quantities
c into primitives (rho1,rho2,u,v,w,p)
c
c =====================================================
subroutine it3euzndrfl(mx,my,mz,meqn,q,qt)
c =====================================================
implicit none
c
common /param/ gamma,gamma1,q0
double precision gamma,gamma1,q0
integer i, j, k, mx, my, mz, meqn
double precision q(meqn,mx,my,mz), qt(meqn,mx,my,mz)
c
do 10 k = 1, mz
do 10 j = 1, my
do 10 i = 1, mx
qt(1,i,j,k) = q(1,i,j,k)
qt(2,i,j,k) = q(2,i,j,k)
qt(3,i,j,k) = q(3,i,j,k)/(q(1,i,j,k)+q(2,i,j,k))
qt(4,i,j,k) = q(4,i,j,k)/(q(1,i,j,k)+q(2,i,j,k))
qt(5,i,j,k) = q(5,i,j,k)/(q(1,i,j,k)+q(2,i,j,k))
qt(6,i,j,k) = gamma1*(q(6,i,j,k) - q(2,i,j,k)*q0 -
& 0.5d0*(q(3,i,j,k)**2+q(4,i,j,k)**2+q(5,i,j,k)**2)/
& (q(1,i,j,k)+q(2,i,j,k)))
10 continue
c
return
end
c
c -----------------------------------------------------
c
c Construction of reflective boundary conditions from
c mirrored primitive values and application in
c conservative form in local patch
c
c =====================================================
subroutine ip3euzndrfl(q,mx,my,mz,lb,ub,meqn,nc,idx,
& qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
implicit none
common /param/ gamma,gamma1,q0
double precision gamma,gamma1,q0
integer mx, my, mz, meqn, maux, nc, idx(3,nc), lb(3),
& ub(3)
double precision q(meqn,mx,my,mz), qex(meqn,nc), xc(3,nc),
& phi(nc), vn(3,nc), auex(maux,nc), dx(3), time
c
c Local variables
c
integer i, j, k, n, stride, getindx
double precision rho1, rho2, u, v, w, p, vl
c
stride = (ub(1) - lb(1))/(mx-1)
c
do 100 n = 1, nc
i = getindx(idx(1,n), lb(1), stride)
j = getindx(idx(2,n), lb(2), stride)
k = getindx(idx(3,n), lb(3), stride)
c
rho1 = qex(1,n)
rho2 = qex(2,n)
u = -qex(3,n)
v = -qex(4,n)
w = -qex(5,n)
p = qex(6,n)
c
c # Add boundary velocities if available
if (maux.ge.3) then
u = u + auex(1,n)
v = v + auex(2,n)
w = w + auex(3,n)
endif
c
vl = 2.d0*(u*vn(1,n)+v*vn(2,n)+w*vn(3,n))
u = qex(3,n) + vl*vn(1,n)
v = qex(4,n) + vl*vn(2,n)
w = qex(5,n) + vl*vn(3,n)
c
q(1,i,j,k) = rho1
q(2,i,j,k) = rho2
q(3,i,j,k) = u*(rho1+rho2)
q(4,i,j,k) = v*(rho1+rho2)
q(5,i,j,k) = w*(rho1+rho2)
q(6,i,j,k) = p/gamma1+0.5d0*(rho1+rho2)*(u**2+v**2+w**2)+
& rho2*q0
c
100 continue
c
return
end
c
c -----------------------------------------------------
c
c Injection of conservative extrapolated values in local patch
c
c =====================================================
subroutine ip3euzndex(q,mx,my,mz,lb,ub,meqn,nc,idx,
& qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
c
implicit none
c
common /param/ gamma,gamma1,q0
double precision gamma,gamma1,q0
integer mx, my, mz, meqn, maux, nc, idx(3,nc), lb(3),
& ub(3)
double precision q(meqn,mx,my,mz), qex(meqn,nc), xc(3,nc),
& phi(nc), vn(3,nc), auex(maux,nc), dx(3), time
c
c Local variables
c
integer i, j, k, n, stride, getindx
double precision rho1, rho2, u, v, w, p, vl
c
stride = (ub(1) - lb(1))/(mx-1)
c
do 100 n = 1, nc
i = getindx(idx(1,n), lb(1), stride)
j = getindx(idx(2,n), lb(2), stride)
k = getindx(idx(3,n), lb(3), stride)
c
rho1 = qex(1,n)
rho2 = qex(2,n)
u = qex(3,n)
v = qex(4,n)
w = qex(5,n)
p = qex(6,n)
c
c # Prescribe normal velocity vector
vl = u*vn(1,n)+v*vn(2,n)+w*vn(3,n)
u = vl*vn(1,n)
v = vl*vn(2,n)
w = vl*vn(3,n)
c
q(1,i,j,k) = rho1
q(2,i,j,k) = rho2
q(3,i,j,k) = u*(rho1+rho2)
q(4,i,j,k) = v*(rho1+rho2)
q(5,i,j,k) = w*(rho1+rho2)
q(6,i,j,k) = p/gamma1+0.5d0*(rho1+rho2)*(u**2+v**2+w**2)+
& rho2*q0
c
100 continue
c
return
end
c