c -----------------------------------------------------
c Predefined internal physical boundary conditions
c for Euler equations in WENO solver
c -----------------------------------------------------
c Transformation of vector of conserved quantities
c into primitives (rho,u,v,w,p,s1,s2,dc)
c =====================================================
SUBROUTINE it3eu(mx,my,mz,meqn,q,qt)
c =====================================================
IMPLICIT NONE
INTEGER mx, my, mz, meqn
DOUBLE PRECISION q(meqn,mx,my,mz)
DOUBLE PRECISION qt(meqn,mx,my,mz)
c ---- Local variables
INTEGER i, j, k, m, nvars, useLES, ierr
DOUBLE PRECISION Temperature(1)
call cles_getiparam('nvars', nvars, ierr)
call cles_getiparam('useles', useLES, ierr)
DO k = 1, mz
DO j = 1, my
DO i = 1, mx
! rho
qt(1,i,j,k) = q(1,i,j,k)
! u, v, w
do m=2, nvars
qt(m,i,j,k) = q(m,i,j,k)/q(1,i,j,k)
enddo
! p
call cles_eqstate(q(1,i,j,k),meqn,qt(1,i,j,k),nvars,1,
$ useLES)
! temperature
qt(nvars+1,i,j,k) = q(nvars+1,i,j,k)
! dcflag
qt(nvars+2,i,j,k) = 0.0
! all others
DO m=nvars+3, meqn
qt(m,i,j,k) = q(m,i,j,k)
END Do
END DO
END DO
END DO
RETURN
END
c -----------------------------------------------------
c Construction of reflective boundary conditions from
c mirrored primitive values and application in
c conservative form in local patch in 3D
c -----------------------------------------------------
c =====================================================
SUBROUTINE ip3eurfl(q,mx,my,mz,lb,ub,meqn,nc,idx,
$ qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
IMPLICIT NONE
INTEGER mx, my, mz, meqn, maux, nc, idx(3,nc), lb(3), ub(3)
DOUBLE PRECISION xc(3,nc),
$ phi(nc), vn(3,nc), auex(maux,nc), dx(3), time
DOUBLE PRECISION q(meqn, mx, my, mz)
DOUBLE PRECISION qex(meqn,nc)
c ---- Local variables
INTEGER i, j, k, n, m, stride, getindx, nvars, useViscous, useLES
INTEGER ierr
DOUBLE PRECISION u(3), ul
call cles_getiparam('nvars', nvars, ierr)
call cles_getiparam('useviscous', useViscous, ierr)
call cles_getiparam('useles', useLES, ierr)
stride = (ub(1) - lb(1))/(mx-1)
DO 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)
u(1) = -qex(2,n)
u(2) = -qex(3,n)
u(3) = -qex(4,n)
c ---- Add boundary velocities if available
if (maux.ge.3) then
u(1) = u(1) + auex(1,n)
u(2) = u(2) + auex(2,n)
u(3) = u(3) + auex(3,n)
endif
u(1) = 2.d0*u(1)
u(2) = 2.d0*u(2)
u(3) = 2.d0*u(3)
c ---- Invert entire velocity vector for Navier-Stokes
IF (useViscous.eq.1.and.useLES.eq.0) THEN
qex(2,n) = qex(2,n) + u(1)
qex(3,n) = qex(3,n) + u(2)
qex(4,n) = qex(4,n) + u(3)
c ---- Invert only normal velocity vector for Euler or LES
ELSE
ul = u(1)*vn(1,n)+u(2)*vn(2,n)+u(3)*vn(3,n)
qex(2,n) = qex(2,n) + ul*vn(1,n)
qex(3,n) = qex(3,n) + ul*vn(2,n)
qex(4,n) = qex(4,n) + ul*vn(3,n)
ENDIF
q(1,i,j,k) = qex(1,n)
do m=2, nvars
q(m,i,j,k) = qex(m,n)*qex(1,n)
enddo
call cles_inveqst(q(1,i,j,k),meqn,qex(1,n),nvars,1,useLES)
! temperature
q(nvars+1,i,j,k) = qex(nvars+1,n)
do m=nvars+3, meqn ! skip dcflag
q(m,i,j,k) = qex(m,n)
enddo
END DO
RETURN
END
c -----------------------------------------------------
c Injection of conservative extrapolated values in local patch
c -----------------------------------------------------
c =====================================================
SUBROUTINE ip3euex(q,mx,my,mz,lb,ub,meqn,nc,idx,
$ qex,xc,phi,vn,maux,auex,dx,time)
c =====================================================
IMPLICIT NONE
INTEGER mx, my, mz, meqn, maux, nc, idx(3,nc), lb(3), ub(3)
DOUBLE PRECISION xc(3,nc),
$ phi(nc), vn(3,nc), auex(maux,nc), dx(3), time
DOUBLE PRECISION q(meqn, mx, my, mz)
DOUBLE PRECISION qex(meqn,nc)
c ---- Local variables
INTEGER i, j, k, n, m, stride, getindx, nvars, useLES, ierr
DOUBLE PRECISION u, v, w, vl
call cles_getiparam('nvars', nvars, ierr)
call cles_getiparam('useles', useLES, ierr)
stride = (ub(1) - lb(1))/(mx-1)
DO 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)
u = qex(2,n)
v = qex(3,n)
w = qex(4,n)
c ---- Prescribe normal velocity vector
vl = u*vn(1,n)+v*vn(2,n)+w*vn(3,n)
qex(2,n) = vl*vn(1,n)
qex(3,n) = vl*vn(2,n)
qex(4,n) = vl*vn(3,n)
! rho
q(1,i,j,k) = qex(1,n)
! rho (u,v,w)
do m=2, nvars
q(m,i,j,k) = qex(m,n)*qex(1,n)
enddo
! E and T
call cles_inveqst(q(1,i,j,k),meqn,qex(1,n),nvars,1,useLES)
! temperature
q(nvars+1,i,j,k) = qex(nvars+1,n)
do m=nvars+3, meqn ! skip dcflag
q(m,i,j,k) = qex(m,n)
enddo
END DO
RETURN
END