fluxes_kt.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: fluxes_kt.f90                                                     #
!#                                                                           #
!# Copyright (C) 2007-2018                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Jannes Klee <jklee@astrophysik.uni-kiel.de>                               #
!# Jubin Lirawi <jlirawi@astrophysik.uni-kiel.de>                            #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE fluxes_kt_mod
  USE fluxes_base_mod
  USE mesh_base_mod
  USE marray_compound_mod
!#ifdef PARALLEL
!       DEFAULT_MPI_REAL
!#endif
  USE physics_base_mod
  USE reconstruction_base_mod
  USE common_dict
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPIF_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  PRIVATE
  TYPE, EXTENDS (fluxes_base) :: fluxes_kt
  CONTAINS
    PROCEDURE         :: initfluxes_kt
    PROCEDURE         :: calculatefluxes
    PROCEDURE         :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       fluxes_kt
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initfluxes_kt(this,Mesh,Physics,config,IO)
    USE mesh_midpoint_mod, ONLY: mesh_midpoint
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_kt),    INTENT(INOUT)       :: this
    CLASS(mesh_base),    INTENT(IN)          :: Mesh
    CLASS(physics_base), INTENT(IN)          :: Physics
    TYPE(Dict_TYP),      INTENT(IN), POINTER :: config,IO
    !------------------------------------------------------------------------!
    CALL this%InitFluxes(mesh,physics,config,io,kt,"KT")

    SELECT TYPE(mesh)
    TYPE IS(mesh_midpoint)
      ! do nothing
!    CASE(TRAPEZOIDAL)
!      CALL this%Error("InitFluxes_kt", "Trapezoidal not supported from > 0.6")
    CLASS DEFAULT
      CALL this%Error("InitFluxes_kt", "Unknown mesh type.")
    END SELECT
  END SUBROUTINE initfluxes_kt


  PURE SUBROUTINE calculatefluxes(this,Mesh,Physics,pvar,cvar, &
                  xfluxdydz,yfluxdzdx,zfluxdxdy)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_kt),    INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)    :: mesh
    CLASS(physics_base), INTENT(INOUT) :: physics
    CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,Physics%vnum), &
                         INTENT(OUT)   :: xfluxdydz,yfluxdzdx,zfluxdxdy
    !------------------------------------------------------------------------!
    INTEGER :: i,j,k,l,m
    !------------------------------------------------------------------------!
    ! execute generic tasks common to all flux types
    CALL this%CalculateFaceData(mesh,physics,pvar,cvar)

    ! compute numerical fluxes along x-direction (west and east) devided by dydz
    IF (mesh%INUM.GT.1) THEN
      m = 1 ! index counting the number of spatial directions for which transport is enabled
      ! physical fluxes
      CALL physics%CalculateFluxesX(mesh,2*m-1,2*m,this%prim,this%cons,this%pfluxes)
!NEC$ SHORTLOOP
      DO l=1,physics%VNUM
!NEC$ IVDEP
        DO k=mesh%KGMIN,mesh%KGMAX
!NEC$ IVDEP
          DO j=mesh%JGMIN,mesh%JGMAX
!NEC$ IVDEP
            DO i=mesh%IMIN+mesh%IM1,mesh%IMAX
                xfluxdydz(i,j,k,l) = mesh%dAxdydz(i+1,j,k,1) / &
                        (this%maxwav%data4d(i,j,k,m) - this%minwav%data4d(i,j,k,m)) * &
                        (this%maxwav%data4d(i,j,k,m) * this%pfluxes%data5d(i,j,k,2*m,l) - &
                        this%minwav%data4d(i,j,k,m) * this%pfluxes%data5d(i+1,j,k,2*m-1,l) + &
                        this%minwav%data4d(i,j,k,m) * this%maxwav%data4d(i,j,k,m) * &
                    (this%cons%data5d(i+1,j,k,2*m-1,l) - this%cons%data5d(i,j,k,2*m,l)))
            END DO
          END DO
        END DO
      END DO
    ELSE
        m = 0
        xfluxdydz(:,:,:,:) = 0.0
    END IF

    ! compute numerical fluxes along y-direction (south and north) devided by dzdx
    IF (mesh%JNUM.GT.1) THEN
      m = m + 1 ! increase wave speed index, may be 1 or 2 now depending
                ! on whether there was transport in x-direction or not
      ! physical fluxes
      CALL physics%CalculateFluxesY(mesh,2*m-1,2*m,this%prim,this%cons,this%pfluxes)
!NEC$ SHORTLOOP
      DO l=1,physics%VNUM
!NEC$ IVDEP
        DO k=mesh%KGMIN,mesh%KGMAX
!NEC$ IVDEP
          DO j=mesh%JMIN+mesh%JM1,mesh%JMAX
!NEC$ IVDEP
            DO i=mesh%IGMIN,mesh%IGMAX
                yfluxdzdx(i,j,k,l) = mesh%dAydzdx(i,j+1,k,1) / &
                        (this%maxwav%data4d(i,j,k,m) - this%minwav%data4d(i,j,k,m)) * &
                        (this%maxwav%data4d(i,j,k,m) * this%pfluxes%data5d(i,j,k,2*m,l) - &
                        this%minwav%data4d(i,j,k,m) * this%pfluxes%data5d(i,j+1,k,2*m-1,l) + &
                        this%minwav%data4d(i,j,k,m) * this%maxwav%data4d(i,j,k,m) * &
                    (this%cons%data5d(i,j+1,k,2*m-1,l) - this%cons%data5d(i,j,k,2*m,l)))
            END DO
          END DO
        END DO
      END DO
    ELSE
        yfluxdzdx(:,:,:,:) = 0.0
    END IF

    ! compute numerical fluxes along z-direction (bottom and top) devided by dxdy
    IF (mesh%KNUM.GT.1) THEN
      m = m + 1 ! increase wave speed index, may be 1, 2 or 3 now depending
                ! on whether there was transport in x- and/or y-direction or not
        ! physical fluxes
      CALL physics%CalculateFluxesZ(mesh,2*m-1,2*m,this%prim,this%cons,this%pfluxes)
!NEC$ SHORTLOOP
      DO l=1,physics%VNUM
!NEC$ IVDEP
        DO k=mesh%KMIN+mesh%KM1,mesh%KMAX
!NEC$ IVDEP
          DO j=mesh%JGMIN,mesh%JGMAX
!NEC$ IVDEP
            DO i=mesh%IGMIN,mesh%IGMAX
                zfluxdxdy(i,j,k,l) = mesh%dAzdxdy(i,j,k+1,1) / &
                        (this%maxwav%data4d(i,j,k,m) - this%minwav%data4d(i,j,k,m)) * &
                        (this%maxwav%data4d(i,j,k,m) * this%pfluxes%data5d(i,j,k,2*m,l) - &
                        this%minwav%data4d(i,j,k,m) * this%pfluxes%data5d(i,j,k+1,2*m-1,l) + &
                        this%minwav%data4d(i,j,k,m) * this%maxwav%data4d(i,j,k,m) * &
                    (this%cons%data5d(i,j,k+1,2*m-1,l) - this%cons%data5d(i,j,k,2*m,l)))
            END DO
          END DO
        END DO
      END DO
    ELSE
        zfluxdxdy(:,:,:,:) = 0.0
    END IF
  END SUBROUTINE calculatefluxes

  ! Here was once a routine for BilinearInterpolation, which was needed for
  ! trapezoidal rule. See (<= 0.6) version of fosite.

  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_kt), INTENT(INOUT)  :: this
    !------------------------------------------------------------------------!
    CALL this%Finalize_base()
  END SUBROUTINE finalize

END MODULE fluxes_kt_mod