Subroutine to calculate state at the face, generalized for all direction : I,J, and K.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in), | 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(inout), | dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var) | :: | f_qp_left | primitive state at faces |
|
| real(kind=wp), | intent(inout), | dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var) | :: | f_qp_right | primitive state at faces |
|
| integer, | intent(in), | dimension(3) | :: | flags | flags for direction switch |
|
| type(celltype), | intent(in), | dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2) | :: | cells | Input cell quantities: volume |
|
| type(extent), | intent(in) | :: | dims | Extent of the domain:imx,jmx,kmx |
subroutine compute_face_states(qp, f_qp_left, f_qp_right, flags, cells, dims)
!< Subroutine to calculate state at the face, generalized for
!< all direction : I,J, and K.
implicit none
type(extent), intent(in) :: dims
!< Extent of the domain:imx,jmx,kmx
integer, dimension(3), intent(in) :: flags
!< flags for direction switch
real(wp), dimension(-2:dims%imx+2, -2:dims%jmx+2, -2:dims%kmx+2, 1:dims%n_var), intent(in):: qp
!< Store primitive variable at cell center
real(wp), dimension(1-flags(1):dims%imx-1+2*flags(1), 1-flags(2):dims%jmx-1+2*flags(2), 1-flags(3):dims%kmx-1+2*flags(3), 1:dims%n_var), intent(inout) :: f_qp_left, f_qp_right
!< primitive state at faces
type(celltype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: cells
!< Input cell quantities: volume
integer :: i, j, k, l
integer :: i_f=0, j_f=0, k_f=0
real(wp), dimension(3) :: P !< polynomial approximation
real(wp), dimension(3) :: B !< smoothness factor
real(wp), dimension(3) :: w !< wieght
real(wp), dimension(3) :: g !< linear wieght
real(wp), dimension(-2:2) :: u !<state_variable
real(wp) :: eps=1e-6
real(wp), dimension(-2:2) :: vol
real(wp) :: U11
real(wp) :: U00
real(wp) :: U21
real(wp) :: U10
real(wp) :: U01
real(wp) :: U12
real(wp) :: alpha12
real(wp) :: alpha01
real(wp) :: alpha10
real(wp) :: alpha21
g(1) = 1.0/10.0
g(2) = 6.0/10.0
g(3) = 3.0/10.0
i_f = flags(1)
j_f = flags(2)
k_f = flags(3)
do l = 1, dims%n_var
do k = 1-k_f, dims%kmx-1+k_f
do j = 1-j_f, dims%jmx-1+j_f
do i = 1-i_f, dims%imx-1+i_f
U(-2) = qp(i-2*i_f,j-2*j_f,k-2*k_f,l) !u_{i-2}
U(-1) = qp(i-1*i_f,j-1*j_f,k-1*k_f,l) !u_{i-1}
u( 0) = qp(i ,j ,k ,l) !u_{i}
U( 1) = qp(i+1*i_f,j+1*j_f,k+1*k_f,l) !u_{i+1}
U( 2) = qp(i+2*i_f,j+2*j_f,k+2*k_f,l) !u_{i+2}
Vol(-2) = cells(i-2*i_f,j-2*j_f,k-2*k_f)%volume !volume_{i-2}
Vol(-1) = cells(i-1*i_f,j-1*j_f,k-1*k_f)%volume !volume_{i-1}
Vol( 0) = cells(i ,j ,k )%volume !volume_{i}
Vol( 1) = cells(i+1*i_f,j+1*j_f,k+1*k_f)%volume !volume_{i+1}
Vol( 2) = cells(i+2*i_f,j+2*j_f,k+2*k_f)%volume !volume_{i+2}
alpha12 = Vol( 2)/(Vol( 1) + Vol( 2))
alpha01 = Vol( 1)/(Vol( 0) + Vol( 1))
alpha10 = Vol( 0)/(Vol(-1) + Vol( 0))
alpha21 = vol(-1)/(Vol(-2) + Vol(-1))
U01 = (1.0-alpha01)*U(0) + alpha01*U(1)
U12 = (1.0-alpha12)*U(1) + alpha12*U(2)
U10 = (1.0-alpha10)*U(-1) + alpha10*U(0)
U21 = (1.0-alpha21)*U(-2) + alpha21*U(-1)
U00 = U(-1) + (1.0-alpha21)*(U(-1) - U(-2))
U11 = U( 1) + alpha12*(U(1) - U(2))
P(1) = ( 6.0*U(0) - 1.0*U10 - 2.0*U00)/3.0
P(2) = (-1.0*U10 + 2.0*U(0) + 2.0*U01)/3.0
P(3) = ( 2.0*U01 + 2.0*U(1) - 1.0*U12)/3.0
B(1) = (13.0/12.0)*(2*U10-2.0*U00 )**2 + (1.0/4.0)*(4*U(0)-2.0*U10 -2.0*U00)**2
B(2) = (13.0/12.0)*(2*U10-4.0*U(0)+2*U01)**2 + (1.0/4.0)*(-2*u10 +2.0*U01)**2
B(3) = (13.0/12.0)*(2*U01-4.0*U(1)+2*U12)**2 + (1.0/4.0)*(-6*U01+8.0*U( 1)-2.0*U12)**2
w(:) = g(:)/(eps + B(:))**2
f_qp_left(i+i_f,j+j_f,k+k_f,l) = SUM(w*P)/SUM(w)
P(1) = ( 6.0*U(0) - 1.0*U01 - 2.0*U11)/3.0
P(2) = (-1.0*U01 + 2.0*U(0) + 2.0*U10)/3.0
P(3) = ( 2.0*U10 + 2.0*U(-1) - 1.0*U21)/3.0
B(1) = (13.0/12.0)*(2*U01-2.0*U11 )**2 + (1.0/4.0)*(4*U(0)-2.0*U01 -2.0*U11)**2
B(2) = (13.0/12.0)*(2*U01-4.0*U( 0)+2*U10)**2 + (1.0/4.0)*(-2*u01 +2.0*U10)**2
B(3) = (13.0/12.0)*(2*U10-4.0*U(-1)+2*U21)**2 + (1.0/4.0)*(-6*U10+8.0*U(-1)-2.0*U21)**2
w(:) = g(:)/(eps + B(:))**2
f_qp_right(i,j,k,l) = SUM(w*P)/SUM(w)
end do
end do
end do
end do
end subroutine compute_face_states