atl_linearEuler_numflux_module.f90 Source File

! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2016 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2017 Daniel Petró <daniel.petro@student.uni-siegen.de>
! Copyright (c) 2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2018 Harald Klimach <harald.klimach@uni-siegen.de>
!
! Permission to use, copy, modify, and distribute this software for any
! purpose with or without fee is hereby granted, provided that the above
! copyright notice and this permission notice appear in all copies.
!
! THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHORS DISCLAIM ALL WARRANTIES
! WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
! ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
! WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
! ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
! OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
! **************************************************************************** !

!> module that holds all routines to calculate the flux for
!! hyperbolic linearzied gas dynamic equations.

module atl_LinearEuler_numflux_module
  ! Treelm modules
  use env_module,         only: rk

  use atl_eqn_LinearEuler_module,      only: atl_LinearEuler_type

  implicit none

  private

  public :: atl_LinearEuler_numflux_subleft
  public :: atl_LinearEuler_numflux_subright
  public :: atl_LinearEuler_numflux_superleft
  public :: atl_LinearEuler_numflux_superright

contains


! ******************************************************************************
subroutine atl_LinearEuler_numflux_subright( nSides, nFaceDofs,  faceRep,      &
    &                                        faceFlux, leftPos, rightPos, var, &
    &                                        LinearEuler, iDir                 )
    ! --------------------------------------------------------------------------
    !> Datatype for LinearEuler equation include all background data
    type(atl_LinearEuler_type), intent(in) :: LinearEuler
    integer, intent(in) :: nFaceDofs, nSides
    real(kind=rk), intent(in) :: faceRep(:, :, :, :)
    real(kind=rk), intent(inout) :: faceFlux(:, :, :, :)
    integer, intent(in) :: leftPos(nSides), rightPos(nsides)
    integer, intent(in) :: var(:)
    !> Direction of the flow, used for background velocity
    integer, intent(in) :: idir
    ! --------------------------------------------------------------------------
    integer :: iSide, left, right, iDof, iter
    real(kind=rk) :: leftstate(5), rightstate(5)
    !> intermediated state on the interface ( in the * region)
    real(kind=rk) :: dens, velX, velY, velZ, p
    real(kind=rk) :: c, c_inv, c_cube_inv, rho0, rho0_inv, u0
    ! --------------------------------------------------------------------------

    c = LinearEuler%SpeedOfSound
    c_inv = 1 / c
    c_cube_inv = 1 / (c ** 2)
    rho0 = LinearEuler%density_0
    rho0_inv = 1 / rho0
    u0 = LinearEuler%velocity_0(iDir)

    ! loop over all dofs
    do iter=1,nFaceDofs*nSides
      iDof = (iter-1)/(nSides)+1
      iSide = mod(iter-1,nSides)+1

      ! The position of the left and right element in the state
      ! vector.
      left = leftPos(iSide)
      right = rightPos(iSide)

      leftstate(1) = faceRep(left, iDof, var(1), 2)
      leftstate(2) = faceRep(left, iDof, var(2), 2)
      leftstate(3) = faceRep(left, iDof, var(3), 2)
      leftstate(4) = faceRep(left, iDof, var(4), 2)
      leftstate(5) = faceRep(left, iDof, var(5), 2)
      rightstate(1) = faceRep(right, iDof, var(1), 1)
      rightstate(2) = faceRep(right, iDof, var(2), 1)
      rightstate(3) = faceRep(right, iDof, var(3), 1)
      rightstate(4) = faceRep(right, iDof, var(4), 1)
      rightstate(5) = faceRep(right, iDof, var(5), 1)

      ! density
      dens = 0.5 * c_cube_inv * ( rightstate(5) + leftstate(5)   &
        &    + c * rho0 * (leftstate(2) - rightstate(2))) &
        &    + rightstate(1) - rightstate(5) * c_cube_inv
      ! velocity
      velX = 0.5 * c_inv * rho0_inv * (leftstate(5) - rightstate(5) &
        &    + c * rho0 * (rightstate(2) + leftstate(2)))
      ! pressure
      p = 0.5 * (rightstate(5) + leftstate(5) &
        & + c * rho0 * (leftstate(2) - rightstate(2)))

      velY = rightstate(3)
      velZ = rightstate(4)

      ! 1....disturbance in density
      faceFlux(left, iDof, var(1), 2) = u0 * dens + rho0 * velX
      ! 2....disturbance velocity in x direction
      faceFlux(left, iDof, var(2), 2) = u0 * velX + p * rho0_inv
      ! 2....disturbance velocity in y direction
      faceFlux(left, iDof, var(3), 2) = u0 * velY
      ! 2....disturbance velocity in z direction
      faceFlux(left, iDof, var(4), 2) = u0 * velZ
      ! 5....distrubance in pressure
      faceFlux(left, iDof, var(5), 2) = u0 * p + LinearEuler%isen_coef &
        &                               * LinearEuler%pressure_0 * velX

      ! Assign the same flux for both adjacent elements
      faceFlux(right, iDof, var(1), 1) = faceFlux(left, iDof, var(1), 2)
      faceFlux(right, iDof, var(2), 1) = faceFlux(left, iDof, var(2), 2)
      faceFlux(right, iDof, var(3), 1) = faceFlux(left, iDof, var(3), 2)
      faceFlux(right, iDof, var(4), 1) = faceFlux(left, iDof, var(4), 2)
      faceFlux(right, iDof, var(5), 1) = faceFlux(left, iDof, var(5), 2)

    end do

end subroutine atl_LinearEuler_numFlux_subright
! ******************************************************************************


! ******************************************************************************
subroutine atl_LinearEuler_numflux_subleft( nSides, nFaceDofs, faceRep,        &
    &                                       faceFlux, leftPos,  rightPos, var, &
    &                                       LinearEuler, iDir                  )
    ! --------------------------------------------------------------------------
    !> Datatype for LinearEuler equation include all background data
    type(atl_LinearEuler_type), intent(in) :: LinearEuler
    integer, intent(in) :: nFaceDofs, nSides
    real(kind=rk), intent(in) :: faceRep(:, :, :, :)
    real(kind=rk), intent(inout) :: faceFlux(:, :, :, :)
    integer, intent(in) :: leftPos(nSides), rightPos(nsides)
    integer, intent(in) :: var(:)
    !> Direction of the flow, used for background velocity
    integer, intent(in) :: idir
    ! --------------------------------------------------------------------------
    integer :: iSide, left, right, iDof, iter
    real(kind=rk) :: leftstate(5), rightstate(5)
    !> intermediated state on the interface ( in the * region)
    real(kind=rk) :: dens, velX, velY, velZ, p
    real(kind=rk) :: c, c_inv, c_cube_inv, rho0, rho0_inv, u0
    ! --------------------------------------------------------------------------

    c = LinearEuler%SpeedOfSound
    c_inv = 1 / c
    c_cube_inv = 1 / (c ** 2)
    rho0 = LinearEuler%density_0
    rho0_inv = 1 / rho0
    u0 = LinearEuler%velocity_0(iDir)

    ! loop over all dofs
    do iter=1,nFaceDofs*nSides
      iDof = (iter-1)/(nSides)+1
      iSide = mod(iter-1,nSides)+1

      ! The position of the left and right element in the state
      ! vector.
      left = leftPos(iSide)
      right = rightPos(iSide)

      leftstate(1) = faceRep(left, iDof, var(1), 2)
      leftstate(2) = faceRep(left, iDof, var(2), 2)
      leftstate(3) = faceRep(left, iDof, var(3), 2)
      leftstate(4) = faceRep(left, iDof, var(4), 2)
      leftstate(5) = faceRep(left, iDof, var(5), 2)
      rightstate(1) = faceRep(right, iDof, var(1), 1)
      rightstate(2) = faceRep(right, iDof, var(2), 1)
      rightstate(3) = faceRep(right, iDof, var(3), 1)
      rightstate(4) = faceRep(right, iDof, var(4), 1)
      rightstate(5) = faceRep(right, iDof, var(5), 1)

      ! density
      dens = 0.5 * c_cube_inv * ( rightstate(5) + leftstate(5) &
        &    + c * rho0 * ( leftstate(2) - rightstate(2) ) )   &
        &    + leftstate(1) - leftstate(5) * c_cube_inv
      ! velocity
      velX = 0.5 * c_inv * rho0_inv * ( leftstate(5) - rightstate(5) &
        &    + c * rho0 * ( rightstate(2) + leftstate(2) ) )
      ! pressure
      p = 0.5 * ( rightstate(5) + leftstate(5) &
        & + c * rho0 * ( leftstate(2) - rightstate(2) ) )

      velY = leftstate(3)
      velZ = leftstate(4)

      ! 1....disturbance in density
      faceFlux(left, iDof, var(1), 2) = u0 * dens + rho0 * velX
      ! 2....disturbance velocity in x direction
      faceFlux(left, iDof, var(2), 2) = u0 * velX + p * rho0_inv
      ! 2....disturbance velocity in y direction
      faceFlux(left, iDof, var(3), 2) = u0 * velY
      ! 2....disturbance velocity in z direction
      faceFlux(left, iDof, var(4), 2) = u0 * velZ
      ! 5....distrubance in pressure
      faceFlux(left, iDof, var(5), 2) = u0 * p + LinearEuler%isen_coef &
        &                               * LinearEuler%pressure_0 * velX

      ! Assign the same flux for both adjacent elements
      faceFlux(right, iDof, var(1), 1) = faceFlux(left, iDof, var(1), 2)
      faceFlux(right, iDof, var(2), 1) = faceFlux(left, iDof, var(2), 2)
      faceFlux(right, iDof, var(3), 1) = faceFlux(left, iDof, var(3), 2)
      faceFlux(right, iDof, var(4), 1) = faceFlux(left, iDof, var(4), 2)
      faceFlux(right, iDof, var(5), 1) = faceFlux(left, iDof, var(5), 2)

    end do

end subroutine atl_LinearEuler_numFlux_subleft
! ******************************************************************************


! ******************************************************************************
subroutine atl_LinearEuler_numflux_superright( nSides, nFaceDofs, faceRep,  &
    &                                          faceFlux, leftPos, rightPos, &
    &                                          var, LinearEuler, iDir       )
    ! --------------------------------------------------------------------------
    !> Datatype for LinearEuler equation include all background data
    type(atl_LinearEuler_type), intent(in) :: LinearEuler
    integer, intent(in) :: nFaceDofs, nSides
    real(kind=rk), intent(in) :: faceRep(:, :, :, :)
    real(kind=rk), intent(inout) :: faceFlux(:, :, :, :)
    integer, intent(in) :: leftPos(nSides), rightPos(nsides)
    integer, intent(in) :: var(:)
    !> Direction of the flow, used for background velocity
    integer, intent(in) :: idir
    ! --------------------------------------------------------------------------
    integer :: iSide, right, left, iDof, iter
    !> intermediated state on the interface ( in the * region)
    real(kind=rk) :: dens, velX, velY, velZ, p
    real(kind=rk) :: rho0, rho0_inv, u0
    ! --------------------------------------------------------------------------

    rho0 = LinearEuler%density_0
    rho0_inv = 1 / rho0
    u0 = LinearEuler%velocity_0(iDir)

    ! loop over all dofs
    do iter=1,nFaceDofs*nSides
      iDof = (iter-1)/(nSides)+1
      iSide = mod(iter-1,nSides)+1

      ! The position of the left and right element in the state
      ! vector.
      right = rightPos(iSide)
      left = leftPos(iSide)

      ! this is the case when all characteristics are going to the left
      ! hence the state is always transport from right to left
      ! which implies the state becomes the state from the right
      dens = faceRep(right, iDof, var(1), 1)
      velX = faceRep(right, iDof, var(2), 1)
      velY = faceRep(right, iDof, var(3), 1)
      velZ = faceRep(right, iDof, var(4), 1)
      p = faceRep(right, iDof, var(5), 1)

      ! 1....disturbance in density
      faceFlux(left, iDof, var(1), 2) = u0 * dens + rho0 * velX
      ! 2....disturbance velocity in x direction
      faceFlux(left, iDof, var(2), 2) = u0 * velX + p * rho0_inv
      ! 2....disturbance velocity in y direction
      faceFlux(left, iDof, var(3), 2) = u0 * velY
      ! 2....disturbance velocity in z direction
      faceFlux(left, iDof, var(4), 2) = u0 * velZ
      ! 5....distrubance in pressure
      faceFlux(left, iDof, var(5), 2) = u0 * p + LinearEuler%isen_coef &
        &                               * LinearEuler%pressure_0 * velX

      ! Assign the same flux for both adjacent elements
      faceFlux(right, iDof, var(1), 1) = faceFlux(left, iDof, var(1), 2)
      faceFlux(right, iDof, var(2), 1) = faceFlux(left, iDof, var(2), 2)
      faceFlux(right, iDof, var(3), 1) = faceFlux(left, iDof, var(3), 2)
      faceFlux(right, iDof, var(4), 1) = faceFlux(left, iDof, var(4), 2)
      faceFlux(right, iDof, var(5), 1) = faceFlux(left, iDof, var(5), 2)

    end do

end subroutine atl_LinearEuler_numFlux_superright
! ******************************************************************************


! ******************************************************************************
subroutine atl_LinearEuler_numflux_superleft( nSides, nFaceDofs, faceRep,  &
    &                                         faceFlux, leftPos, rightPos, &
    &                                         var, LinearEuler, iDir       )
    ! --------------------------------------------------------------------------
    !> Datatype for LinearEuler equation include all background data
    type(atl_LinearEuler_type), intent(in) :: LinearEuler
    integer, intent(in) :: nFaceDofs, nSides
    real(kind=rk), intent(in) :: faceRep(:, :, :, :)
    real(kind=rk), intent(inout) :: faceFlux(:, :, :, :)
    integer, intent(in) :: leftPos(nSides), rightPos(nsides)
    integer, intent(in) :: var(:)
    !> Direction of the flow, used for background velocity
    integer, intent(in) :: idir
    ! --------------------------------------------------------------------------
    integer :: iSide, left, right, iDof, iter
    !> intermediated state on the interface ( in the * region)
    real(kind=rk) :: dens, velX, velY, velZ, p
    real(kind=rk) :: rho0, rho0_inv, u0
    ! --------------------------------------------------------------------------

    rho0 = LinearEuler%density_0
    rho0_inv = 1 / rho0
    u0 = LinearEuler%velocity_0(iDir)

    ! loop over all dofs
    do iter=1,nFaceDofs*nSides
      iDof = (iter-1)/(nSides)+1
      iSide = mod(iter-1,nSides)+1

      ! The position of the left and right element in the state
      ! vector.
      right = rightPos(iSide)
      left = leftPos(iSide)

      ! this is the case when all characteristics are going to the right
      ! hence the state is always transport from left to right
      ! which implies the state becomes the state from the left
      dens = faceRep(left, iDof, var(1), 2)
      velX = faceRep(left, iDof, var(2), 2)
      velY = faceRep(left, iDof, var(3), 2)
      velZ = faceRep(left, iDof, var(4), 2)
      p = faceRep(left, iDof, var(5), 2)

      ! 1....disturbance in density
      faceFlux(left, iDof, var(1), 2) = u0 * dens + rho0 * velX
      ! 2....disturbance velocity in x direction
      faceFlux(left, iDof, var(2), 2) = u0 * velX + p * rho0_inv
      ! 2....disturbance velocity in y direction
      faceFlux(left, iDof, var(3), 2) = u0 * velY
      ! 2....disturbance velocity in z direction
      faceFlux(left, iDof, var(4), 2) = u0 * velZ
      ! 5....distrubance in pressure
      faceFlux(left, iDof, var(5), 2) = u0 * p + LinearEuler%isen_coef &
        &                               * LinearEuler%pressure_0 * velX

      ! Assign the same flux for both adjacent elements
      faceFlux(right, iDof, var(1), 1) = faceFlux(left, iDof, var(1), 2)
      faceFlux(right, iDof, var(2), 1) = faceFlux(left, iDof, var(2), 2)
      faceFlux(right, iDof, var(3), 1) = faceFlux(left, iDof, var(3), 2)
      faceFlux(right, iDof, var(4), 1) = faceFlux(left, iDof, var(4), 2)
      faceFlux(right, iDof, var(5), 1) = faceFlux(left, iDof, var(5), 2)

    end do

end subroutine atl_LinearEuler_numFlux_superleft
! ******************************************************************************

end module atl_LinearEuler_numflux_module