!********************************************************************************
!* A collection of three-dimensional Euler numerical fluxes, Version 3 (2012),
!*
!* written by Dr. Katate Masatsuka (info[at]cfdbooks.com),
!*
!* the author of useful CFD books, "I do like CFD" (http://www.cfdbooks.com).
!*
!*------------------------
!* List of Flux Functions:
!*
!* inviscid_roe : Roe flux (using face normal and tangent vectors)
!* inviscid_roe_n : Roe flux (using face normal vector only)
!* inviscid_rotated_rhll : Rotated-RHLL flux (one of the robust fluxes for carbuncle)
!*
!*------------------------
!*
!* These F90 routines were written and made available for an educational purpose.
!* Detailed descripstion of each numerical flux can be found in the original paper,
!* or in popular textbooks, or in the (will-be-available) second volume of
!* "I do like CFD".
!*
!* Note: These subroutines have been prepared for an educational purpose.
!* It is not at all efficient. Think about how you can optimize it.
!* One way to make it efficient is to reduce the number of local variables,
!* by re-using temporary variables as many times as possible.
!*
!* Note: Please let me know if you find bugs. I'll greatly appreciate it and
!* fix the bugs.
!*
!* Version 01
!* 05-06-11: Minor bugs fixed.
!*
!* Version 02
!* 04-06-12: Major bugs fixed for Rotated-RHLL flux.
!* - qL and qR were not computed with the correct normals
!* - Tangent vector (2nd direction) is now computed more carefully.
!* (Bugs reported by Dr. Yoshitaka Nakashima. Thank you!)
!*
!* Version 03
!* 04-16-12: Equal sign has been added in the if-statement for tangent vector to
!* avoid a possible failure (reported by Dr. Yoshitaka Nakashima. Thank you!)
!* Roe flux without tangent vectors (inviscid_roe_n) added.
!* Many more comments added.
!*
!* This file may be updated in future. (to add more flux functions.)
!********************************************************************************
!********************************************************************************
!* -- 3D Roe's Flux Function with an entropy fix --
!*
!* This subroutine computes the Roe flux for the Euler equations
!* in the direction, njk=[nx,ny,nz].
!*
!* P. L. Roe, Approximate Riemann Solvers, Parameter Vectors and Difference
!* Schemes, Journal of Computational Physics, 43, pp. 357-372.
!*
!* Conservative form of the Euler equations:
!*
!* dU/dt + dF/dx + dG/dy + dH/dz = 0
!*
!* This subroutine computes the numerical flux for the flux in the direction,
!* njk=[nx,ny,nz]:
!*
!* Fn = F*nx + G*ny + H*nz = | rho*qn |
!* | rho*qn*u + p*nx |
!* | rho*qn*v + p*ny |
!* | rho*qn*w + p*nz |
!* | rho*qn*H | (qn = u*nx + v*ny + w*nz)
!*
!* The Roe flux is implemented in the following form:
!*
!* Numerical flux = 1/2 [ Fn(UR) + Fn(UL) - |An|dU ],
!*
!* where
!*
!* An = dFn/dU, |An| = R|Lambda|L, dU = UR - UL.
!*
!* The dissipation term, |An|dU, is actually computed as
!*
!* sum_{k=1,5} |lambda_k| * (LdU)_k * r_k,
!*
!* where lambda_k is the k-th eigenvalue, (LdU)_k is the k-th wave strength,
!* and r_k is the k-th right-eigenvector.
!*
!* ------------------------------------------------------------------------------
!* Input: primL(1:5) = Left state (rhoL, uL, vL, wR, pL)
!* primR(1:5) = Right state (rhoR, uR, vR, wR, pR)
!* njk(1:3) = unit face normal vector (nx, ny, nz), pointing from Left to Right.
!*
!* njk
!* Face normal ^ o Right data point
!* | .
!* | .
!* |.
!* -------x-------- Face
!* . Left and right states are
!* . 1. Values at data points for 1st-order accuracy
!* . 2. Extrapolated values at the face midpoint 'x'
!* o Left data point for 2nd/higher-order accuracy.
!*
!*
!* Output: flux(1:5) = The Roe flux with an entropy fix
!* ------------------------------------------------------------------------------
!*
!* Note: This subroutine has been prepared for an educational purpose.
!* It is not at all efficient. Think about how you can optimize it.
!* One way to make it efficient is to reduce the number of local variables,
!* by re-using temporary variables as many times as possible.
!*
!* Note: Please let me know if you find bugs. I'll greatly appreciate it and
!* fix the bugs.
!*
!* Katate Masatsuka, March 2011. http://www.cfdbooks.com
!********************************************************************************
subroutine inviscid_roe(primL, primR, njk, num_flux)
implicit none
integer , parameter :: p2 = selected_real_kind(15) ! Double precision
!Input
real(p2), intent( in) :: primL(5), primR(5) ! Input: primitive variables
real(p2), intent( in) :: njk(3) ! Input: face normal vector
!Output
real(p2), intent(out) :: num_flux(5) ! Output: numerical flux
!Some constants
real(p2) :: zero = 0.0_p2
real(p2) :: one = 1.0_p2
real(p2) :: two = 2.0_p2
real(p2) :: half = 0.5_p2
real(p2) :: fifth = 0.2_p2
!Local variables
real(p2) :: nx, ny, nz ! Normal vector
real(p2) :: mx, my, mz ! Orthogonal tangent vector
real(p2) :: lx, ly, lz ! Another orthogonal tangent vector
real(p2) :: abs_n_cross_l ! Magnitude of n x l
real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components.
real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables.
real(p2) :: qnL, qnR, qmL, qmR, qlL, qlR ! Normal and tangent velocities
real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy
real(p2) :: RT,rho,u,v,w,H,a,qn, ql, qm ! Roe-averages
real(p2) :: drho,dqn,dql,dqm,dp,LdU(5) ! Wave strengths
real(p2) :: ws(5), R(5,5) ! Wave speeds and right-eigenvectors
real(p2) :: dws(5) ! Width of a parabolic fit for entropy fix
real(p2) :: fL(5), fR(5), diss(5) ! Fluxes ad dissipation term
real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats
real(p2) :: temp, tempx, tempy, tempz ! Temoprary variables
! Face normal vector (unit vector)
nx = njk(1)
ny = njk(2)
nz = njk(3)
!Define face tangent vectors, m and l, such that the three vectors, n, m, and l
!are mutually orthogonal.
! l = (lx.ly,lz)
!
! We first carefully choose one tangent vector: l = (lx,ly,lz).
! E.g., if n = (1,0,0), then, (0,-nz,ny) = (0,0,0), which is a problem.
! To avoid such a situation, we choose the 2D tangential vector
! having the maximum length.
tempx = ny*ny + nz*nz
tempy = nz*nz + nx*nx
tempz = nx*nx + ny*ny
if ( tempx >= tempy .and. tempx >= tempz ) then
lx = zero
ly = -nz
lz = ny
elseif ( tempy >= tempx .and. tempy >= tempz ) then
lx = -nz
ly = zero
lz = nx
elseif ( tempz >= tempx .and. tempz >= tempy ) then
lx = -ny
ly = nx
lz = zero
else
! Impossible to happen
write(*,*) "subroutine inviscid_roe: Impossible to happen. Please report the problem."
stop
endif
! Make it the unit vector.
temp = sqrt( lx*lx + ly*ly + lz*lz )
lx = lx/temp
ly = ly/temp
lz = lz/temp
! m = (mx,my,mz)
!
! The other one, m = (mx,my,mz), is chosen as a vector orthogonal to both n and l
! defined by the vector product: m = n x l / |n x l|
mx = ny*lz - nz*ly
my = nz*lx - nx*lz
mz = nx*ly - ny*lx
abs_n_cross_l = sqrt(mx**2 + my**2 + mz**2)
mx = mx / abs_n_cross_l
my = my / abs_n_cross_l
mz = mz / abs_n_cross_l
!(Do you like such ambiguous tangent vectors? Actually, the Roe flux can
! be implemented without any tangent vector. See "I do like CFD, VOL.1",
! or the subroutine "inviscid_roe_n" or "inviscid_rotated_rhll" for details.
! In fact, the resulting flux is independent of the choice of these tangent vectors
! as it should be.)
!Primitive and other variables.
! Left state
rhoL = primL(1)
uL = primL(2)
vL = primL(3)
wL = primL(4)
qnL = uL*nx + vL*ny + wL*nz
qlL = uL*lx + vL*ly + wL*lz
qmL = uL*mx + vL*my + wL*mz
pL = primL(5)
aL = sqrt(gamma*pL/rhoL)
HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
! Right state
rhoR = primR(1)
uR = primR(2)
vR = primR(3)
wR = primR(4)
qnR = uR*nx + vR*ny + wR*nz
qlR = uR*lx + vR*ly + wR*lz
qmR = uR*mx + vR*my + wR*mz
pR = primR(5)
aR = sqrt(gamma*pR/rhoR)
HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)
!First compute the Roe-averaged quantities
! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for
! the Roe-averaged density.
RT = sqrt(rhoR/rhoL)
rho = RT*rhoL !Roe-averaged density
u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity
v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity
w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity
H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy
a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound
qn = u*nx + v*ny + w*nz !Roe-averaged face-normal velocity
ql = u*lx + v*ly + w*lz !Roe-averaged face-tangent velocity
qm = u*mx + v*my + w*mz !Roe-averaged face-tangent velocity
!Wave Strengths
drho = rhoR - rhoL !Density difference
dp = pR - pL !Pressure difference
dqn = qnR - qnL !Normal velocity difference
dql = qlR - qlL !Tangent velocity difference in l
dqm = qmR - qmL !Tangent velocity difference in m
LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength
LdU(2) = drho - dp/(a*a) !Entropy wave strength
LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength
LdU(4) = rho*dql !Shear wave strength
LdU(5) = rho*dqm !Shear wave strength
!Absolute values of the wave speeds
ws(1) = abs(qn-a) !Left-moving acoustic wave speed
ws(2) = abs(qn) !Entropy wave speed
ws(3) = abs(qn+a) !Right-moving acoustic wave speed
ws(4) = abs(qn) !Shear wave speed
ws(5) = abs(qn) !Shear wave speed
!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.
!NOTE: It avoids vanishing wave speeds by making a parabolic fit near ws = 0.
dws(1) = fifth
if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
dws(3) = fifth
if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )
!Right Eigenvectors
! Left-moving acoustic wave
R(1,1) = one
R(2,1) = u - a*nx
R(3,1) = v - a*ny
R(4,1) = w - a*nz
R(5,1) = H - a*qn
! Entropy wave
R(1,2) = one
R(2,2) = u
R(3,2) = v
R(4,2) = w
R(5,2) = half*(u*u + v*v + w*w)
! Right-moving acoustic wave
R(1,3) = one
R(2,3) = u + a*nx
R(3,3) = v + a*ny
R(4,3) = w + a*nz
R(5,3) = H + a*qn
! Shear wave
R(1,4) = zero
R(2,4) = lx
R(3,4) = ly
R(4,4) = lz
R(5,4) = ql
! Shear wave
R(1,5) = zero
R(2,5) = mx
R(3,5) = my
R(4,5) = mz
R(5,5) = qm
!Dissipation Term: |An|(UR-UL) = R|Lambda|L*dU = sum_k of [ ws(k) * R(:,k) * L*dU(k) ]
diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) + ws(3)*LdU(3)*R(:,3) &
+ ws(4)*LdU(4)*R(:,4) + ws(5)*LdU(5)*R(:,5)
!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)
fL(1) = rhoL*qnL
fL(2) = rhoL*qnL * uL + pL*nx
fL(3) = rhoL*qnL * vL + pL*ny
fL(4) = rhoL*qnL * wL + pL*nz
fL(5) = rhoL*qnL * HL
fR(1) = rhoR*qnR
fR(2) = rhoR*qnR * uR + pR*nx
fR(3) = rhoR*qnR * vR + pR*ny
fR(4) = rhoR*qnR * wR + pR*nz
fR(5) = rhoR*qnR * HR
! This is the numerical flux: Roe flux = 1/2 *[ Fn(UL)+Fn(UR) - |An|(UR-UL) ]
num_flux = half * (fL + fR - diss)
!Normal max wave speed in the normal direction.
! wsn = abs(qn) + a
end subroutine inviscid_roe
!--------------------------------------------------------------------------------
!********************************************************************************
!* -- 3D Roe's Flux Function with an entropy fix and without tangent vectors --
!*
!* NOTE: This version does not use any tangent vector.
!* See "I do like CFD, VOL.1" about how tangent vectors are eliminated.
!*
!* This subroutine computes the Roe flux for the Euler equations
!* in the direction, njk=[nx,ny,nz].
!*
!* P. L. Roe, Approximate Riemann Solvers, Parameter Vectors and Difference
!* Schemes, Journal of Computational Physics, 43, pp. 357-372.
!*
!* Conservative form of the Euler equations:
!*
!* dU/dt + dF/dx + dG/dy + dH/dz = 0
!*
!* This subroutine computes the numerical flux for the flux in the direction,
!* njk=[nx,ny,nz]:
!*
!* Fn = F*nx + G*ny + H*nz = | rho*qn |
!* | rho*qn*u + p*nx |
!* | rho*qn*v + p*ny |
!* | rho*qn*w + p*nz |
!* | rho*qn*H | (qn = u*nx + v*ny + w*nz)
!*
!* The Roe flux is implemented in the following form:
!*
!* Numerical flux = 1/2 [ Fn(UR) + Fn(UL) - |An|dU ],
!*
!* where
!*
!* An = dFn/dU, |An| = R|Lambda|L, dU = UR - UL.
!*
!* The dissipation term, |An|dU, is actually computed as
!*
!* sum_{k=1,5} |lambda_k| * (LdU)_k * r_k,
!*
!* where lambda_k is the k-th eigenvalue, (LdU)_k is the k-th wave strength,
!* and r_k is the k-th right-eigenvector.
!*
!* ------------------------------------------------------------------------------
!* Input: primL(1:5) = Left state (rhoL, uL, vL, wR, pL)
!* primR(1:5) = Right state (rhoR, uR, vR, wR, pR)
!* njk(1:3) = unit face normal vector (nx, ny, nz), pointing from Left to Right.
!*
!* njk
!* Face normal ^ o Right data point
!* | .
!* | .
!* |.
!* -------x-------- Face
!* . Left and right states are
!* . 1. Values at data points for 1st-order accuracy
!* . 2. Extrapolated values at the face midpoint 'x'
!* o Left data point for 2nd/higher-order accuracy.
!*
!*
!* Output: flux(1:5) = The Roe flux with an entropy fix
!* ------------------------------------------------------------------------------
!*
!* Note: This subroutine has been prepared for an educational purpose.
!* It is not at all efficient. Think about how you can optimize it.
!* One way to make it efficient is to reduce the number of local variables,
!* by re-using temporary variables as many times as possible.
!*
!* Note: Please let me know if you find bugs. I'll greatly appreciate it and
!* fix the bugs.
!*
!* Katate Masatsuka, April 2012. http://www.cfdbooks.com
!********************************************************************************
subroutine inviscid_roe_n(primL, primR, njk, num_flux)
implicit none
integer , parameter :: p2 = selected_real_kind(15) ! Double precision
!Input
real(p2), intent( in) :: primL(5), primR(5) ! Input: primitive variables
real(p2), intent( in) :: njk(3) ! Input: face normal vector
!Output
real(p2), intent(out) :: num_flux(5) ! Output: numerical flux
!Some constants
real(p2) :: zero = 0.0_p2
real(p2) :: one = 1.0_p2
real(p2) :: two = 2.0_p2
real(p2) :: half = 0.5_p2
real(p2) :: fifth = 0.2_p2
!Local variables
real(p2) :: nx, ny, nz ! Normal vector
real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components.
real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables.
real(p2) :: qnL, qnR ! Normal velocities
real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy
real(p2) :: RT,rho,u,v,w,H,a,qn ! Roe-averages
real(p2) :: drho,dqn,dp,LdU(4) ! Wave strengths
real(p2) :: du, dv, dw ! Velocity differences
real(p2) :: ws(4), R(5,4) ! Wave speeds and right-eigenvectors
real(p2) :: dws(4) ! Width of a parabolic fit for entropy fix
real(p2) :: fL(5), fR(5), diss(5) ! Fluxes ad dissipation term
real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats
! Face normal vector (unit vector)
nx = njk(1)
ny = njk(2)
nz = njk(3)
!Primitive and other variables.
! Left state
rhoL = primL(1)
uL = primL(2)
vL = primL(3)
wL = primL(4)
qnL = uL*nx + vL*ny + wL*nz
pL = primL(5)
aL = sqrt(gamma*pL/rhoL)
HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
! Right state
rhoR = primR(1)
uR = primR(2)
vR = primR(3)
wR = primR(4)
qnR = uR*nx + vR*ny + wR*nz
pR = primR(5)
aR = sqrt(gamma*pR/rhoR)
HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)
!First compute the Roe-averaged quantities
! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for
! the Roe-averaged density.
RT = sqrt(rhoR/rhoL)
rho = RT*rhoL !Roe-averaged density
u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity
v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity
w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity
H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy
a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound
qn = u*nx + v*ny + w*nz !Roe-averaged face-normal velocity
!Wave Strengths
drho = rhoR - rhoL !Density difference
dp = pR - pL !Pressure difference
dqn = qnR - qnL !Normal velocity difference
LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength
LdU(2) = drho - dp/(a*a) !Entropy wave strength
LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength
LdU(4) = rho !Shear wave strength (not really, just a factor)
!Absolute values of the wave Speeds
ws(1) = abs(qn-a) !Left-moving acoustic wave
ws(2) = abs(qn) !Entropy wave
ws(3) = abs(qn+a) !Right-moving acoustic wave
ws(4) = abs(qn) !Shear waves
!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.
!NOTE: It avoids vanishing wave speeds by making a parabolic fit near ws = 0.
dws(1) = fifth
if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
dws(3) = fifth
if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )
!Right Eigenvectors
!Note: Two shear wave components are combined into one, so that tangent vectors
! are not required. And that's why there are only 4 vectors here.
! See "I do like CFD, VOL.1" about how tangent vectors are eliminated.
! Left-moving acoustic wave
R(1,1) = one
R(2,1) = u - a*nx
R(3,1) = v - a*ny
R(4,1) = w - a*nz
R(5,1) = H - a*qn
! Entropy wave
R(1,2) = one
R(2,2) = u
R(3,2) = v
R(4,2) = w
R(5,2) = half*(u*u + v*v + w*w)
! Right-moving acoustic wave
R(1,3) = one
R(2,3) = u + a*nx
R(3,3) = v + a*ny
R(4,3) = w + a*nz
R(5,3) = H + a*qn
! Two shear wave components combined into one (wave strength incorporated).
du = uR - uL
dv = vR - vL
dw = wR - wL
R(1,4) = zero
R(2,4) = du - dqn*nx
R(3,4) = dv - dqn*ny
R(4,4) = dw - dqn*nz
R(5,4) = u*du + v*dv + w*dw - qn*dqn
!Dissipation Term: |An|(UR-UL) = R|Lambda|L*dU = sum_k of [ ws(k) * R(:,k) * L*dU(k) ]
diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) &
+ ws(3)*LdU(3)*R(:,3) + ws(4)*LdU(4)*R(:,4)
!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)
fL(1) = rhoL*qnL
fL(2) = rhoL*qnL * uL + pL*nx
fL(3) = rhoL*qnL * vL + pL*ny
fL(4) = rhoL*qnL * wL + pL*nz
fL(5) = rhoL*qnL * HL
fR(1) = rhoR*qnR
fR(2) = rhoR*qnR * uR + pR*nx
fR(3) = rhoR*qnR * vR + pR*ny
fR(4) = rhoR*qnR * wR + pR*nz
fR(5) = rhoR*qnR * HR
! This is the numerical flux: Roe flux = 1/2 *[ Fn(UL)+Fn(UR) - |An|(UR-UL) ]
num_flux = half * (fL + fR - diss)
!Normal max wave speed in the normal direction.
! wsn = abs(qn) + a
end subroutine inviscid_roe_n
!--------------------------------------------------------------------------------
!********************************************************************************
!* -- 3D Rotated-Roe-HLL (Rotated-RHLL) Flux Function ---
!*
!* This subroutine computes the Rotated-RHLL flux for the Euler equations
!* in the direction, njk=[nx,ny,nz].
!*
!* H. Nishikawa and K. Kitamura, Very Simple, Carbuncle-Free, Boundary-Layer
!* Resolving, Rotated-Hybrid Riemann Solvers,
!* Journal of Computational Physics, 227, pp. 2560-2581, 2008.
!*
!* Very robust for nonlinear instability (carbuncle).
!* Recommended for high-speed flows involving strong shocks.
!* It is also capable of resolving boundary layers.
!*
!*
!* Conservative form of the Euler equations:
!*
!* dU/dt + dF/dx + dG/dy + dH/dz = 0
!*
!* This subroutine computes the numerical flux for the flux in the direction,
!* njk=[nx,ny,nz]:
!*
!* Fn = F*nx + G*ny + H*nz = | rho*qn |
!* | rho*qn*u + p*nx |
!* | rho*qn*v + p*ny |
!* | rho*qn*w + p*nz |
!* | rho*qn*H | (qn = u*nx + v*ny + w*nz)
!*
!* The Rotated-RHLL flux is defined as
!*
!* Numerical flux = alpha1*HLL(n1) + alpha2*Roe(n2), if |dq| > eps
!* Roe(n) , if |dq| < eps
!*
!* where n1 is taken as the velocity difference vector (normal to shock or
!* tangent to shear waves), n2 is a vector perpendicular to n1,
!* alpha1 = n*n1, alpha2 = n*n2. That is, we decompose the normal vector, n,
!* in the two directions, n1 and n2, and apply the HLL in n1 and the Roe in n2.
!* However, the resulting flux can be simplified and can be made to look like
!* a single Roe-type flux. The additional cost over the Roe flux is reported
!* to be as small as 14% of the Roe flux.
!*
!* Note: HLL is introduced only near discontinuities and only by a fraction, alpha1.
!* Note: For full convergence to machine zero for staedy state computations,
!* the vectors, n1 and n2, may have to be frozen. For time-accurate
! calculation, it is not necessary.
!* Note: Here, the Roe flux is computed without tangent vectors.
!*
!* ------------------------------------------------------------------------------
!* Input: primL(1:5) = Left state (rhoL, uL, vL, wR, pL)
!* primR(1:5) = Right state (rhoR, uR, vR, wR, pR)
!* njk(1:3) = unit face normal vector (nx, ny, nz), pointing from Left to Right
!*
!* njk
!* Face normal ^ o Right data point
!* | .
!* | .
!* |.
!* -------x-------- Face
!* . Left and right states are
!* . 1. Values at data points for 1st-order accuracy
!* . 2. Extrapolated values at the face midpoint 'x'
!* o Left data point for 2nd/higher-order accuracy.
!*
!* Output: flux(1:5) = the Roe flux
!* ------------------------------------------------------------------------------
!*
!* Note: This subroutine has been prepared for an educational purpose.
!* It is not at all efficient. Think about how you can optimize it.
!*
!* Note: Please let me know if you find bugs. I'll greatly appreciate it and
!* fix the bugs.
!*
!* Katate Masatsuka, March 2011. http://www.cfdbooks.com
!********************************************************************************
subroutine inviscid_rotated_rhll(primL, primR, njk, num_flux)
implicit none
integer , parameter :: p2 = selected_real_kind(15) ! Double precision
!Input
real(p2), intent( in) :: primL(5), primR(5) ! Input: primitive variables
real(p2), intent( in) :: njk(3) ! Input: face normal vector
!Output
real(p2), intent(out) :: num_flux(5) ! Output: numerical flux
!Some constants
real(p2) :: zero = 0.0_p2
real(p2) :: one = 1.0_p2
real(p2) :: two = 2.0_p2
real(p2) :: half = 0.5_p2
real(p2) :: fifth = 0.2_p2
!Local variables
real(p2) :: nx, ny, nz ! Face normal vector
real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components.
real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables.
real(p2) :: qnL, qnR ! Normal velocity
real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy
real(p2) :: RT,rho,u,v,w,H,a,qn ! Roe-averages
real(p2) :: drho,dqn,dp,LdU(4) ! Wave strengths
real(p2) :: du, dv, dw ! Velocity conponent differences
real(p2) :: eig(4) ! Eigenvalues
real(p2) :: ws(4), R(5,4) ! Absolute Wave speeds and right-eigenvectors
real(p2) :: dws(4) ! Width of a parabolic fit for entropy fix
real(p2) :: fL(5), fR(5), diss(5) ! Fluxes ad dissipation term
real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats
real(p2) :: SRp,SLm ! Wave speeds for the HLL part
real(p2) :: nx1, ny1, nz1 ! Vector along which HLL is applied
real(p2) :: nx2, ny2, nz2 ! Vector along which Roe is applied
real(p2) :: alpha1, alpha2 ! Projections of the new normals
real(p2) :: abs_dq ! Magnitude of the velocity difference
real(p2) :: temp, tempx, tempy, tempz ! Temporary variables
! Face normal vector (unit vector)
nx = njk(1)
ny = njk(2)
nz = njk(3)
!Primitive and other variables.
! Left state
rhoL = primL(1)
uL = primL(2)
vL = primL(3)
wL = primL(4)
qnL = uL*nx + vL*ny + wL*nz
pL = primL(5)
aL = sqrt(gamma*pL/rhoL)
HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
! Right state
rhoR = primR(1)
uR = primR(2)
vR = primR(3)
wR = primR(4)
qnR = uR*nx + vR*ny + wR*nz
pR = primR(5)
aR = sqrt(gamma*pR/rhoR)
HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)
!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)
fL(1) = rhoL*qnL
fL(2) = rhoL*qnL * uL + pL*nx
fL(3) = rhoL*qnL * vL + pL*ny
fL(4) = rhoL*qnL * wL + pL*nz
fL(5) = rhoL*qnL * HL
fR(1) = rhoR*qnR
fR(2) = rhoR*qnR * uR + pR*nx
fR(3) = rhoR*qnR * vR + pR*ny
fR(4) = rhoR*qnR * wR + pR*nz
fR(5) = rhoR*qnR * HR
!--------------------------------------------------------------------------------
!Define n1 and n2, and compute alpha1 and alpha2: (4.2) in the original paper.
! Note: n1 and n2 may need to be frozen at some point during
! a steady calculation to fully make it converge. For time-accurate
! calculation, it is not necessary.
! Note: For a boundary face, you may want to set (nx2,ny2,nz2)=(nx,ny,nz),
! (nx1,ny1)=tangent vector, i.e., use the Roe flux.
abs_dq = sqrt( (uR-uL)**2 + (vR-vL)**2 + (wR-wL)**2 )
!n1 = Velocity difference vector: normal to shock or tangent to shear
if ( abs_dq > 1.0e-12_p2) then
nx1 = (uR-uL)/abs_dq
ny1 = (vR-vL)/abs_dq
nz1 = (wR-wL)/abs_dq
!n1 = Face tangent vector if abs_dq is too small.
! Note: There are infinitely many choices for the tangent vector.
! The best choice may be discovered in future.
! Here, we choose a vector in a plane (essentially 2D vector).
! Note that we must be careful and make sure that the vector is not a zero vector.
else
! E.g., if n = (1,0,0), then, (0,-nz,ny) = (0,0,0). This is a problem.
! To avoid such a situation, we choose the 2D tangential vector
! having the maximum length.
tempx = ny*ny + nz*nz
tempy = nz*nz + nx*nx
tempz = nx*nx + ny*ny
if ( tempx >= tempy .and. tempx >= tempz ) then
nx1 = zero
ny1 = -nz
nz1 = ny
elseif ( tempy >= tempx .and. tempy >= tempz ) then
nx1 = -nz
ny1 = zero
nz1 = nx
elseif ( tempz >= tempx .and. tempz >= tempy ) then
nx1 = -ny
ny1 = nx
nz1 = zero
else
! Impossible to happen
write(*,*) "inviscid_rotated_rhll: Impossible to happen. Please report the problem."
stop
endif
! Make it the unit vector.
temp = sqrt( nx1*nx1 + ny1*ny1 + nz1*nz1 )
nx1 = nx1/temp
ny1 = ny1/temp
nz1 = nz1/temp
endif
alpha1 = nx*nx1 + ny*ny1 + nz*nz1
! Make alpha1 always positive.
temp = sign(one,alpha1)
nx1 = temp * nx1
ny1 = temp * ny1
nz1 = temp * nz1
alpha1 = temp * alpha1
!n2 = direction perpendicular to n1.
! Note: There are infinitely many choices for this vector.
! The best choice may be discovered in future.
! Here, we employ the formula (4.4) in the paper:
! (nx2,ny2,nz2) = (n1xn)xn1 / |(n1xn)xn1| ('x' is the vector product.)
! (tempx,tempy,tempz) = n1xn
tempx = ny1*nz - nz1*ny
tempy = nz1*nx - nx1*nz
tempz = nx1*ny - ny1*nx
! (nx2,ny2,nz2) = (n1xn)xn1
nx2 = tempy*nz1 - tempz*ny1
ny2 = tempz*nx1 - tempx*nz1
nz2 = tempx*ny1 - tempy*nx1
! Make n2 the unit vector
temp = sqrt( nx2*nx2 + ny2*ny2 + nz2*nz2 )
nx2 = nx2/temp
ny2 = ny2/temp
nz2 = nz2/temp
alpha2 = nx*nx2 + ny*ny2 + nz*nz2
! Make alpha2 always positive.
temp = sign(one,alpha2)
nx2 = temp * nx2
ny2 = temp * ny2
nz2 = temp * nz2
alpha2 = temp * alpha2
!--------------------------------------------------------------------------------
!Now we are going to compute the Roe flux with n2 as the normal with modified
!wave speeds (5.12). NOTE: the Roe flux here is computed without tangent vectors.
!See "I do like CFD, VOL.1" for details: page 57, Equation (3.6.31).
!First compute the Roe-averaged quantities
! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for
! the Roe-averaged density.
RT = sqrt(rhoR/rhoL)
rho = RT*rhoL !Roe-averaged density.
u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity
v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity
w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity
H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy
a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound
!----------------------------------------------------
!Compute the wave speed estimates for the HLL part,
!following Einfeldt:
!
! B. Einfeldt, On Godunov-type methods for gas dynamics,
! SIAM Journal on Numerical Analysis 25 (2) (1988) 294–318.
!
! Note: HLL is actually applied to n1, but this is
! all we need to incorporate HLL. See JCP2008 paper.
qn = u *nx1 + v *ny1 + w *nz1
qnL = uL*nx1 + vL*ny1 + wL*nz1
qnR = uR*nx1 + vR*ny1 + wR*nz1
SLm = min( zero, qn - a, qnL - aL ) !Minimum wave speed estimate
SRp = max( zero, qn + a, qnR + aR ) !Maximum wave speed estimate
! This is the only place where n1=(nx1,ny1,nz1) is used.
! n1=(nx1,ny1,nz1) is never used below.
!----------------------------------------------------
!Wave Strengths
qn = u *nx2 + v *ny2 + w *nz2
qnL = uL*nx2 + vL*ny2 + wL*nz2
qnR = uR*nx2 + vR*ny2 + wR*nz2
drho = rhoR - rhoL !Density difference
dp = pR - pL !Pressure difference
dqn = qnR - qnL !Normal velocity difference
LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength
LdU(2) = drho - dp/(a*a) !Entropy wave strength
LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength
LdU(4) = rho !Shear wave strength (not really, just a factor)
!Wave Speed (Eigenvalues)
eig(1) = qn-a !Left-moving acoustic wave velocity
eig(2) = qn !Entropy wave velocity
eig(3) = qn+a !Right-moving acoustic wave velocity
eig(4) = qn !Shear wave velocity
!Absolute values of the wave speeds (Eigenvalues)
ws(1) = abs(qn-a) !Left-moving acoustic wave speed
ws(2) = abs(qn) !Entropy wave speed
ws(3) = abs(qn+a) !Right-moving acoustic wave speed
ws(4) = abs(qn) !Shear wave speed
!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.
!NOTE: It avoids vanishing wave speeds by making a parabolic fit near ws = 0.
dws(1) = fifth
if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
dws(3) = fifth
if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )
!Combine the wave speeds for Rotated-RHLL: Eq.(5.12) in the original JCP2008 paper.
ws = alpha2*ws - (alpha1*two*SRp*SLm + alpha2*(SRp+SLm)*eig)/(SRp-SLm)
!Below, we compute the Roe dissipation term in the direction n2
!with the above modified wave speeds. HLL wave speeds act something like
!the entropy fix or eigenvalue limiting; they contribute only by the amount
!given by the fraction, alpha1 (less than or equal to 1.0). See JCP2008 paper.
!Right Eigenvectors:
!Note: Two shear wave components are combined into one, so that tangent vectors
! are not required. And that's why there are only 4 vectors here.
! Left-moving acoustic wave
R(1,1) = one
R(2,1) = u - a*nx2
R(3,1) = v - a*ny2
R(4,1) = w - a*nz2
R(5,1) = H - a*qn
! Entropy wave
R(1,2) = one
R(2,2) = u
R(3,2) = v
R(4,2) = w
R(5,2) = half*(u*u + v*v + w*w)
! Right-moving acoustic wave
R(1,3) = one
R(2,3) = u + a*nx2
R(3,3) = v + a*ny2
R(4,3) = w + a*nz2
R(5,3) = H + a*qn
! Two shear wave components combined into one (wave strength incorporated).
du = uR - uL
dv = vR - vL
dw = wR - wL
R(1,4) = zero
R(2,4) = du - dqn*nx2
R(3,4) = dv - dqn*ny2
R(4,4) = dw - dqn*nz2
R(5,4) = u*du + v*dv + w*dw - qn*dqn
!Dissipation Term: Roe dissipation with the modified wave speeds.
! |An|dU = R|Lambda|L*dU = sum_k of [ ws(k) * R(:,k) * L*dU(k) ], where n=n2.
diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) &
+ ws(3)*LdU(3)*R(:,3) + ws(4)*LdU(4)*R(:,4)
!Compute the Rotated-RHLL flux. (It looks like the HLL flux with Roe dissipation.)
num_flux = (SRp*fL - SLm*fR)/(SRp-SLm) - half*diss
!Normal max wave speed in the normal direction.
! wsn = abs(qn) + a
end subroutine inviscid_rotated_rhll
!--------------------------------------------------------------------------------