compute_flux Subroutine

private 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

Arguments

Type IntentOptional AttributesName
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 variable 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 variable 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


Calls

proc~~compute_flux~6~~CallsGraph proc~compute_flux~6 compute_flux debugcall debugcall proc~compute_flux~6->debugcall

Called by

proc~~compute_flux~6~~CalledByGraph proc~compute_flux~6 compute_flux proc~compute_fluxes~6 compute_fluxes proc~compute_fluxes~6->proc~compute_flux~6

Contents

Source Code


Source Code

        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 variable 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) :: M_ldfss, M_plus_ldfss, M_minus_ldfss

            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)
                
                ! LDFSS0 modification             
                M_ldfss = 0.25 * beta_left * beta_right * &
                    (sqrt((M_perp_left ** 2 + M_perp_right ** 2) * 0.5) &
                     - 1) ** 2
                M_plus_ldfss = M_ldfss * &
                        (1 - (f_qp_left(i, j, k,5) - f_qp_right(i, j, k,5)) / &
                            (2 * f_qp_left(i, j, k,1) * (sound_speed_avg ** 2)))
                M_minus_ldfss = M_ldfss * &
                        (1 - (f_qp_left(i, j, k,5) - f_qp_right(i, j, k,5)) / &
                            (2 * f_qp_right(i, j, k,1) * (sound_speed_avg ** 2)))
                c_plus = c_plus - M_plus_ldfss
                c_minus = c_minus + M_minus_ldfss

                ! 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