sources_gravity.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_gravity.f90                                               #
!#                                                                           #
!# Copyright (C) 2014-2018                                                   #
!# Björn Sperling   <sperling@astrophysik.uni-kiel.de>                       #
!# Tobias Illenseer <tillense@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 sources_gravity_mod
  USE logging_base_mod
  USE mesh_base_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE physics_base_mod
  USE sources_base_mod
  USE sources_c_accel_mod
  USE gravity_base_mod
  USE gravity_generic_mod
  USE fluxes_base_mod
  USE boundary_base_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER(LEN=32), PARAMETER :: source_name = "gravity"
  TYPE, EXTENDS(sources_c_accel) :: sources_gravity
    CLASS(gravity_base),    POINTER :: glist => null() 
     LOGICAL                        :: addtoenergy
  CONTAINS
    PROCEDURE :: initsources_gravity
    PROCEDURE :: infosources
    PROCEDURE :: updategravity
    PROCEDURE :: externalsources_single
    PROCEDURE :: calcdiskheight
    PROCEDURE :: finalize
  END TYPE sources_gravity
  abstract INTERFACE
  END INTERFACE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       sources_gravity
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initsources_gravity(this,Mesh,Physics,Fluxes,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_gravity), INTENT(INOUT) :: this
    CLASS(mesh_base),       INTENT(IN)    :: Mesh
    CLASS(fluxes_base),     INTENT(IN)    :: Fluxes
    CLASS(physics_base),    INTENT(IN)    :: Physics
    TYPE(Dict_TYP),         POINTER       :: config, IO
    !------------------------------------------------------------------------!
    INTEGER :: err, stype,i, valwrite
    !------------------------------------------------------------------------!
    CALL getattr(config,"stype",stype)
    CALL this%InitLogging(stype,source_name)
    CALL this%InitSources(mesh,fluxes,physics,config,io)

    this%accel = marray_base(physics%VDIM)
    this%accel%data1d(:) = 0.
    this%pot = marray_base(4)
    this%pot%data1d(:) = 0.

    ! Add source terms to energy equation?
    ! Set this to zero, if a potential is defined in physics_euler2Diamt
    CALL getattr(config, "energy", i, 1)
    IF(i.EQ.0) THEN
      this%addtoenergy = .false.
    ELSE
      this%addtoenergy = .true.
    END IF

    ! enable update of disk scale height if requested
    CALL getattr(config, "update_disk_height", i, 0)
    IF (i.EQ.1) THEN
      IF (physics%VDIM.EQ.2) THEN
        this%update_disk_height = .true.

        this%height = marray_base()
        this%h_ext= marray_base()
        this%invheight2 = marray_base()
        this%height%data1d = 0.0
        this%h_ext%data1d = 0.0
        this%invheight2%data1d = 0.0
      ELSE
        CALL this%Error("InitGravity", "DiskHeight is only supported in 2D")
      END IF
    ELSE
       this%update_disk_height = .false.
    END IF

    ! initialize gravity
    CALL new_gravity(this%glist,mesh,fluxes,physics,config,io)

    CALL getattr(config, "output/height", valwrite, 0)
    IF (valwrite .EQ. 1) THEN
       CALL setattr(config, "update_disk_height", 1)
       IF (.NOT.ASSOCIATED(this%height%data1d)) THEN
          ALLOCATE(this%height%data3d(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
                   stat=err)
          IF (err.NE.0) CALL this%Error("SetOutput", &
                                   "Memory allocation failed for this%height!")
       END IF
       CALL setattr(io, "height", &
         this%height%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
    END IF

  END SUBROUTINE initsources_gravity

  SUBROUTINE externalsources_single(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    USE physics_euler_mod, ONLY : statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_gravity), 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
    !------------------------------------------------------------------------!
    ! go through all gravity terms in the list
    CALL this%UpdateGravity(mesh,physics,fluxes,pvar,time,dt)

    ! update disk scale height if requested
    IF (this%update_disk_height) THEN
      SELECT TYPE(phys => physics)
      CLASS IS (physics_eulerisotherm)
        CALL this%CalcDiskHeight(mesh,phys,pvar)
      END SELECT
    END IF

    ! gravitational source terms
    CALL physics%ExternalSources(this%accel,pvar,cvar,sterm)

    ! Set src term in energy equation to zero, if it is handeled in the physics
    ! module
    SELECT TYPE(s => sterm)
    TYPE IS (statevector_euler)
      IF (.NOT.this%addtoenergy) s%energy%data1d(:) = 0.
    END SELECT

  END SUBROUTINE externalsources_single

  SUBROUTINE updategravity(this,Mesh,Physics,Fluxes,pvar,time,dt)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_gravity), TARGET, 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
    REAL,                INTENT(IN)    :: time, dt
    !------------------------------------------------------------------------!
    CLASS(gravity_base), POINTER       :: gravptr
    !------------------------------------------------------------------------!
    ! update acceleration of all gravity sources
    ! reset gterm
    this%pot%data1d(:) = 0.
    this%accel%data1d(:) = 0.

    gravptr => this%glist
    DO WHILE (ASSOCIATED(gravptr))
      ! get gravitational acceleration of each gravity module
      CALL gravptr%UpdateGravity_single(mesh,physics,fluxes,pvar,time,dt)

      ! add contribution to the overall gravitational acceleration
      this%accel = this%accel + gravptr%accel

      ! add potential if available
      IF(ASSOCIATED(gravptr%pot)) &
        this%pot%data4d(:,:,:,:) = this%pot%data4d(:,:,:,:) + gravptr%pot(:,:,:,:)

      ! next source term
      gravptr => gravptr%next
    END DO
  END SUBROUTINE

  SUBROUTINE calcdiskheight(this,Mesh,Physics,pvar)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    USE gravity_pointmass_mod, ONLY : gravity_pointmass
    USE gravity_spectral_mod, ONLY : gravity_spectral
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_gravity),TARGET,INTENT(INOUT) :: this
    CLASS(mesh_base),             INTENT(IN)    :: Mesh
    CLASS(physics_eulerisotherm), INTENT(INOUT) :: Physics
    CLASS(marray_compound),       INTENT(INOUT) :: pvar
    !------------------------------------------------------------------------!
    CLASS(gravity_base), POINTER :: grav_ptr,selfgrav_ptr => null()
    LOGICAL                      :: has_external_potential = .false.
    !------------------------------------------------------------------------!
    ! reset inverse scale height^2
    ! go through all gravity terms in the list
    grav_ptr => this%glist
    DO WHILE(ASSOCIATED(grav_ptr))
      SELECT TYPE (grav => grav_ptr)
      CLASS IS (gravity_pointmass)
!CDIR IEXPAND
        CALL grav%CalcDiskHeight_single(mesh,physics,pvar,physics%bccsound,this%h_ext,this%height)
        this%invheight2%data1d(:) = this%invheight2%data1d(:) + 1./this%h_ext%data1d(:)**2
        has_external_potential = .true.
      CLASS IS (gravity_spectral)
        selfgrav_ptr => grav_ptr
      END SELECT

      ! next gravity term
      grav_ptr => grav_ptr%next
    END DO

    ! self-gravity of the disk needs special treatment
    IF (ASSOCIATED(selfgrav_ptr)) THEN
      IF (has_external_potential) THEN
        ! compute the resultant height due to all external gravitational forces
        this%h_ext%data1d(:) = 1./sqrt(this%invheight2%data1d(:))
      END IF
      CALL selfgrav_ptr%CalcDiskHeight_single(mesh,physics,pvar,physics%bccsound,this%h_ext,this%height)
    ELSE
      ! non-selfgravitating disk
      this%height%data1d(:) = 1./sqrt(this%invheight2%data1d(:))
    END IF

  END SUBROUTINE calcdiskheight



  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_gravity), INTENT(IN) :: this
    CLASS(mesh_base),       INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
  END SUBROUTINE infosources

  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_gravity), INTENT(INOUT) :: this
    CLASS(gravity_base),    POINTER       :: gravptr
    !------------------------------------------------------------------------!
    CALL this%pot%Destroy()
    CALL this%accel%Destroy()
    CALL this%height%Destroy()
    CALL this%invheight2%Destroy()
    CALL this%h_ext%Destroy()

    gravptr => this%glist
    DO WHILE (ASSOCIATED(gravptr))

      CALL gravptr%Finalize()

      gravptr => gravptr%next
    END DO

    CALL this%sources_c_accel%Finalize()
  END SUBROUTINE finalize

END MODULE sources_gravity_mod