Update the RANS (SA) equation with LU-SGS
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(inout), | dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var) | :: | qp | Store primitive variable at cell center |
|
| real(kind=wp), | intent(in), | dimension(:, :, :, :) | :: | residue | Store residue at each cell-center |
|
| real(kind=wp), | intent(in), | dimension(1:dims%imx-1, 1:dims%jmx-1, 1:dims%kmx-1) | :: | delta_t | Local time increment value at each cell center |
|
| type(celltype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2) | :: | cells | Input cell quantities: volume |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+3,-2:dims%jmx+2,-2:dims%kmx+2) | :: | Ifaces | Input varaible which stores I faces' area and unit normal |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2) | :: | Jfaces | Input varaible which stores J faces' area and unit normal |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3) | :: | Kfaces | Input varaible which stores K faces' area and unit normal |
|
| type(extent), | intent(in) | :: | dims | Extent of the domain:imx,jmx,kmx |
subroutine update_SA_variables(qp, residue, delta_t, cells, Ifaces, Jfaces, Kfaces, dims)
!< Update the RANS (SA) equation with LU-SGS
implicit none
type(extent), intent(in) :: dims
!< Extent of the domain:imx,jmx,kmx
real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(inout) :: qp
!< Store primitive variable at cell center
real(wp), dimension(:, :, :, :), intent(in) :: residue
!< Store residue at each cell-center
real(wp) , dimension(1:dims%imx-1, 1:dims%jmx-1, 1:dims%kmx-1), intent(in) :: delta_t
!< Local time increment value at each cell center
type(celltype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: cells
!< Input cell quantities: volume
type(facetype), dimension(-2:dims%imx+3,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: Ifaces
!< Input varaible which stores I faces' area and unit normal
type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2), intent(in) :: Jfaces
!< Input varaible which stores J faces' area and unit normal
type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3), intent(in) :: Kfaces
!< Input varaible which stores K faces' area and unit normal
integer :: i,j,k
real(wp), dimension(1:6) :: deltaU
real(wp), dimension(1:6) :: D
real(wp), dimension(1:6) :: conservativeQ
real(wp), dimension(1:6) :: OldIminusFlux
real(wp), dimension(1:6) :: OldJminusFlux
real(wp), dimension(1:6) :: OldKminusFlux
real(wp), dimension(1:6) :: NewIminusFlux
real(wp), dimension(1:6) :: NewJminusFlux
real(wp), dimension(1:6) :: NewKminusFlux
real(wp), dimension(1:6) :: DelIminusFlux
real(wp), dimension(1:6) :: DelJminusFlux
real(wp), dimension(1:6) :: DelKminusFlux
real(wp), dimension(1:6) :: LambdaTimesArea
real(wp), dimension(1:6) :: Q0 ! state at cell
real(wp), dimension(1:6) :: Q1 ! state at neighbours
real(wp), dimension(1:6) :: Q2
real(wp), dimension(1:6) :: Q3
real(wp), dimension(1:6) :: Q4
real(wp), dimension(1:6) :: Q5
real(wp), dimension(1:6) :: Q6
real(wp), dimension(1:6) :: DQ0! change in state
real(wp), dimension(1:6) :: DQ1
real(wp), dimension(1:6) :: DQ2
real(wp), dimension(1:6) :: DQ3
real(wp), dimension(1:6) :: DQ4
real(wp), dimension(1:6) :: DQ5
real(wp), dimension(1:6) :: DQ6
real(wp), dimension(1:7) :: Flist1
real(wp), dimension(1:7) :: Flist2
real(wp), dimension(1:7) :: Flist3
real(wp), dimension(1:7) :: Flist4
real(wp), dimension(1:7) :: Flist5
real(wp), dimension(1:7) :: Flist6
real(wp), dimension(1:3) :: C0
real(wp), dimension(1:3) :: C1
real(wp), dimension(1:3) :: C2
real(wp), dimension(1:3) :: C3
real(wp), dimension(1:3) :: C4
real(wp), dimension(1:3) :: C5
real(wp), dimension(1:3) :: C6
real(wp) :: eps
real(wp) :: M
real(wp) :: VMag
real(wp) :: SoundMag
real(wp) :: u,v,w,p,H,tv!r
real(wp) :: factor
real(wp), dimension(1:6,1:6) :: PrecondInv
real(wp) :: fv1
real(wp) :: fv2
real(wp) :: fw
real(wp) :: g
real(wp) :: r
real(wp) :: dist_i
real(wp) :: dist_i_2
real(wp) :: Ji
real(wp) :: Ji_2
real(wp) :: Ji_3
real(wp) :: S
real(wp) :: Omega
real(wp) :: k2
real(wp) :: inv_k2_d2
real(wp) :: Shat
real(wp) :: inv_Shat
real(wp) :: nu
real(wp) :: glim
real(wp) :: g_6
real(wp) :: dfv1
real(wp) :: dfv2
real(wp) :: dfw
real(wp) :: dShat
real(wp) :: dr
real(wp) :: dg
real(wp) :: density
!intialize delQ
delQstar = 0.0
!forward sweep
do k=1,dims%kmx-1
do j=1,dims%jmx-1
do i=1,dims%imx-1
density = qp(i,j,k,1)
C0 = (/Cells(i ,j ,k )%Centerx,Cells(i ,j ,k )%Centery,Cells(i ,j ,k )%Centerz/)
C1 = (/Cells(i-1,j ,k )%Centerx,Cells(i-1,j ,k )%Centery,Cells(i-1,j ,k )%Centerz/)
C2 = (/Cells(i ,j-1,k )%Centerx,Cells(i ,j-1,k )%Centery,Cells(i ,j-1,k )%Centerz/)
C3 = (/Cells(i ,j ,k-1)%Centerx,Cells(i ,j ,k-1)%Centery,Cells(i ,j ,k-1)%Centerz/)
C4 = (/Cells(i+1,j ,k )%Centerx,Cells(i+1,j ,k )%Centery,Cells(i+1,j ,k )%Centerz/)
C5 = (/Cells(i ,j+1,k )%Centerx,Cells(i ,j+1,k )%Centery,Cells(i ,j+1,k )%Centerz/)
C6 = (/Cells(i ,j ,k+1)%Centerx,Cells(i ,j ,k+1)%Centery,Cells(i ,j ,k+1)%Centerz/)
Q0 = qp(i , j , k , 1:6)
Q1 = qp(i-1, j , k , 1:6)
Q2 = qp(i , j-1, k , 1:6)
Q3 = qp(i , j , k-1, 1:6)
Q4 = qp(i+1, j , k , 1:6)
Q5 = qp(i , j+1, k , 1:6)
Q6 = qp(i , j , k+1, 1:6)
DQ0 = 0.0
DQ1 = delQstar(i-1, j , k , 1:6)
DQ2 = delQstar(i , j-1, k , 1:6)
DQ3 = delQstar(i , j , k-1, 1:6)
Flist1(1) = Ifaces(i,j,k)%A
Flist1(2) = -Ifaces(i,j,k)%nx
Flist1(3) = -Ifaces(i,j,k)%ny
Flist1(4) = -Ifaces(i,j,k)%nz
Flist1(5) = 0.5*(cells(i-1, j , k )%volume + cells(i,j,k)%volume)
Flist1(6) = 0.5*( mmu(i-1, j , k ) + mmu(i,j,k))
Flist1(7) = 0.5*( tmu(i-1, j , k ) + tmu(i,j,k))
Flist2(1) = Jfaces(i,j,k)%A
Flist2(2) = -Jfaces(i,j,k)%nx
Flist2(3) = -Jfaces(i,j,k)%ny
Flist2(4) = -Jfaces(i,j,k)%nz
Flist2(5) = 0.5*(cells(i , j-1, k )%volume + cells(i,j,k)%volume)
Flist2(6) = 0.5*( mmu(i , j-1, k ) + mmu(i,j,k))
Flist2(7) = 0.5*( tmu(i , j-1, k ) + tmu(i,j,k))
Flist3(1) = Kfaces(i,j,k)%A
Flist3(2) = -Kfaces(i,j,k)%nx
Flist3(3) = -Kfaces(i,j,k)%ny
Flist3(4) = -Kfaces(i,j,k)%nz
Flist3(5) = 0.5*(cells(i , j , k-1)%volume + cells(i,j,k)%volume)
Flist3(6) = 0.5*( mmu(i , j , k-1) + mmu(i,j,k))
Flist3(7) = 0.5*( tmu(i , j , k-1) + tmu(i,j,k))
Flist4(1) = Ifaces(i+1,j,k)%A
Flist4(2) = +Ifaces(i+1,j,k)%nx
Flist4(3) = +Ifaces(i+1,j,k)%ny
Flist4(4) = +Ifaces(i+1,j,k)%nz
Flist4(5) = 0.5*(cells(i+1, j , k )%volume + cells(i,j,k)%volume)
Flist4(6) = 0.5*( mmu(i+1, j , k ) + mmu(i,j,k))
Flist4(7) = 0.5*( tmu(i+1, j , k ) + tmu(i,j,k))
Flist5(1) = Jfaces(i,j+1,k)%A
Flist5(2) = +Jfaces(i,j+1,k)%nx
Flist5(3) = +Jfaces(i,j+1,k)%ny
Flist5(4) = +Jfaces(i,j+1,k)%nz
Flist5(5) = 0.5*(cells(i , j+1, k )%volume + cells(i,j,k)%volume)
Flist5(6) = 0.5*( mmu(i , j+1, k ) + mmu(i,j,k))
Flist5(7) = 0.5*( tmu(i , j+1, k ) + tmu(i,j,k))
Flist6(1) = Kfaces(i,j,k+1)%A
Flist6(2) = +Kfaces(i,j,k+1)%nx
Flist6(3) = +Kfaces(i,j,k+1)%ny
Flist6(4) = +Kfaces(i,j,k+1)%nz
Flist6(5) = 0.5*(cells(i , j , k+1)%volume + cells(i,j,k)%volume)
Flist6(6) = 0.5*( mmu(i , j , k+1) + mmu(i,j,k))
Flist6(7) = 0.5*( tmu(i , j , k+1) + tmu(i,j,k))
NewIminusFlux = SAFlux(Q1, Q0, DQ1, Flist1)
NewJminusFlux = SAFlux(Q2, Q0, DQ2, Flist2)
NewKminusFlux = SAFlux(Q3, Q0, DQ3, Flist3)
OldIminusFlux = SAFlux(Q1, Q0, DQ0, Flist1)
OldJminusFlux = SAFlux(Q2, Q0, DQ0, Flist2)
OldKminusFlux = SAFlux(Q3, Q0, DQ0, Flist3)
!---preconditioning---
r = Q0(1)
u = Q0(2)
v = Q0(3)
w = Q0(4)
p = Q0(5)
tv = Q0(6)
VMag = sqrt(u*u + v*v + w*w)
SoundMag = sqrt(gm*p/r)
M = VMag/SoundMag
H = (gm*p/(r*(gm-1.0))) + 0.5*(VMag)
eps = min(1.0, max(M*M, Minf*Minf))
factor = (1.0-eps)*(gm-1.0)/(SoundMag*SoundMag)
LambdaTimesArea(1)= SpectralRadius(Q1, Q0, Flist1, C1, C0, eps)
LambdaTimesArea(2)= SpectralRadius(Q2, Q0, Flist2, C2, C0, eps)
LambdaTimesArea(3)= SpectralRadius(Q3, Q0, Flist3, C3, C0, eps)
LambdaTimesArea(4)= SpectralRadius(Q4, Q0, Flist4, C4, C0, eps)
LambdaTimesArea(5)= SpectralRadius(Q5, Q0, Flist5, C5, C0, eps)
LambdaTimesArea(6)= SpectralRadius(Q6, Q0, Flist6, C6, C0, eps)
PrecondInv(1,1) = 1.0 - factor*1*VMag*VMag/2.0
PrecondInv(2,1) = 0.0 - factor*u*VMag*VMag/2.0
PrecondInv(3,1) = 0.0 - factor*v*VMag*VMag/2.0
PrecondInv(4,1) = 0.0 - factor*w*VMag*VMag/2.0
PrecondInv(5,1) = 0.0 - factor*H*VMag*VMag/2.0
PrecondInv(6,1) = 0.0 - factor*tv*VMag*VMag/2.0
PrecondInv(1,2) = 0.0 - factor*1*(-u)
PrecondInv(2,2) = 1.0 - factor*u*(-u)
PrecondInv(3,2) = 0.0 - factor*v*(-u)
PrecondInv(4,2) = 0.0 - factor*w*(-u)
PrecondInv(5,2) = 0.0 - factor*H*(-u)
PrecondInv(6,2) = 0.0 - factor*tv*(-u)
PrecondInv(1,3) = 0.0 - factor*1*(-v)
PrecondInv(2,3) = 0.0 - factor*u*(-v)
PrecondInv(3,3) = 1.0 - factor*v*(-v)
PrecondInv(4,3) = 0.0 - factor*w*(-v)
PrecondInv(5,3) = 0.0 - factor*H*(-v)
PrecondInv(6,3) = 0.0 - factor*tv*(-v)
PrecondInv(1,4) = 0.0 - factor*1*(-w)
PrecondInv(2,4) = 0.0 - factor*u*(-w)
PrecondInv(3,4) = 0.0 - factor*v*(-w)
PrecondInv(4,4) = 1.0 - factor*w*(-w)
PrecondInv(5,4) = 0.0 - factor*H*(-w)
PrecondInv(6,4) = 0.0 - factor*tv*(-w)
PrecondInv(1,5) = 0.0 - factor*1*(1.)
PrecondInv(2,5) = 0.0 - factor*u*(1.)
PrecondInv(3,5) = 0.0 - factor*v*(1.)
PrecondInv(4,5) = 0.0 - factor*w*(1.)
PrecondInv(5,5) = 1.0 - factor*H*(1.)
PrecondInv(6,5) = 0.0 - factor*tv*(1.)
PrecondInv(1,6) = 0.0 - factor*1*(-1.)
PrecondInv(2,6) = 0.0 - factor*u*(-1.)
PrecondInv(3,6) = 0.0 - factor*v*(-1.)
PrecondInv(4,6) = 0.0 - factor*w*(-1.)
PrecondInv(5,6) = 0.0 - factor*H*(-1.)
PrecondInv(6,6) = 1.0 - factor*tv*(-1.)
!---end preconditioning
! multiply above flux with area to get correct values
DelIminusFlux = NewIminusFlux - OldIminusFlux
DelJminusFlux = NewJminusFlux - OldJminusFlux
DelKminusFlux = NewKminusFlux - OldKminusFlux
D = (cells(i,j,k)%volume/delta_t(i,j,k)) + 0.5*SUM(LambdaTimesArea)
!storing D in Iflux array for backward sweep
!F_p(i,j,k,1) = D
! -- source term derivatives -- !
Omega = sqrt( ((gradw_y(i,j,k)- gradv_z(i,j,k))**2 &
+ (gradu_z(i,j,k)- gradw_x(i,j,k))**2 &
+ (gradv_x(i,j,k)- gradu_y(i,j,k))**2 &
)&
)
dist_i = dist(i,j,k)
dist_i_2 = dist_i*dist_i
k2 = kappa_sa*kappa_sa
nu = mu(i,j,k)/density
Ji = Q0(6)/nu
Ji_2 = Ji*Ji
Ji_3 = Ji_2*ji
! ___ functions ___
fv1 = (Ji_3)/((Ji_3) + (cv1_3))
fv2 = 1.0 - Ji/(1.0 + (Ji*fv1))
! ___ Shear stress for production ___
S = Omega
inv_k2_d2 = 1.0/(k2*dist_i_2)
Shat = S + Q0(6)*fv2*inv_k2_d2
Shat = max(Shat, 1.0e-10)
inv_Shat = 1.0/Shat
dfv1 = 3.0*Ji_2*cv1_3/(nu*(Ji_3+cv1_3)**2)
dfv2 = -((1.0/nu) - Ji_2*dfv1)/((1.0+Ji*fv1)**2)
dShat = (fv2+Q0(6)*dfv2)*inv_k2_d2
D = D - cb1*(Q0(6)*dShat+Shat)*cells(i,j,k)%volume
! ___ Destruction term___ !
r = min(Q0(6)*inv_Shat*inv_k2_d2, 10.0)
g = r + cw2*((r**6) - r)
g_6 = g**6
glim = ((1.0+cw3_6)/(g_6+cw3_6))**(1.0/6.0)
fw = g*glim
dr = (Shat-Q0(6)*dShat)*inv_Shat*inv_Shat*inv_k2_d2
dg = dr*(1.0+cw2*(6.0*(r**5)-1.0))
dfw= dg*glim*(1.0-g_6/(g_6+cw3_6))
D = D+cw1*(dfw*Q0(6) + 2*fw)*Q0(6)/dist_i_2*cells(i,j,k)%volume
! -- end of source term -- !
deltaU(1:6) = -matmul(PrecondInv,residue(i,j,k,1:6)) &
- 0.5*((matmul(PrecondInv,DelIminusFlux) - LambdaTimesArea(1)*delQstar(i-1,j,k,1:6)) &
+ (matmul(PrecondInv,DelJminusFlux) - LambdaTimesArea(2)*delQstar(i,j-1,k,1:6)) &
+ (matmul(PrecondInv,DelKminusFlux) - LambdaTimesArea(3)*delQstar(i,j,k-1,1:6)) )
delQstar(i,j,k,1:6) = deltaU(1:6)/D
end do
end do
end do
delQ=0.0
!backward sweep
do i=dims%imx-1,1,-1
do j=dims%jmx-1,1,-1
do k=dims%kmx-1,1,-1
density = qp(i,j,k,1)
C0 = (/Cells(i ,j ,k )%Centerx,Cells(i ,j ,k )%Centery,Cells(i ,j ,k )%Centerz/)
C1 = (/Cells(i-1,j ,k )%Centerx,Cells(i-1,j ,k )%Centery,Cells(i-1,j ,k )%Centerz/)
C2 = (/Cells(i ,j-1,k )%Centerx,Cells(i ,j-1,k )%Centery,Cells(i ,j-1,k )%Centerz/)
C3 = (/Cells(i ,j ,k-1)%Centerx,Cells(i ,j ,k-1)%Centery,Cells(i ,j ,k-1)%Centerz/)
C4 = (/Cells(i+1,j ,k )%Centerx,Cells(i+1,j ,k )%Centery,Cells(i+1,j ,k )%Centerz/)
C5 = (/Cells(i ,j+1,k )%Centerx,Cells(i ,j+1,k )%Centery,Cells(i ,j+1,k )%Centerz/)
C6 = (/Cells(i ,j ,k+1)%Centerx,Cells(i ,j ,k+1)%Centery,Cells(i ,j ,k+1)%Centerz/)
Q0 = qp(i , j , k , 1:6)
Q1 = qp(i-1, j , k , 1:6)
Q2 = qp(i , j-1, k , 1:6)
Q3 = qp(i , j , k-1, 1:6)
Q4 = qp(i+1, j , k , 1:6)
Q5 = qp(i , j+1, k , 1:6)
Q6 = qp(i , j , k+1, 1:6)
DQ0 = 0.0
DQ4 = delQ(i+1, j , k , 1:6)
DQ5 = delQ(i , j+1, k , 1:6)
DQ6 = delQ(i , j , k+1, 1:6)
Flist1(1) = Ifaces(i,j,k)%A
Flist1(2) = -Ifaces(i,j,k)%nx
Flist1(3) = -Ifaces(i,j,k)%ny
Flist1(4) = -Ifaces(i,j,k)%nz
Flist1(5) = 0.5*(cells(i-1, j , k )%volume + cells(i,j,k)%volume)
Flist1(6) = 0.5*( mmu(i-1, j , k ) + mmu(i,j,k))
Flist1(7) = 0.5*( tmu(i-1, j , k ) + tmu(i,j,k))
Flist2(1) = Jfaces(i,j,k)%A
Flist2(2) = -Jfaces(i,j,k)%nx
Flist2(3) = -Jfaces(i,j,k)%ny
Flist2(4) = -Jfaces(i,j,k)%nz
Flist2(5) = 0.5*(cells(i , j-1, k )%volume + cells(i,j,k)%volume)
Flist2(6) = 0.5*( mmu(i , j-1, k ) + mmu(i,j,k))
Flist2(7) = 0.5*( tmu(i , j-1, k ) + tmu(i,j,k))
Flist3(1) = Kfaces(i,j,k)%A
Flist3(2) = -Kfaces(i,j,k)%nx
Flist3(3) = -Kfaces(i,j,k)%ny
Flist3(4) = -Kfaces(i,j,k)%nz
Flist3(5) = 0.5*(cells(i , j , k-1)%volume + cells(i,j,k)%volume)
Flist3(6) = 0.5*( mmu(i , j , k-1) + mmu(i,j,k))
Flist3(7) = 0.5*( tmu(i , j , k-1) + tmu(i,j,k))
Flist4(1) = Ifaces(i+1,j,k)%A
Flist4(2) = +Ifaces(i+1,j,k)%nx
Flist4(3) = +Ifaces(i+1,j,k)%ny
Flist4(4) = +Ifaces(i+1,j,k)%nz
Flist4(5) = 0.5*(cells(i+1, j , k )%volume + cells(i,j,k)%volume)
Flist4(6) = 0.5*( mmu(i+1, j , k ) + mmu(i,j,k))
Flist4(7) = 0.5*( tmu(i+1, j , k ) + tmu(i,j,k))
Flist5(1) = Jfaces(i,j+1,k)%A
Flist5(2) = +Jfaces(i,j+1,k)%nx
Flist5(3) = +Jfaces(i,j+1,k)%ny
Flist5(4) = +Jfaces(i,j+1,k)%nz
Flist5(5) = 0.5*(cells(i , j+1, k )%volume + cells(i,j,k)%volume)
Flist5(6) = 0.5*( mmu(i , j+1, k ) + mmu(i,j,k))
Flist5(7) = 0.5*( tmu(i , j+1, k ) + tmu(i,j,k))
Flist6(1) = Kfaces(i,j,k+1)%A
Flist6(2) = +Kfaces(i,j,k+1)%nx
Flist6(3) = +Kfaces(i,j,k+1)%ny
Flist6(4) = +Kfaces(i,j,k+1)%nz
Flist6(5) = 0.5*(cells(i , j , k+1)%volume + cells(i,j,k)%volume)
Flist6(6) = 0.5*( mmu(i , j , k+1) + mmu(i,j,k))
Flist6(7) = 0.5*( tmu(i , j , k+1) + tmu(i,j,k))
NewIminusFlux = SAFlux(Q4, Q0, DQ4, Flist4)
NewJminusFlux = SAFlux(Q5, Q0, DQ5, Flist5)
NewKminusFlux = SAFlux(Q6, Q0, DQ6, Flist6)
OldIminusFlux = SAFlux(Q4, Q0, DQ0, Flist4)
OldJminusFlux = SAFlux(Q5, Q0, DQ0, Flist5)
OldKminusFlux = SAFlux(Q6, Q0, DQ0, Flist6)
!---preconditioning---
r = Q0(1)
u = Q0(2)
v = Q0(3)
w = Q0(4)
p = Q0(5)
tv = Q0(6)
VMag = sqrt(u*u + v*v + w*w)
SoundMag = sqrt(gm*p/r)
M = VMag/SoundMag
H = (gm*p/(r*(gm-1.0))) + 0.5*(VMag)
eps = min(1.0, max(M*M, Minf*Minf))
factor = (1.0-eps)*(gm-1.0)/(SoundMag*SoundMag)
LambdaTimesArea(1)= SpectralRadius(Q1, Q0, Flist1, C1, C0, eps)
LambdaTimesArea(2)= SpectralRadius(Q2, Q0, Flist2, C2, C0, eps)
LambdaTimesArea(3)= SpectralRadius(Q3, Q0, Flist3, C3, C0, eps)
LambdaTimesArea(4)= SpectralRadius(Q4, Q0, Flist4, C4, C0, eps)
LambdaTimesArea(5)= SpectralRadius(Q5, Q0, Flist5, C5, C0, eps)
LambdaTimesArea(6)= SpectralRadius(Q6, Q0, Flist6, C6, C0, eps)
PrecondInv(1,1) = 1.0 - factor*1*VMag*VMag/2.0
PrecondInv(2,1) = 0.0 - factor*u*VMag*VMag/2.0
PrecondInv(3,1) = 0.0 - factor*v*VMag*VMag/2.0
PrecondInv(4,1) = 0.0 - factor*w*VMag*VMag/2.0
PrecondInv(5,1) = 0.0 - factor*H*VMag*VMag/2.0
PrecondInv(6,1) = 0.0 - factor*tv*VMag*VMag/2.0
PrecondInv(1,2) = 0.0 - factor*1*(-u)
PrecondInv(2,2) = 1.0 - factor*u*(-u)
PrecondInv(3,2) = 0.0 - factor*v*(-u)
PrecondInv(4,2) = 0.0 - factor*w*(-u)
PrecondInv(5,2) = 0.0 - factor*H*(-u)
PrecondInv(6,2) = 0.0 - factor*tv*(-u)
PrecondInv(1,3) = 0.0 - factor*1*(-v)
PrecondInv(2,3) = 0.0 - factor*u*(-v)
PrecondInv(3,3) = 1.0 - factor*v*(-v)
PrecondInv(4,3) = 0.0 - factor*w*(-v)
PrecondInv(5,3) = 0.0 - factor*H*(-v)
PrecondInv(6,3) = 0.0 - factor*tv*(-v)
PrecondInv(1,4) = 0.0 - factor*1*(-w)
PrecondInv(2,4) = 0.0 - factor*u*(-w)
PrecondInv(3,4) = 0.0 - factor*v*(-w)
PrecondInv(4,4) = 1.0 - factor*w*(-w)
PrecondInv(5,4) = 0.0 - factor*H*(-w)
PrecondInv(6,4) = 0.0 - factor*tv*(-w)
PrecondInv(1,5) = 0.0 - factor*1*(1.)
PrecondInv(2,5) = 0.0 - factor*u*(1.)
PrecondInv(3,5) = 0.0 - factor*v*(1.)
PrecondInv(4,5) = 0.0 - factor*w*(1.)
PrecondInv(5,5) = 1.0 - factor*H*(1.)
PrecondInv(6,5) = 0.0 - factor*tv*(1.)
PrecondInv(1,6) = 0.0 - factor*1*(-1.)
PrecondInv(2,6) = 0.0 - factor*u*(-1.)
PrecondInv(3,6) = 0.0 - factor*v*(-1.)
PrecondInv(4,6) = 0.0 - factor*w*(-1.)
PrecondInv(5,6) = 0.0 - factor*H*(-1.)
PrecondInv(6,6) = 1.0 - factor*tv*(-1.)
!---end preconditioning
! multiply above flux with area to get correct values
DelIminusFlux = NewIminusFlux - OldIminusFlux
DelJminusFlux = NewJminusFlux - OldJminusFlux
DelKminusFlux = NewKminusFlux - OldKminusFlux
D = (cells(i,j,k)%volume/delta_t(i,j,k)) + 0.5*SUM(LambdaTimesArea)
! -- source term derivatives -- !
Omega = sqrt( ((gradw_y(i,j,k)- gradv_z(i,j,k))**2 &
+ (gradu_z(i,j,k)- gradw_x(i,j,k))**2 &
+ (gradv_x(i,j,k)- gradu_y(i,j,k))**2 &
)&
)
dist_i = dist(i,j,k)
dist_i_2 = dist_i*dist_i
k2 = kappa_sa*kappa_sa
nu = mu(i,j,k)/density
Ji = Q0(6)/nu
Ji_2 = Ji*Ji
Ji_3 = Ji_2*ji
! ___ functions ___
fv1 = (Ji_3)/((Ji_3) + (cv1_3))
fv2 = 1.0 - Ji/(1.0 + (Ji*fv1))
! ___ Shear stress for production ___
S = Omega
inv_k2_d2 = 1.0/(k2*dist_i_2)
Shat = S + Q0(6)*fv2*inv_k2_d2
Shat = max(Shat, 1.0e-10)
inv_Shat = 1.0/Shat
dfv1 = 3.0*Ji_2*cv1_3/(nu*(Ji_3+cv1_3)**2)
dfv2 = -((1.0/nu) - Ji_2*dfv1)/((1.0+Ji*fv1)**2)
dShat = (fv2+Q0(6)*dfv2)*inv_k2_d2
D = D - cb1*(Q0(6)*dShat+Shat)*cells(i,j,k)%volume
! ___ Destruction term___ !
r = min(Q0(6)*inv_Shat*inv_k2_d2, 10.0)
g = r + cw2*((r**6) - r)
g_6 = g**6
glim = ((1.0+cw3_6)/(g_6+cw3_6))**(1.0/6.0)
fw = g*glim
dr = (Shat-Q0(6)*dShat)*inv_Shat*inv_Shat*inv_k2_d2
dg = dr*(1.0+cw2*(6.0*(r**5)-1.0))
dfw= dg*glim*(1.0-g_6/(g_6+cw3_6))
D = D+cw1*(dfw*Q0(6) + 2*fw)*Q0(6)/dist_i_2*cells(i,j,k)%volume
! -- end of source term -- !
delQ(i,j,k,1:6) = delQstar(i,j,k,1:6) &
- 0.5*((matmul(PrecondInv,DelIminusFlux) - LambdaTimesArea(4)*delQ(i+1,j,k,1:6)) &
+ (matmul(PrecondInv,DelJminusFlux) - LambdaTimesArea(5)*delQ(i,j+1,k,1:6)) &
+ (matmul(PrecondInv,DelKminusFlux) - LambdaTimesArea(6)*delQ(i,j,k+1,1:6)) )/D
end do
end do
end do
do k=1,dims%kmx-1
do j = 1,dims%jmx-1
do i = 1,dims%imx-1
conservativeQ(1) = qp(i,j,k,1)
conservativeQ(2) = qp(i,j,k,1) * qp(i,j,k,2)
conservativeQ(3) = qp(i,j,k,1) * qp(i,j,k,3)
conservativeQ(4) = qp(i,j,k,1) * qp(i,j,k,4)
conservativeQ(5) = (qp(i,j,k,5) / (gm-1.0)) + ( 0.5 * qp(i,j,k,1) * sum( qp(i,j,k,2:4)**2) )
conservativeQ(6) = qp(i,j,k,1) * qp(i,j,k,6)
! add new change into conservative solution
conservativeQ(1:6) = conservativeQ(1:6) + delQ(i,j,k,1:6)
! convert back conservative to primitive
qp(i,j,k,1) = conservativeQ(1)
qp(i,j,k,2) = conservativeQ(2) / conservativeQ(1)
qp(i,j,k,3) = conservativeQ(3) / conservativeQ(1)
qp(i,j,k,4) = conservativeQ(4) / conservativeQ(1)
qp(i,j,k,5) = (gm-1.0) * ( conservativeQ(5) - (0.5 * sum(conservativeQ(2:4)**2) / conservativeQ(1)) )
qp(i,j,k,6) = conservativeQ(6) / conservativeQ(1)
qp(i,j,k,6) = max(qp(i,j,k,6), 1.e-8)
end do
end do
end do
end subroutine update_SA_variables