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 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~4~~CallsGraph proc~compute_flux~4 compute_flux debugcall debugcall proc~compute_flux~4->debugcall

Called by

proc~~compute_flux~4~~CalledByGraph proc~compute_flux~4 compute_flux proc~compute_fluxes~4 compute_fluxes proc~compute_fluxes~4->proc~compute_flux~4

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 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 for DO loop
            integer :: i_f, j_f, k_f 
            !< Flags to determine face direction
            real(wp), dimension(1:dims%n_var) :: F_plus
            !< Right flux through the face
            real(wp), dimension(1:dims%n_var) ::F_minus
            !< Left flux through  the face
            real(wp) :: pbar
            real(wp) :: mass
            real(wp) :: HL, HR
            !< Enthalpy
            real(wp) :: uL, uR
            !< X-component of velocity
            real(wp) :: vL, vR
            !< Y-component of velocity
            real(wp) :: wL, wR
            !< Z-component of velocity
            real(wp) :: pL, pR
            !< Pressure
            real(wp) :: rL, rR
            !< Density
            real(wp) :: cL, cR
            !< Speed sound left/right
            real(wp) :: C
            !< Speed of sound at face
            real(wp) :: ML, MR
            !< Mach number left/right
            real(wp) :: VnL, VnR
            !< Face normal velocity left/right
            real(wp) :: betaL, betaR
            real(wp) :: alphaL, alphaR
            real(wp) :: FmL, FmR
            real(wp) :: Mface
            real(wp) :: Cs


            DebugCall('compute_flux '//trim(f_dir))
            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

                ! -- primitve face state assignment --
                ! ---- left face quantities ----
                rL = f_qp_left(i,j,k,1)
                uL = f_qp_left(i,j,k,2)
                vL = f_qp_left(i,j,k,3)
                wL = f_qp_left(i,j,k,4)
                pL = f_qp_left(i,j,k,5)

                ! ---- right face quantities ----
                rR = f_qp_right(i,j,k,1)
                uR = f_qp_right(i,j,k,2)
                vR = f_qp_right(i,j,k,3)
                wR = f_qp_right(i,j,k,4)
                pR = f_qp_right(i,j,k,5)

                !-- calculated quntaties --
                ! ---- total enthalpy ----
                HL = (0.5*(uL*uL + vL*vL + wL*wL)) + ((flow%gm/(flow%gm - 1.))*pL/rL)
                HR = (0.5*(uR*uR + vR*vR + wR*wR)) + ((flow%gm/(flow%gm - 1.))*pR/rR)

                ! ---- face normal velocity ----
                VnL = uL*faces(i, j, k)%nx + vL*faces(i, j, k)%ny + wL*faces(i, j, k)%nz
                VnR = uR*faces(i, j, k)%nx + vR*faces(i, j, k)%ny + wR*faces(i, j, k)%nz

                ! ---- speed of sound ----
                cs = sqrt(2.0*(flow%gm-1.0)*(0.5*(HL + HR))/(flow%gm+1.0))
                cL = cs*cs/(max(cs, abs(VnL)))
                cR = cs*cs/(max(cs, abs(VnR)))
                C  = min(cL, CR)

                ! ---- Mach at face ----
                ML = VnL/C
                MR = VnR/C

                ! ---- switch for supersonic flow ----
                alphaL= max(0, 1-floor(abs(ML)))
                alphaR= max(0, 1-floor(abs(MR)))

                
                ! Compute '+' direction quantities
                FmL   = (0.5*(1.0+sign(1.0,ML))*(1.0-alphaL)*ML) + alphaL*0.25*((1.0+ML)**2)
                betaL = (0.5*(1.0+sign(1.0,ML))*(1.0-alphaL))    + alphaL*0.25*((1.0+ML)**2) * (2.0 - ML)

                ! Compute '-' direction quantities
                FmR   = (0.5*(1.0-sign(1.0,MR))*(1.0-alphaR)*MR) - alphaR*0.25*((1.0-MR)**2)
                betaR = (0.5*(1.0-sign(1.0,MR))*(1.0-alphaR))    + alphaR*0.25*((1.0-MR)**2)*(2.0 + MR)

                !AUSM+modification
                ! Compute '+' direction quantities
                FmL   = FmL   + alphaL*0.1250*((ML**2-1.0)**2)
                betaL = betaL + alphaL*0.1875*((ML**2-1.0)**2)*ML

                ! Compute '-' direction quantities
                FmR   = FmR   - alphaR*0.1250*((MR**2-1.0)**2)
                betaR = betaR - alphaR*0.1875*((MR**2-1.0)**2)*MR

                
                ! mass coefficient
                Mface = FmL + FmR
                ! -- Pressure coeffient--
                pbar = betaL*pL + betaR*pR


                ! -- mass --
                if(Mface>0.0)then
                    mass = Mface*C*rL
                else
                    mass = Mface*c*rR
                end if

                mass = mass *(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 plus mass flux
                ! Construct other fluxes in terms of the F mass flux
                F_plus(1) = 0.5*(mass + abs(mass))
                F_plus(2) = (F_plus(1) * uL)
                F_plus(3) = (F_plus(1) * vL)
                F_plus(4) = (F_plus(1) * wL)
                F_plus(5) = (F_plus(1) * HL)
                
                ! F minus mass flux
                ! Construct other fluxes in terms of the F mass flux
                F_minus(1) = 0.5*(mass - abs(mass))
                F_minus(2) = (F_minus(1) * uR)
                F_minus(3) = (F_minus(1) * vR)
                F_minus(4) = (F_minus(1) * wR)
                F_minus(5) = (F_minus(1) * HR)

                !! -- Turbulence variables mass flux --
                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

                ! total flux
                Flux(i, j, k, :) = F_plus(:) + F_minus(:)

                ! Get the total flux for a face
                ! -- Pressure flux addition --
                Flux(i, j, K, 2) = Flux(i, j, k, 2) + (pbar * faces(i, j, k)%nx)
                Flux(i, j, K, 3) = Flux(i, j, k, 3) + (pbar * faces(i, j, k)%ny)
                Flux(i, j, K, 4) = Flux(i, j, k, 4) + (pbar * faces(i, j, k)%nz)

                Flux(i, j, k, :) = Flux(i, j, k, :)*faces(i,j,k)%A
              end do
             end do
            end do 

        end subroutine compute_flux