c
c Boundary conditions for ghost-fluid methods.
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
c -----------------------------------------------------
c Internal reflecting physical boundary conditions
c for Euler equations for multiple thermally perfect species
c -----------------------------------------------------
c
c Transformation of vector of conserved quantities
c into (rho1,...,rhoK,u,v,p,T,...)
c
c =====================================================
subroutine it2eurhokrfl(mx,my,meqn,q,qt)
c =====================================================
c
implicit double precision(a-h,o-z)
include "ck.i"
c
integer i, j, k, mx, my, meqn
double precision q(meqn,mx,my), qt(meqn,mx,my)
c
do 10 j = 1, my
do 10 i = 1, mx
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rho = rho + q(k,i,j)
rhoW = rhoW + q(k,i,j)/Wk(k)
qt(k,i,j) = q(k,i,j)
enddo
qt(Nsp+1,i,j) = q(Nsp+1,i,j)/rho
qt(Nsp+2,i,j) = q(Nsp+2,i,j)/rho
qt(Nsp+3,i,j) = rhoW*RU*q(Nsp+4,i,j)
do k = Nsp+4, meqn
qt(k,i,j) = q(k,i,j)
enddo
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 ip2eurhokrfl(q,mx,my,lb,ub,meqn,nc,idx,
& qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
c
implicit double precision(a-h,o-z)
include "ck.i"
c
integer mx, my, meqn, maux, nc, idx(2,nc), lb(2), ub(2)
double precision q(meqn, mx, my), qex(meqn,nc), xc(2,nc),
& phi(nc), vn(2,nc), auex(maux,nc), dx(2), time
c
c Local variables
c
integer i, j, k, n, stride, getindx
double precision rho, rhoW, u, v, vl, p, T
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)
c
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rho = rho + qex(k,n)
rhoW = rhoW + qex(k,n)/Wk(k)
q(k,i,j) = qex(k,n)
enddo
c
u = -qex(Nsp+1,n)
v = -qex(Nsp+2,n)
p = qex(Nsp+3,n)
T = p/(rhoW*RU)
c
c # Add boundary velocities if available
if (maux.ge.2) then
u = u + auex(1,n)
v = v + auex(2,n)
endif
c
c # Construct normal velocity vector
c # Tangential velocity remains unchanged
vl = 2.d0*(u*vn(1,n)+v*vn(2,n))
u = qex(Nsp+1,n) + vl*vn(1,n)
v = qex(Nsp+2,n) + vl*vn(2,n)
c
q(Nsp+1,i,j) = u*rho
q(Nsp+2,i,j) = v*rho
q(Nsp+3,i,j) = rho*0.5d0*(u**2+v**2) +
& avgtabip(T,q(1,i,j),hms,Nsp) - p
q(Nsp+4,i,j) = T
c
do k = Nsp+5, meqn
q(k,i,j) = qex(k,n)
enddo
c
100 continue
c
return
end
c
c
c -----------------------------------------------------
c
c Injection of conservative extrapolated values in local patch
c
c =====================================================
subroutine ip2eurhokex(q,mx,my,lb,ub,meqn,nc,idx,
& qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
c
implicit double precision(a-h,o-z)
include "ck.i"
c
integer mx, my, meqn, maux, nc, idx(2,nc), lb(2), ub(2)
double precision q(meqn, mx, my), qex(meqn,nc), xc(2,nc),
& phi(nc), vn(2,nc), auex(maux,nc), dx(2), time
c
c Local variables
c
integer i, j, n, stride, getindx
double precision rho, rhoW, u, v, vl, p, T
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)
c
rho = 0.d0
rhoW = 0.d0
do k = 1, Nsp
rho = rho + qex(k,n)
rhoW = rhoW + qex(k,n)/Wk(k)
q(k,i,j) = qex(k,n)
enddo
c
u = qex(Nsp+1,n)
v = qex(Nsp+2,n)
p = qex(Nsp+3,n)
T = p/(rhoW*RU)
c
c # Prescribe normal velocity vector
vl = u*vn(1,n)+v*vn(2,n)
u = vl*vn(1,n)
v = vl*vn(2,n)
c
q(Nsp+1,i,j) = u*rho
q(Nsp+2,i,j) = v*rho
q(Nsp+3,i,j) = rho*0.5d0*(u**2+v**2) +
& avgtabip(T,q(1,i,j),hms,Nsp) - p
q(Nsp+4,i,j) = T
c
do k = Nsp+5, meqn
q(k,i,j) = qex(k,n)
enddo
c
100 continue
c
return
end
c