fluxes_base.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: fluxes_base.f90                                                   #
!#                                                                           #
!# Copyright (C) 2007-2017                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Jannes Klee <jklee@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_base_mod
  USE logging_base_mod
  USE mesh_base_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE reconstruction_generic_mod
  USE physics_base_mod
  USE physics_generic_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, ABSTRACT, EXTENDS (logging_base) ::  fluxes_base
     CLASS(reconstruction_base), ALLOCATABLE &
                                  :: reconstruction
     CLASS(marray_base), ALLOCATABLE &
                                  :: minwav,maxwav
     CLASS(marray_compound), POINTER :: prim, &      !< primitive/conservative
                                        cons, &      !< state vectors on cell faces
                                        pfluxes
     REAL, DIMENSION(:,:,:,:), POINTER, CONTIGUOUS &
                                  :: bxflux,byflux, &
                                     bzflux,bxfold, &
                                     byfold,bzfold
  CONTAINS
    PROCEDURE                               :: initfluxes
    procedure(calculatefluxes), DEFERRED   :: calculatefluxes
    PROCEDURE                               :: calculatefacedata
    PROCEDURE                               :: getboundaryflux
    procedure(finalize), DEFERRED          :: finalize
    PROCEDURE                               :: finalize_base
  END TYPE fluxes_base

  abstract INTERFACE
    SUBROUTINE calculatefluxes(this,Mesh,Physics,pvar,cvar, &
                xfluxdydz,yfluxdzdx,zfluxdxdy)
      IMPORT fluxes_base,mesh_base,physics_base,marray_compound
      IMPLICIT NONE
      CLASS(fluxes_base),   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
    END SUBROUTINE
    SUBROUTINE finalize(this)
      IMPORT fluxes_base
      IMPLICIT NONE
      CLASS(fluxes_base), INTENT(INOUT)   :: this
    END SUBROUTINE
  END INTERFACE
  INTEGER, PARAMETER :: kt       = 1
!  INTEGER, PARAMETER :: HLL      = 2
!  INTEGER, PARAMETER :: HLLC     = 3
!  INTEGER, PARAMETER :: EXACT    = 4
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       fluxes_base, &
       ! constants
       kt!, HLL, HLLC, EXACT
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initfluxes(this,Mesh,Physics,config,IO,ftype,fname)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_base),  INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)    :: Mesh
    CLASS(physics_base), INTENT(IN)    :: Physics
    TYPE(Dict_TYP),POINTER             :: config,IO
    INTEGER                            :: err
    INTEGER                            :: ftype
    CHARACTER(LEN=*)                   :: fname
    !------------------------------------------------------------------------!
    TYPE(Dict_TYP),POINTER             :: IOrec => null()
    INTEGER                            :: valwrite,i
    CHARACTER(LEN=60)                  :: key
    !------------------------------------------------------------------------!
    CALL this%InitLogging(ftype,fname)

    ! check initialization of Mesh and Physics
    IF (.NOT.mesh%Initialized().OR..NOT.physics%Initialized()) &
         CALL this%Error("InitFluxes","mesh and/or physics module uninitialized")

    ! allocate memory for all arrays used in fluxes
    ALLOCATE(this%minwav,this%maxwav, &
      this%bxflux(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,2,physics%VNUM),      &
      this%byflux(mesh%KGMIN:mesh%KGMAX,mesh%IGMIN:mesh%IGMAX,2,physics%VNUM),      &
      this%bzflux(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,2,physics%VNUM),      &
      this%bxfold(mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,2,physics%VNUM),      &
      this%byfold(mesh%KGMIN:mesh%KGMAX,mesh%IGMIN:mesh%IGMAX,2,physics%VNUM),      &
      this%bzfold(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,2,physics%VNUM),      &
      stat = err)
    IF (err.NE.0) THEN
       CALL this%Error("InitFluxes", "Unable to allocate memory.")
    END IF

    ! create RANK 1 mesh arrays
    this%minwav = marray_base(mesh%NDIMS)
    this%maxwav = marray_base(mesh%NDIMS)

    ! initialize state vectors on cell interfaces
    CALL physics%new_statevector(this%prim,primitive,mesh%NFACES)
    CALL physics%new_statevector(this%cons,conservative,mesh%NFACES)

    ! initialize state vector for physical fluxes
    CALL physics%new_statevector(this%pfluxes,conservative,mesh%NFACES)

    ! print some information
    CALL this%Info(" FLUXES---> fluxes type        " // trim(this%GetName()))

    ! enable output of reconstructed cell face data if requested
    CALL getattr(config, "output/pstates", valwrite,0)
    IF(valwrite.EQ.1) THEN
      DO i=1, physics%VNUM
        key = trim(physics%pvarname(i)) // "_pstates"
        CALL setattr(io, trim(key), &
                     this%prim%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:,i))
      END DO
    END IF
    CALL getattr(config, "output/cstates", valwrite, 0)
    IF(valwrite.EQ.1) THEN
      DO i=1, physics%VNUM
        key = trim(physics%cvarname(i)) // "_cstates"
        CALL setattr(io, trim(key), &
                     this%cons%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:,i))
      END DO
    END IF
    CALL getattr(config, "output/pfluxes", valwrite, 0)
    IF(valwrite.EQ.1) THEN
      DO i=1, physics%VNUM
        key = trim(physics%pvarname(i)) // "_pfluxes"
        CALL setattr(io, trim(key), &
                     this%pfluxes%data5d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,:,i))
      END DO
    END IF

    CALL getattr(config, "output/wave_speeds", valwrite, 0)
    IF(valwrite.EQ.1) THEN
      CALL setattr(io, "minwav", &
                   this%minwav%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1:mesh%NDIMS))
      CALL setattr(io, "maxwav", &
                   this%maxwav%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1:mesh%NDIMS))
      END IF

    ! initialize reconstruction modules
    CALL new_reconstruction(this%Reconstruction,mesh,physics,config,iorec)

    ! add output data of reconstruction modules to IO dictionary
    IF(ASSOCIATED(iorec)) CALL setattr(io,"reconstruction",iorec)


    ! initialize boundary fluxes
    this%bxflux(:,:,:,:) = 0.
    this%byflux(:,:,:,:) = 0.
    this%bzflux(:,:,:,:) = 0.
    this%bxfold(:,:,:,:) = 0.
    this%byfold(:,:,:,:) = 0.
    this%bzfold(:,:,:,:) = 0.
  END SUBROUTINE initfluxes

  FUNCTION getboundaryflux(this,Mesh,Physics,direction,comm) RESULT(bflux)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_base),  INTENT(IN) :: this
    CLASS(mesh_base),    INTENT(IN) :: mesh
    CLASS(physics_base), INTENT(IN) :: physics
    INTEGER                         :: direction
    REAL, DIMENSION(Physics%VNUM)   :: bflux
    INTEGER, OPTIONAL               :: comm        ! communicator for MPI
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    REAL, DIMENSION(Physics%VNUM) :: bflux_local
    INTEGER                       :: sender_rank(1),dest_ranks(1),rank0(1)
    INTEGER                       :: dest_comm,union_comm
    INTEGER                       :: world_group,dest_group,union_group,sender_group
    INTEGER                       :: ierror
#endif
    !------------------------------------------------------------------------!
    INTENT(IN)                    :: direction
    !------------------------------------------------------------------------!
#ifdef PARALLEL
    IF (PRESENT(comm)) THEN
       dest_comm = comm
    ELSE
       dest_comm = mpi_comm_null
    END IF
#endif
    SELECT CASE(direction)
    CASE(west) ! western boundary flux
#ifdef PARALLEL
       IF (mesh%mycoords(1).EQ.0) THEN
          bflux_local(:) = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1,:),1),1)
       ELSE
          bflux_local(:) = 0.0
       END IF
#else
       bflux = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1,:),1),1)
#endif
    CASE(east) ! eastern boundary flux
#ifdef PARALLEL
       IF (mesh%mycoords(1).EQ.mesh%dims(1)-1) THEN
          bflux_local(:) = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2,:),1),1)
       ELSE
          bflux_local(:) = 0.0
       END IF
#else
       bflux = sum(sum(this%bxflux(mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,2,:),1),1)
#endif
    CASE(south) ! southern boundary flux
#ifdef PARALLEL
       IF (mesh%mycoords(2).EQ.0) THEN
          bflux_local(:) = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,1,:),1),1)
       ELSE
          bflux_local(:) = 0.0
       END IF
#else
       bflux = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,1,:),1),1)
#endif
    CASE (north) ! northern boundary flux
#ifdef PARALLEL
       IF (mesh%mycoords(2).EQ.mesh%dims(2)-1) THEN
          bflux_local(:) = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,2,:),1),1)
       ELSE
          bflux_local(:) = 0.0
       END IF
#else
       bflux = sum(sum(this%byflux(mesh%KMIN:mesh%KMAX,mesh%IMIN:mesh%IMAX,2,:),1),1)
#endif
    CASE (bottom) ! bottom boundary flux
#ifdef PARALLEL
       IF (mesh%mycoords(3).EQ.0) THEN
          bflux_local(:) = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1,:),1),1)
       ELSE
          bflux_local(:) = 0.0
       END IF
#else
       bflux = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,1,:),1),1)
#endif
    CASE (top) ! topper boundary flux
#ifdef PARALLEL
       IF (mesh%mycoords(3).EQ.mesh%dims(3)-1) THEN
          bflux_local(:) = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2,:),1),1)
       ELSE
          bflux_local(:) = 0.0
       END IF
#else
       bflux = sum(sum(this%bzflux(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,2,:),1),1)
#endif
    CASE DEFAULT
       CALL this%Error("GetBoundaryFlux","wrong direction")
    END SELECT

#ifdef PARALLEL
    ! if dest_comm is the world comm, use simpler and faster AllReduce
    IF(dest_comm.EQ.mpi_comm_world) THEN
      CALL mpi_allreduce(bflux_local,bflux,physics%VNUM,default_mpi_real, &
           mpi_sum,mpi_comm_world,ierror)
    ELSE
      ! create world group
      CALL mpi_comm_group(mpi_comm_world,world_group,ierror)
      ! determine the destination for the result
      IF(dest_comm.NE.mpi_comm_null) THEN
        dest_comm = comm
        ! set destination group
        CALL mpi_comm_group(dest_comm,dest_group,ierror)
      ELSE
        ! if no communicator is given, send the result
        ! to the rank0 process
        dest_ranks(1) = 0
        CALL mpi_group_incl(world_group,1,dest_ranks,dest_group,ierror)
      END IF
      ! collect and sum up the result in process with rank 0 with respect to the
      ! subset of MPI processes at the boundary under consideration
      IF (mesh%comm_boundaries(direction).NE.mpi_comm_null) THEN
         CALL mpi_reduce(bflux_local,bflux,physics%VNUM,default_mpi_real, &
              mpi_sum,0,mesh%comm_boundaries(direction),ierror)
      ELSE
         bflux(:) = 0.0
      END IF
      ! get sender group
      rank0(1) = mesh%rank0_boundaries(direction)
      CALL mpi_group_incl(world_group,1,rank0,sender_group,ierror)
      ! merge sender with destination
      CALL mpi_group_union(sender_group,dest_group,union_group,ierror)
      ! create a communicator for the union group
      CALL mpi_comm_create(mpi_comm_world,union_group,union_comm,ierror)
      IF (union_comm.NE.mpi_comm_null) THEN
         ! get rank of sender in union group
         rank0(1) = 0
         CALL mpi_group_translate_ranks(sender_group,1,rank0,union_group,sender_rank,ierror)
         IF (sender_rank(1).EQ.mpi_undefined) &
             CALL this%Error("GetBoundaryFlux","sender rank undefined")
         ! send result to all processes in communicator 'union_comm'
         CALL mpi_bcast(bflux,physics%VNUM,default_mpi_real,sender_rank(1),union_comm,ierror)
         ! free union communicator
         CALL mpi_comm_free(union_comm,ierror)
      END IF
      ! free all groups
      CALL mpi_group_free(union_group,ierror)
      CALL mpi_group_free(sender_group,ierror)
      CALL mpi_group_free(world_group,ierror)
      CALL mpi_group_free(dest_group,ierror)
    END IF

#endif
  END FUNCTION getboundaryflux

  PURE SUBROUTINE calculatefacedata(this,Mesh,Physics,pvar,cvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_base),  INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)    :: mesh
    CLASS(physics_base), INTENT(INOUT) :: physics
    CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
    !------------------------------------------------------------------------!
    ! reconstruct data on cell faces
    IF (this%Reconstruction%PrimRecon()) THEN
       CALL this%Reconstruction%CalculateStates(mesh,physics,pvar,this%prim)
       CALL physics%Convert2Conservative(this%prim,this%cons)
    ELSE
       CALL this%Reconstruction%CalculateStates(mesh,physics,cvar,this%cons)
       CALL physics%Convert2Primitive(this%cons,this%prim)
    END IF

    ! update the speed of sound on cell faces (non-isotherml physics only)
    SELECT TYPE(phys => physics)
    CLASS IS(physics_euler)
      SELECT TYPE(prim => this%prim)
      CLASS IS(statevector_euler)
        CALL phys%UpdateSoundSpeed(prim)
      END SELECT
    END SELECT

    ! get minimal & maximal wave speeds on cell interfaces
    CALL physics%CalculateWaveSpeeds(mesh,this%prim%data5d,this%cons%data5d,this%minwav,this%maxwav)
  END SUBROUTINE calculatefacedata

  SUBROUTINE finalize_base(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fluxes_base), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    IF (.NOT.this%Initialized()) &
        CALL this%Error("CloseFluxes","not initialized")
    CALL this%minwav%Destroy()
    CALL this%maxwav%Destroy()
    CALL this%prim%Destroy()
    CALL this%cons%Destroy()
    CALL this%pfluxes%Destroy()
    DEALLOCATE(this%cons,this%prim,this%pfluxes,this%minwav,this%maxwav, &
         this%bxflux,this%byflux,this%bzflux,this%bxfold,this%byfold,this%bzfold)
    CALL this%Reconstruction%Finalize()
    DEALLOCATE(this%Reconstruction)
  END SUBROUTINE finalize_base

END MODULE fluxes_base_mod