sources_generic.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: mesh_generic.f90                                                  #
!#                                                                           #
!# Copyright (C) 2016-2024                                                   #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Tobias Illenseer <tillense@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 sources_generic_mod
  USE sources_base_mod
  USE sources_c_accel_mod
  USE sources_cooling_mod
  USE sources_diskcooling_mod
  USE sources_gravity_mod
  USE sources_planetheating_mod
  USE sources_planetcooling_mod
  USE sources_rotframe_mod
  USE sources_shearbox_mod
  USE sources_viscosity_mod
  USE marray_compound_mod
  USE marray_base_mod
  USE mesh_base_mod
  USE fluxes_base_mod
  USE physics_base_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  TYPE, EXTENDS(sources_base) :: sources_list
    CLASS(marray_compound), POINTER :: sterm => null() 
 
  CONTAINS
    PROCEDURE :: initsources
    PROCEDURE :: externalsources
    PROCEDURE :: calctimestep
    final     :: finalize
  END TYPE sources_list
  !--------------------------------------------------------------------------!
  ! flags for source terms
  INTEGER, PARAMETER :: list             = 0
  INTEGER, PARAMETER :: gravity          = 1
!  INTEGER, PARAMETER :: DISK_THOMSON     = 2
  INTEGER, PARAMETER :: viscosity        = 3
  INTEGER, PARAMETER :: c_accel          = 4
  INTEGER, PARAMETER :: cooling          = 5
  INTEGER, PARAMETER :: rotating_frame   = 20
!  INTEGER, PARAMETER :: SGS              = 23
  INTEGER, PARAMETER :: disk_cooling     = 24
!  INTEGER, PARAMETER :: WAVE_DAMPING     = 25
!  INTEGER, PARAMETER :: FORCING          = 26
  INTEGER, PARAMETER :: planet_heating   = 27
  INTEGER, PARAMETER :: planet_cooling   = 28
!  INTEGER, PARAMETER :: STELLAR_HEATING  = 29
  INTEGER, PARAMETER :: shearbox         = 30
  !--------------------------------------------------------------------------!
  PUBLIC :: &
    ! types
    sources_list, &
 
    ! constants
    viscosity, c_accel, shearbox, gravity, disk_cooling, rotating_frame, &
    planet_heating, planet_cooling, cooling, &
    molecular,alpha,beta,powerlaw,alpha_alt,viscosity_name, &
    gray, gammie, gammie_sb
 
 
CONTAINS
 
  SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_list), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN)    :: Mesh
    CLASS(physics_base), INTENT(IN)    :: Physics
    CLASS(fluxes_base),  INTENT(IN)    :: Fluxes
    TYPE(dict_typ), POINTER            :: config, IO
    !------------------------------------------------------------------------!
    CLASS(sources_base), POINTER :: newsrc => null(), tmpsrc => null()
    TYPE(dict_typ),      POINTER :: dir,next,src,IOsrc,gsrc => null(),gdir => null()
    INTEGER                      :: update_disk_height = 0
    INTEGER                      :: stype
    !------------------------------------------------------------------------!
    IF (.NOT.physics%Initialized().OR..NOT.mesh%Initialized()) &
         CALL this%Error("sources_generic::InitSources","physics and/or mesh module uninitialized")
    IF (this%Initialized()) &
      CALL this%Error("sources_generic::InitSources","list of source term already initialized")
 
    CALL this%InitLogging(list,"sources list")
 
    ALLOCATE(this%sterm,stat=this%err)
    IF (this%err.NE.0) &
      CALL this%Error("sources_generic::InitSources","memory allocation failed")
 
    CALL physics%new_statevector(this%sterm,conservative)
 
    dir => config
    DO WHILE(ASSOCIATED(dir))
      NULLIFY(iosrc)
      IF(haschild(dir)) THEN
        src => getchild(dir)
        CALL getattr(src, "stype", stype)
 
        ! object creation
        SELECT CASE(stype)
        CASE(gravity)
           ! skip initialization of gravity modules here and initialize them
           ! at the end to make sure gravity is the first source term in the list
           gsrc => src
           gdir => dir
           NULLIFY(newsrc)
        CASE(c_accel)
          ALLOCATE(sources_c_accel::newsrc)
        CASE(cooling)
          ALLOCATE(sources_cooling::newsrc)
        CASE(disk_cooling)
          ALLOCATE(sources_diskcooling::newsrc)
        CASE(planet_heating)
          ALLOCATE(sources_planetheating::newsrc)
        CASE(planet_cooling)
          ALLOCATE(sources_planetcooling::newsrc)
        CASE(rotating_frame)
          tmpsrc => this%GetSourcesPointer(rotating_frame)
          IF (ASSOCIATED(tmpsrc)) &
            CALL this%Error("sources_generic::InitSources","only one rotating frame source term allowed")
          ALLOCATE(sources_rotframe::newsrc)
        CASE(shearbox)
          tmpsrc => this%GetSourcesPointer(shearbox)
          IF (ASSOCIATED(tmpsrc)) &
            CALL this%Error("sources_generic::InitSources","only one shearing box source term allowed")
          ALLOCATE(sources_shearbox::newsrc)
        CASE(viscosity)
          ALLOCATE(sources_viscosity::newsrc)
        CASE DEFAULT
          CALL this%Error("sources_generic::InitSources","Unknown source type")
        END SELECT
 
        ! basic initialization of all source terms except gravity
        IF (ASSOCIATED(newsrc)) THEN
          CALL newsrc%InitSources(mesh,physics,fluxes,src,iosrc)
          SELECT TYPE(obj => newsrc)
          TYPE IS (sources_diskcooling)
            IF (obj%cooling%GetType().EQ.gray) update_disk_height = 1
          TYPE IS (sources_viscosity)
            IF (obj%viscosity%GetType().EQ.alpha_alt) update_disk_height = 1
          END SELECT
 
          ! prepend new source term to list
          newsrc%next => this%next
          this%next => newsrc
        END IF
 
        ! process the output dictionary
        IF(ASSOCIATED(iosrc)) CALL setattr(io, getkey(dir), iosrc)
 
      END IF
      dir => getnext(dir)
    END DO
 
    ! finally initialize gravity
    IF(ASSOCIATED(gsrc)) THEN
      NULLIFY(iosrc)
 
      ALLOCATE(sources_gravity::newsrc)
      IF (ASSOCIATED(newsrc)) THEN
        CALL setattr(gsrc,"update_disk_height", update_disk_height)
        CALL newsrc%InitSources(mesh,physics,fluxes,gsrc,iosrc)
        ! prepend new source term to list
        newsrc%next => this%next
        this%next => newsrc
        ! check for disk cooling source term and set pointer to gravity source
        tmpsrc => this%GetSourcesPointer(disk_cooling)
        IF (ASSOCIATED(tmpsrc)) THEN
          SELECT TYPE(sp => tmpsrc)
          CLASS IS(sources_diskcooling)
            SELECT TYPE(grav => newsrc)
            CLASS IS(sources_gravity)
              sp%grav => grav
            END SELECT
          END SELECT
        END IF
      END IF
 
      IF(ASSOCIATED(iosrc)) CALL setattr(io, getkey(gdir), iosrc)
    END IF
 
  END SUBROUTINE initsources
 
 
  SUBROUTINE externalsources(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_list), INTENT(INOUT)      :: this
    CLASS(mesh_base),    INTENT(IN)         :: Mesh
    CLASS(physics_base), INTENT(INOUT)      :: Physics
    CLASS(fluxes_base),  INTENT(IN)         :: Fluxes
    CLASS(sources_base), INTENT(INOUT)      :: Sources
    REAL,                INTENT(IN)         :: time,dt
    CLASS(marray_compound), INTENT(INOUT)   :: pvar,cvar,sterm
    !------------------------------------------------------------------------!
    CLASS(sources_base), POINTER :: srcptr
    !------------------------------------------------------------------------!
    ! reset sterm
    sterm%data1d(:) = 0.0
    ! go through all source terms in the list
    srcptr => this%next
    DO WHILE (ASSOCIATED(srcptr))
 
      CALL srcptr%ExternalSources(mesh,physics,fluxes,this,time,dt,pvar,cvar,this%sterm)
 
      ! add to the sources
      sterm%data1d(:) = sterm%data1d(:) + this%sterm%data1d(:)
 
      ! next source term
      srcptr => srcptr%next
 
    END DO
    ! reset ghost cell data
    IF (mesh%GINUM.GT.0) THEN
      sterm%data4d(mesh%IGMIN:mesh%IMIN-mesh%IP1,:,:,:) = 0.0
      sterm%data4d(mesh%IMAX+mesh%IP1:mesh%IGMAX,:,:,:) = 0.0
    END IF
    IF (mesh%GJNUM.GT.0) THEN
      sterm%data4d(:,mesh%JGMIN:mesh%JMIN-mesh%JP1,:,:) = 0.0
      sterm%data4d(:,mesh%JMAX+mesh%JP1:mesh%JGMAX,:,:) = 0.0
    END IF
    IF (mesh%GKNUM.GT.0) THEN
      sterm%data4d(:,:,mesh%KGMIN:mesh%KMIN-mesh%KP1,:) = 0.0
      sterm%data4d(:,:,mesh%KMAX+mesh%KP1:mesh%KGMAX,:) = 0.0
    END IF
  END SUBROUTINE externalsources
 
 
  SUBROUTINE calctimestep(this,Mesh,Physics,Fluxes,pvar,cvar,time,dt,dtcause)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_list),         INTENT(INOUT) :: this
    CLASS(mesh_base),            INTENT(IN)    :: Mesh
    CLASS(physics_base),         INTENT(INOUT) :: Physics
    CLASS(fluxes_base),          INTENT(IN)    :: Fluxes
    CLASS(marray_compound),      INTENT(INOUT) :: pvar,cvar
    REAL, INTENT(IN)              :: time
    REAL, INTENT(INOUT)           :: dt
    INTEGER, INTENT(OUT)          :: dtcause
    !------------------------------------------------------------------------!
    CLASS(sources_base), POINTER :: srcptr
    REAL              :: dt_new
    !------------------------------------------------------------------------!
    dt_new = dt
 
    ! go through all source terms in the list
    srcptr => this%next
    DO WHILE(ASSOCIATED(srcptr))
 
       CALL srcptr%CalcTimestep(mesh,physics,fluxes,pvar,cvar,time,dt_new,dtcause)
 
 
       IF (dt_new .LT. dt) dtcause=srcptr%GetType()
       dt = min(dt,dt_new)
       ! next source term
       srcptr => srcptr%next
    END DO
  END SUBROUTINE calctimestep
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(sources_list) :: this
    !------------------------------------------------------------------------!
    CLASS(sources_base), POINTER :: srcptr,next
    !------------------------------------------------------------------------!
    ! loop over all source terms and call finalization subroutines
    srcptr => this%next
    DO WHILE (ASSOCIATED(srcptr))
      next => srcptr%next
      DEALLOCATE(srcptr)
      srcptr => next
    END DO
    IF (ASSOCIATED(this%sterm)) DEALLOCATE(this%sterm)
  END SUBROUTINE finalize
 
END MODULE sources_generic_mod