ausm.f90 Source File

Flux splitting scheme: AUSM


This file depends on

sourcefile~~ausm.f90~~EfferentGraph sourcefile~ausm.f90 ausm.f90 sourcefile~vartypes.f90 vartypes.f90 sourcefile~ausm.f90->sourcefile~vartypes.f90

Files dependent on this one

sourcefile~~ausm.f90~~AfferentGraph sourcefile~ausm.f90 ausm.f90 sourcefile~scheme.f90 scheme.f90 sourcefile~scheme.f90->sourcefile~ausm.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 ausm
    !<
    !< Reference: Liou, M.S. and Steffen Jr, C.J., A new flux splitting scheme, 
    !< Journal of Computational physics, vol. 107, no. 1, pp.23-39, 1993
    !-------------------------------------------------------------------
#include "../../../debug.h"
#include "../../../error.h"
     use vartypes
    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-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

        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 ausm