ausmP.f90 Source File

Flux splitting scheme: AUSM+


This file depends on

sourcefile~~ausmp.f90~~EfferentGraph sourcefile~ausmp.f90 ausmP.f90 sourcefile~face_interpolant.f90 face_interpolant.f90 sourcefile~ausmp.f90->sourcefile~face_interpolant.f90 sourcefile~vartypes.f90 vartypes.f90 sourcefile~ausmp.f90->sourcefile~vartypes.f90 sourcefile~face_interpolant.f90->sourcefile~vartypes.f90 sourcefile~weno.f90 weno.f90 sourcefile~face_interpolant.f90->sourcefile~weno.f90 sourcefile~weno_nm.f90 weno_NM.f90 sourcefile~face_interpolant.f90->sourcefile~weno_nm.f90 sourcefile~utils.f90 utils.f90 sourcefile~face_interpolant.f90->sourcefile~utils.f90 sourcefile~ppm.f90 ppm.f90 sourcefile~face_interpolant.f90->sourcefile~ppm.f90 sourcefile~muscl.f90 muscl.f90 sourcefile~face_interpolant.f90->sourcefile~muscl.f90 sourcefile~weno.f90->sourcefile~vartypes.f90 sourcefile~weno_nm.f90->sourcefile~vartypes.f90 sourcefile~ppm.f90->sourcefile~vartypes.f90 sourcefile~muscl.f90->sourcefile~vartypes.f90

Files dependent on this one

sourcefile~~ausmp.f90~~AfferentGraph sourcefile~ausmp.f90 ausmP.f90 sourcefile~scheme.f90 scheme.f90 sourcefile~scheme.f90->sourcefile~ausmp.f90 sourcefile~solver.f90 solver.f90 sourcefile~solver.f90->sourcefile~scheme.f90 sourcefile~update.f90 update.f90 sourcefile~solver.f90->sourcefile~update.f90 sourcefile~update.f90->sourcefile~scheme.f90 sourcefile~main.f90 main.f90 sourcefile~main.f90->sourcefile~solver.f90

Contents

Source Code


Source Code

    !< Flux splitting scheme: AUSM+
module ausmP
    !< Flux splitting scheme: AUSM+ module ausmP !< 
    !< Reference: Liou, M. S., “A sequel to AUSM: AUSM+,” 
    !< Journal of Computational Physics, vol. 129, pp. 364–382, 1996 
    !------------------------------------------------------------------- 
#include "../../../debug.h"
#include "../../../error.h"    
    use vartypes
    use face_interpolant, only: x_qp_left, x_qp_right 
    use face_interpolant, only: y_qp_left, y_qp_right
    use face_interpolant, only:  z_qp_left, z_qp_right


    implicit none
    private

    public :: compute_fluxes
    
    contains

        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


        subroutine compute_fluxes(F,G,H, x_qp_l, x_qp_r, y_qp_l, y_qp_r, z_qp_l, z_qp_r, Ifaces, Jfaces, Kfaces, flow, bc, dims)
        !subroutine compute_fluxes(F,G,H, flow, dims)
          !< Call to compute fluxes throught faces in each direction
            
            implicit none
            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.
            real(wp), dimension(:, :, :, :), intent(inout) :: F
            !< Store fluxes throught the I faces
            real(wp), dimension(:, :, :, :), intent(inout) :: G
            !< Store fluxes throught the J faces
            real(wp), dimension(:, :, :, :), intent(inout) :: H
            !< Store fluxes throught the K faces
            type(facetype), dimension(-2:dims%imx+3,-2:dims%jmx+2,-2:dims%kmx+2), intent(in) :: Ifaces
            !< Store face quantites for I faces 
            type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+3,-2:dims%kmx+2), intent(in) :: Jfaces
            !< Store face quantites for J faces 
            type(facetype), dimension(-2:dims%imx+2,-2:dims%jmx+2,-2:dims%kmx+3), intent(in) :: Kfaces
            !< Store face quantites for K faces 
            real(wp), dimension(0:dims%imx+1,1:dims%jmx-1,1:dims%kmx-1,1:dims%n_var), intent(inout) :: x_qp_l, x_qp_r
            !< Store primitive state at the I-face 
            real(wp), dimension(1:dims%imx-1,0:dims%jmx+1,1:dims%kmx-1,1:dims%n_var), intent(inout) :: y_qp_l, y_qp_r
            !< Store primitive state at the J-face 
            real(wp), dimension(1:dims%imx-1,1:dims%jmx-1,0:dims%kmx+1,1:dims%n_var), intent(inout) :: z_qp_l, z_qp_r
            !< Store primitive state at the K-face 
            type(boundarytype), intent(in) :: bc
            !< boundary conditions and fixed values
            integer, dimension(3) :: flags

            
            DebugCall('compute_fluxes')

            flags=(/1,0,0/)
            call compute_flux(F, x_qp_l, x_qp_r, Ifaces, flags, flow, bc, dims)
            if (any(isnan(F))) then
              Fatal_error
            end if    

            flags=(/0,1,0/)
            call compute_flux(G, y_qp_l, y_qp_r, Jfaces, flags, flow, bc, dims)
            if (any(isnan(G))) then 
              Fatal_error
            end if    
            
            if(dims%kmx==2) then
              H = 0.
            else
              flags=(/0,0,1/)
              call compute_flux(H, z_qp_l, z_qp_r, Kfaces, flags, flow, bc, dims)
            end if
            if (any(isnan(H))) then
              Fatal_error
            end if

        end subroutine compute_fluxes

end module ausmP