A generalized subroutine to calculate flux through the input-argument direction, :x,y, or z which corresponds to the I,J, or K direction respectively
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(inout), | dimension(:, :, :, :) | :: | Flux | Store fluxes throught the any(I,J,K) 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_left | primitve state variables at face |
|
| 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 | primitve state variables at face |
|
| type(facetype), | intent(in), | dimension(-2:dims%imx+2+flags(1),-2:dims%jmx+2+flags(2),-2:dims%kmx+2+flags(3)) | :: | faces | Face quantities: area and unit normal |
|
| integer, | intent(in), | dimension(3) | :: | flags | flags for direction switch |
|
| type(flowtype), | intent(in) | :: | flow | Information about fluid flow: freestream-speed, ref-viscosity,etc. |
||
| type(boundarytype), | intent(in) | :: | bc | boundary conditions and fixed values |
||
| type(extent), | intent(in) | :: | dims | Extent of the domain:imx,jmx,kmx |
subroutine compute_flux(Flux, f_qp_left, f_qp_right, faces, flags, flow, bc, dims)
!< A generalized subroutine to calculate
!< flux through the input-argument direction, :x,y, or z
!< which corresponds to the I,J, or K direction respectively
!------------------------------------------------------------
implicit none
integer, dimension(3), intent(in) :: flags
!< flags for direction switch
type(extent), intent(in) :: dims
!< Extent of the domain:imx,jmx,kmx
type(flowtype), intent(in) :: flow
!< Information about fluid flow: freestream-speed, ref-viscosity,etc.
type(boundarytype), intent(in) :: bc
!< boundary conditions and fixed values
real(wp), dimension(:, :, :, :), intent(inout) :: Flux
!< Store fluxes throught the any(I,J,K) faces
type(facetype), dimension(-2:dims%imx+2+flags(1),-2:dims%jmx+2+flags(2),-2:dims%kmx+2+flags(3)), intent(in) :: faces
!< Face quantities: area and unit normal
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
!< primitve state variables at face
integer :: i, j, k
integer :: i_f, j_f, k_f ! Flags to determine face direction
real(wp), dimension(1:dims%n_var) :: F_plus, F_minus
real(wp) :: M_perp_left, M_perp_right
real(wp) :: alpha_plus, alpha_minus
real(wp) :: beta_left, beta_right
real(wp) :: M_plus, M_minus
real(wp) :: D_plus, D_minus
real(wp) :: c_plus, c_minus
real(wp) :: scrD_plus, scrD_minus
real(wp) :: sound_speed_avg, face_normal_speeds
real(wp) :: temp_c
DebugCall('compute_flux')
i_f = flags(1)
j_f = flags(2)
k_f = flags(3)
do k = 1, dims%kmx - 1 + k_f
do j = 1, dims%jmx - 1 + j_f
do i = 1, dims%imx - 1 + i_f
sound_speed_avg = 0.5 * (sqrt(flow%gm * f_qp_left(i, j, k,5) / &
f_qp_left(i, j, k,1) ) + &
sqrt(flow%gm * f_qp_right(i, j, k,5) / &
f_qp_right(i, j, k,1) ) )
! Compute '+' direction quantities
face_normal_speeds = f_qp_left(i, j, k,2) * faces(i, j, k)%nx + &
f_qp_left(i, j, k,3) * faces(i, j, k)%ny + &
f_qp_left(i, j, k,4) * faces(i, j, k)%nz
M_perp_left = face_normal_speeds / sound_speed_avg
alpha_plus = 0.5 * (1.0 + sign(1.0, M_perp_left))
beta_left = -max(0, 1 - floor(abs(M_perp_left)))
M_plus = 0.25 * ((1. + M_perp_left) ** 2.)
D_plus = 0.25 * ((1. + M_perp_left) ** 2.) * (2. - M_perp_left)
c_plus = (alpha_plus * (1.0 + beta_left) * M_perp_left) - &
beta_left * M_plus
scrD_plus = (alpha_plus * (1. + beta_left)) - &
(beta_left * D_plus)
! Compute '-' direction quantities
face_normal_speeds = f_qp_right(i, j, k, 2) * faces(i, j, k)%nx + &
f_qp_right(i, j, k, 3) * faces(i, j, k)%ny + &
f_qp_right(i, j, k, 4) * faces(i, j, k)%nz
M_perp_right = face_normal_speeds / sound_speed_avg
alpha_minus = 0.5 * (1.0 - sign(1.0, M_perp_right))
beta_right = -max(0, 1 - floor(abs(M_perp_right)))
M_minus = -0.25 * ((1. - M_perp_right) ** 2.)
D_minus = 0.25 * ((1. - M_perp_right) ** 2.) * (2. + M_perp_right)
c_minus = (alpha_minus * (1.0 + beta_right) * M_perp_right) - &
beta_right * M_minus
scrD_minus = (alpha_minus * (1. + beta_right)) - &
(beta_right * D_minus)
! AUSM modification
temp_c = c_plus + c_minus
c_plus = max(0., temp_c)
c_minus = min(0., temp_c)
! F plus mass flux
F_plus(1) = f_qp_left(i, j, k,1) * sound_speed_avg * c_plus
! F minus mass flux
F_minus(1) = f_qp_right(i, j, k,1) * sound_speed_avg * c_minus
F_plus(1) = F_plus(1) *(i_f*bc%make_F_flux_zero(i) &
+ j_f*bc%make_G_flux_zero(j) &
+ k_f*bc%make_H_flux_zero(k))
F_minus(1) = F_minus(1)*(i_f*bc%make_F_flux_zero(i) &
+ j_f*bc%make_G_flux_zero(j) &
+ k_f*bc%make_H_flux_zero(k))
! Construct other fluxes in terms of the F mass flux
F_plus(2) = (F_plus(1) * f_qp_left(i, j, k,2)) + &
(scrD_plus * f_qp_left(i, j, k,5) * faces(i, j, k)%nx)
F_plus(3) = (F_plus(1) * f_qp_left(i, j, k,3)) + &
(scrD_plus * f_qp_left(i, j, k,5) * faces(i, j, k)%ny)
F_plus(4) = (F_plus(1) * f_qp_left(i, j, k,4)) + &
(scrD_plus * f_qp_left(i, j, k,5) * faces(i, j, k)%nz)
F_plus(5) = F_plus(1) * &
((0.5 * (f_qp_left(i, j, k,2) ** 2. + &
f_qp_left(i, j, k,3) ** 2. + &
f_qp_left(i, j, k,4) ** 2.)) + &
((flow%gm / (flow%gm - 1.)) * f_qp_left(i, j, k,5) / &
f_qp_left(i, j, k,1)))
! Construct other fluxes in terms of the F mass flux
F_minus(2) = (F_minus(1) * f_qp_right(i, j, k,2)) + &
(scrD_minus * f_qp_right(i, j, k,5) * faces(i, j, k)%nx)
F_minus(3) = (F_minus(1) * f_qp_right(i, j, k,3)) + &
(scrD_minus * f_qp_right(i, j, k,5) * faces(i, j, k)%ny)
F_minus(4) = (F_minus(1) * f_qp_right(i, j, k,4)) + &
(scrD_minus * f_qp_right(i, j, k,5) * faces(i, j, k)%nz)
F_minus(5) = F_minus(1) * &
((0.5 * (f_qp_right(i, j, k,2) ** 2. + &
f_qp_right(i, j, k,3) ** 2. + &
f_qp_right(i, j, k,4) ** 2.)) + &
((flow%gm / (flow%gm - 1.)) * f_qp_right(i, j, k,5) / &
f_qp_right(i, j, k,1)))
!turbulent fluxes
if(dims%n_var>5) then
F_plus(6:) = F_Plus(1) * f_qp_left(i,j,k,6:)
F_minus(6:) = F_minus(1) * f_qp_right(i,j,k,6:)
end if
! Multiply in the face areas
F_plus(:) = F_plus(:) * faces(i, j, k)%A
F_minus(:) = F_minus(:) * faces(i, j, k)%A
! Get the total flux for a face
Flux(i, j, k, :) = F_plus(:) + F_minus(:)
end do
end do
end do
end subroutine compute_flux