! ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
MODULE Generic_EigenSystem
SAVE
INTEGER CLESLOG_ROE
INTEGER CLESLOG_ROE_STATES
PARAMETER (CLESLOG_ROE=2)
PARAMETER (CLESLOG_ROE_STATES=1)
INTERFACE EigenSystem
MODULE PROCEDURE EigenSystem_Roe
END INTERFACE
CONTAINS
SUBROUTINE EigenSystem_Roe(ux0, ux1, vx0, vx1, lc, rc, lamda)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
use cles_interfaces
! ----
IMPLICIT NONE
include 'cles.i'
! ---- Eigen System
DOUBLE PRECISION, INTENT(OUT) :: lc(nvars,nvars)
DOUBLE PRECISION, INTENT(OUT) :: lamda(nvars)
DOUBLE PRECISION, INTENT(OUT) :: rc(nvars,nvars)
! ---- Primative variables
DOUBLE PRECISION, INTENT(IN) :: vx0(nvars), vx1(nvars)
DOUBLE PRECISION, INTENT(IN) :: ux0(ncomps), ux1(ncomps)
DOUBLE PRECISION :: vroe(nvars)
DOUBLE PRECISION :: rho, u, v, w, z, u2, v2, w2
! ---- for 3D
! ---- Left (l) and Right (r) states
DOUBLE PRECISION :: c2, gamma, R, T, c, gamma1
DOUBLE PRECISION :: rhol, rhor, rholr
DOUBLE PRECISION :: xi, muz, cg1, qsq2, qsq, s, H, c2g, xil ,xir
double precision :: mu(nscal+1)
character*(nscal*12+10) slog
INTEGER :: i, j, is
! Roe averages
rhol=sqrt(vx0(1))
rhor=sqrt(vx1(1))
rholr = rhol+rhor
! ---- assemble the rho-averaged states
vroe(1) = rhol*rhor
vroe(2:nvars)=(vx0(2:nvars)*rhol+vx1(2:nvars)*rhor)/rholr
H = ((ux0(5)+vx0(5))/rhol+(ux1(5)+vx1(5))/rhor)/rholr
T = (ux0(nvars+1)*rhol+ux1(nvars+1)*rhor)/rholr
rho = vroe(1)
u = vroe(2)
v = vroe(3)
w = vroe(4)
u2 = u*u
v2 = v*v
w2 = w*w
! chemistry parameters
call cles_roe(ux0, ux1, ncomps, vx0, vx1, nvars, gamma, R, mu, nscal)
if ( cleslog_active(CLESLOG_ROE, CLESLOG_ROE_STATES) &
.eq. CLES_TRUE ) then
write(slog, 111) 'Gamma = ', gamma
call cleslog_log(slog//char(0))
write(slog, 111) 'R = ', R
call cleslog_log(slog//char(0))
write(slog, 111) 'Mu = ', (mu(i), i=1,nscal)
call cleslog_log(slog//char(0))
call cleslog_log_flush()
111 format(A8,1(E12.4))
endif
gamma1 = gamma-1.0d0
! Sound speed
c2 = gamma * R * T
qsq = u2 + v2 + w2
qsq2 = ((vx0(2)**2+vx0(3)**2+vx0(4)**2)*rhol &
+(vx1(2)**2+vx1(3)**2+vx1(4)**2)*rhor)/rholr
c2 = c2 + gamma1*(qsq2-qsq)*0.5d0
c = SQRT(c2)
cg1 = c/gamma1
c2g = c2/gamma1
if ( nscal .gt. 0 ) then
muz = SUM(vroe(6:nvars)*mu(1:nscal))
else
muz = 0.0d0
endif
s = -H+v2+w2+c2g+muz
! ---- The right eigenvectors
rc(:,:) = 0.0d0
rc(1,1) = 1.d0
rc(2,1) = u-c
rc(3,1) = v
rc(4,1) = w
rc(5,1) = H-c*u
if ( nscal .gt. 0 ) rc(6:nvars,1) = vroe(6:nvars)
rc(3,2) = 1.0d0
rc(5,2) = v
rc(4,3) = 1.0d0
rc(5,3) = w
rc(1,4) = 1.0d0
rc(2,4) = u
rc(3,4) = v
rc(4,4) = w
rc(5,4) = H-c2g
if ( nscal .gt. 0 ) rc(6:nvars,4) = vroe(6:nvars)
do is=1,nscal
rc(5,4+is) = mu(is)
rc(5+is,4+is) = 1.0d0
enddo
rc(1,nvars) = 1.0d0
rc(2,nvars) = u+c
rc(3,nvars) = v
rc(4,nvars) = w
rc(5,nvars) = H + c * u
if ( nscal .gt. 0 ) rc(6:nvars,nvars) = vroe(6:nvars)
! ---- The left eigenvectors
lc(:,:) = 0.0d0
lc(1,1) = (s+u*(cg1+u))
lc(1,2) = -(cg1+u)
lc(1,3) = -v
lc(1,4) = -w
lc(1,5) = 1.0d0
if ( nscal .gt. 0 ) lc(1,6:nvars) = -mu(1:nscal)
! normalize
lc(1,:) = lc(1,:)/(2.0d0*c2g)
lc(2,1) = -v
lc(2,3) = 1.0d0
lc(3,1) = -w
lc(3,4) = 1.0d0
lc(4,1) = -s-u2+c2g
lc(4,2) = u
lc(4,3) = v
lc(4,4) = w
lc(4,5) = -1.0d0
if ( nscal .gt. 0 ) lc(4,6:nvars) = mu(1:nscal)
! normalize
lc(4,:) = lc(4,:)/c2g
do is=1, nscal
xi = vroe(5+is)
lc(4+is,1) = -xi
lc(4+is,5+is) = 1.0d0
enddo
lc(nvars,1) = (s+u*(u-cg1))
lc(nvars,2) = (cg1-u)
lc(nvars,3) = -v
lc(nvars,4) = -w
lc(nvars,5) = 1.0d0
if ( nscal .gt. 0 ) lc(nvars,6:nvars) = -mu(1:nscal)
! normalize
lc(nvars,:) = lc(nvars,:)/(2.0d0*c2g)
! ---- the eigenvalues ordered from smaller to larger values
lamda(1) = u-c
lamda(2:nvars-1) = u
lamda(nvars) = u+c
END SUBROUTINE EigenSystem_Roe
END MODULE Generic_EigenSystem