sources_gravity.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_gravity.f90                                               #
!#                                                                           #
!# Copyright (C) 2014-2024                                                   #
!# 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() 
    TYPE(marray_base)               :: pot
    TYPE(marray_base), ALLOCATABLE  :: height
    TYPE(marray_base), ALLOCATABLE  :: invheight2
    TYPE(marray_base), ALLOCATABLE  :: h_ext
     LOGICAL                        :: addtoenergy
     LOGICAL                        :: update_disk_height
  CONTAINS
    PROCEDURE :: initsources
    PROCEDURE :: updategravity
    PROCEDURE :: externalsources
    PROCEDURE :: calcdiskheight
    final :: finalize
  END TYPE sources_gravity
  abstract INTERFACE
  END INTERFACE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       sources_gravity
  !--------------------------------------------------------------------------!
 
CONTAINS
 
  SUBROUTINE initsources(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%InitSources_base(stype,source_name)
 
    ALLOCATE(this%accel)
    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
 
    ! initialize gravity
    CALL new_gravity(this%glist,mesh,fluxes,physics,config,io)
 
    ! activate output of scale height if requested (default: disabled)
    CALL getattr(config, "output/height", valwrite, 0)
 
    ! this parameter is usually set automatically depending on other
    ! settings (see sources_generic), but if output of height
    ! is enabled, computation of disk height is activated as well
    ! (works only 1D/2D accretion disk simulations with flat geometry)
    IF (valwrite.EQ.1) CALL setattr(config, "update_disk_height", 1)
    CALL getattr(config, "update_disk_height",valwrite, 0)
 
    IF (valwrite.EQ.1) THEN
      IF (physics%VDIM.EQ.2) THEN
        this%update_disk_height = .true.
        ALLOCATE(this%height,this%h_ext,this%invheight2)
        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
 
        IF (valwrite.EQ.1) THEN
           CALL setattr(io, "height", &
             this%height%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
        END IF
      ELSE
        CALL this%Error("sources_gravity::InitSources", "Computation of disk scale height is only supported in 1D/2D flat geometries")
      END IF
    ELSE
       this%update_disk_height = .false.
    END IF
 
    ! activate output of gravitational acceleration
    CALL getattr(config, "output/accel", valwrite, 0)
    IF (valwrite .EQ. 1) THEN
       CALL setattr(io, "accel", &
          this%accel%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1:physics%VDIM))
    END IF
 
  END SUBROUTINE initsources
 
  SUBROUTINE externalsources(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
 
  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
    !------------------------------------------------------------------------!
    INTEGER                            :: i,j,k,l
    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 (ALLOCATED(gravptr%pot)) THEN
        this%pot = this%pot + gravptr%pot
      END IF
 
      ! 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
    this%invheight2%data1d(:) = 0.0
    ! 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)
        CALL grav%CalcDiskHeight_single(mesh,physics,pvar,physics%bccsound,this%h_ext,this%height)
        this%invheight2%data1d(:) = this%invheight2%data1d(:) + 1./this%height%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 finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(sources_gravity),  INTENT(INOUT) :: this
    CLASS(gravity_base),    POINTER       :: gravptr,tmpptr
    !------------------------------------------------------------------------!
    gravptr => this%glist
    DO WHILE (ASSOCIATED(gravptr))
      tmpptr => gravptr
      gravptr => gravptr%next
      DEALLOCATE(tmpptr)
    END DO
    IF (ALLOCATED(this%height)) DEALLOCATE(this%height)
    IF (ALLOCATED(this%h_ext)) DEALLOCATE(this%h_ext)
    IF (ALLOCATED(this%invheight2)) DEALLOCATE(this%invheight2)
    ! deallocation of this%accel is done in inherited destructor
    ! which is called automatically
  END SUBROUTINE finalize
 
END MODULE sources_gravity_mod