gravity_pointmass.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: gravity_pointmass.f90                                             #
!#                                                                           #
!# Copyright (C) 2007-2019                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Björn Sperling   <sperling@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 gravity_pointmass_mod
  USE gravity_base_mod
  USE boundary_base_mod
  USE fluxes_base_mod
  USE physics_base_mod
  USE mesh_base_mod
  USE logging_base_mod
  USE marray_base_mod
  USE marray_compound_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
  CHARACTER(LEN=16), PARAMETER :: potential_name(2) = (/ &
                                  "Newton          ", &
                                  "Paczynski-Wiita " /)
  INTEGER, PARAMETER :: newton = 1
  INTEGER, PARAMETER :: wiita  = 2

  TYPE, EXTENDS(gravity_base) :: gravity_pointmass
    CLASS(logging_base), ALLOCATABLE    :: potential
    INTEGER                             :: outbound
    REAL, POINTER                       :: mass
    REAL, POINTER                       :: accrate
    REAL, POINTER                       :: massloss
    REAL                                :: mdot
    REAL                                :: acclimit
    REAL                                :: switchon
    REAL, DIMENSION(:,:),       POINTER :: pos
    REAL, DIMENSION(:),         POINTER :: r0
    REAL, DIMENSION(:,:,:),     POINTER :: r_prim
    REAL, DIMENSION(:,:,:,:),   POINTER :: fr_prim
    REAL, DIMENSION(:,:,:,:),   POINTER :: posvec_prim
    REAL, DIMENSION(:,:,:,:),   POINTER :: posvec_prim_tmp
    REAL, DIMENSION(:,:,:,:,:), POINTER :: fposvec_prim
    REAL, DIMENSION(:,:,:,:), POINTER   :: pot_prim
    REAL, DIMENSION(:,:,:),   POINTER   :: omega
    REAL, DIMENSION(:,:,:,:), POINTER   :: omega2
  CONTAINS
    PROCEDURE :: initgravity_pointmass
    PROCEDURE :: calcpotential
    PROCEDURE :: updategravity_single
    PROCEDURE :: infogravity
    PROCEDURE :: calcdiskheight_single
    PROCEDURE :: finalize
  END TYPE

  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       gravity_pointmass, &
       ! constants
       newton, wiita, &
       ! elemental functions
       getdiskheight
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initgravity_pointmass(this,Mesh,Physics,config,IO)
    USE geometry_logcylindrical_mod, ONLY: geometry_logcylindrical
    USE geometry_spherical_mod, ONLY: geometry_spherical
    USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_pointmass), INTENT(INOUT) :: this
    CLASS(mesh_base),         INTENT(IN)    :: Mesh
    CLASS(physics_base),      INTENT(IN)    :: Physics
    TYPE(Dict_TYP),           POINTER       :: config,IO
    !------------------------------------------------------------------------!
    INTEGER :: potential_number, valwrite, gtype
    INTEGER :: err,i,j,k,l
    REAL    :: r_schwarzschild=0.0,eps
    !------------------------------------------------------------------------!
    CALL this%InitGravity(mesh,physics,"central point mass",config,io)
    !\todo last dimensions are absolutely not clear if right.. What are the standing for?
    ALLOCATE(this%potential, &
             this%pot(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,4), &
             this%omega(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
             this%omega2(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,1), &
             this%r_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
             this%fr_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
             this%posvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
             this%posvec_prim_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
             this%fposvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),& ! last two entries (EAST,NORTH,TOP)x(dim1,dim2,dim3)
             this%mass, this%accrate, this%massloss, this%pos(1,3), &
         stat = err)
    IF (err.NE.0) CALL this%Error("InitGravity_pointmass", "Unable allocate memory!")

    ! type of the potential
    CALL getattr(config, "potential", potential_number, newton)
    CALL this%potential%InitLogging(potential_number,potential_name(potential_number))

    ! set (initial) mass
    CALL getattr(config, "mass", this%mass, 1.0)
    CALL setattr(io, "mass", ref(this%mass))

    ! cartesian position of the point mass
    CALL getattr(config, "x", this%pos(1,1), 0.0)
    CALL getattr(config, "y", this%pos(1,2), 0.0)
    CALL getattr(config, "z", this%pos(1,3), 0.0)

    ! Softening (e.g. for planets inside the computational domain)
    CALL getattr(config, "softening", eps, 0.0)

    ! soft switch on
    CALL getattr(config, "switchon", this%switchon, -1.0)

    ! accretion limit (e.g. eddington limit) divided by central mass
    ! the units are therefore (e.g. SI): kg/s / central_mass = 1 / s
    ! -1.0 == off
    CALL getattr(config, "acclimit", this%acclimit, -1.0)

    CALL setattr(io, "accrate", ref(this%accrate))
    IF (this%acclimit.GT.0.0) &
       CALL setattr(io, "massloss", ref(this%massloss))

    SELECT CASE(potential_number)
    CASE(newton)
       r_schwarzschild = 0.0
    CASE(wiita)
       ! Schwarzschild radius
       r_schwarzschild = 2.*physics%constants%GN * this%mass / physics%constants%C**2
    CASE DEFAULT
       CALL this%Error("InitGravity_pointmass", "potential must be either NEWTON or WIITA")
    END SELECT

    ! output cartesian positions of the point mass
    CALL getattr(config, "output/position", valwrite, 0)
    IF (valwrite .EQ. 1) &
       CALL setattr(io, "position", this%pos)

    ! reset mass flux and massloss and register for output
    this%mdot = 0.0
    this%massloss = 0.0
    this%accrate = 0.0

    ! define position vector from the central object to all cell bary centers
    ! and the corresponding distances
    ! TODO Here should be a different evaluation for either 2D or 3D simulation
    IF (abs(this%pos(1,1)).LE.tiny(this%pos(1,1)).AND.abs(this%pos(1,2)).LE.tiny(this%pos(1,2))) THEN
       ! no shift: point mass is located at the origin of the mesh
       ! set position vector and inverse radius using appropriate mesh arrays
       this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
       this%r_prim(:,:,:) = mesh%radius%bcenter(:,:,:)
       this%fr_prim(:,:,:,1) = mesh%radius%faces(:,:,:,east)
       this%fr_prim(:,:,:,2) = mesh%radius%faces(:,:,:,north)
       this%fr_prim(:,:,:,3) = mesh%radius%faces(:,:,:,top)
    ELSE
       ! shifted point mass position:
       ! compute curvilinear components of shift vector
       ! PLEASE NOTE: We need the shifted curvilinear components at the bary_centers and
       ! faces in EASTERN and NORTHERN direction. This is the case because we want to
       ! calculate the potential at these face positions further below.
       this%posvec_prim_tmp(:,:,:,1) = this%pos(1,1)
       this%posvec_prim_tmp(:,:,:,2) = this%pos(1,2)
       this%posvec_prim_tmp(:,:,:,3) = this%pos(1,3)

       CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,east,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,1,:))
       CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,north,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,2,:))
       CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,top,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,3,:))

       CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_prim_tmp,this%posvec_prim)

       ! subtract the result from the position vector:
       ! this gives you the curvilinear components of all vectors pointing
       ! from the point mass to the bary center of any cell on the mesh
       this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_prim(:,:,:,:)
       this%fposvec_prim(:,:,:,1,:) = mesh%posvec%faces(:,:,:,east,:) - this%fposvec_prim(:,:,:,1,:)  ! EAST
       this%fposvec_prim(:,:,:,2,:) = mesh%posvec%faces(:,:,:,north,:) - this%fposvec_prim(:,:,:,2,:) ! NORTH
       this%fposvec_prim(:,:,:,3,:) = mesh%posvec%faces(:,:,:,top,:) - this%fposvec_prim(:,:,:,3,:)   ! TOP

       this%r_prim(:,:,:) = sqrt(this%posvec_prim(:,:,:,1)**2+this%posvec_prim(:,:,:,2)**2+this%posvec_prim(:,:,:,3)**2)
       this%fr_prim(:,:,:,1) = sqrt(this%fposvec_prim(:,:,:,1,1)**2+this%fposvec_prim(:,:,:,1,2)**2+this%fposvec_prim(:,:,:,1,3)**2) ! shifted EAST-faces
       this%fr_prim(:,:,:,2) = sqrt(this%fposvec_prim(:,:,:,2,1)**2+this%fposvec_prim(:,:,:,2,2)**2+this%fposvec_prim(:,:,:,2,3)**2) ! shifted NORTH-faces
       this%fr_prim(:,:,:,3) = sqrt(this%fposvec_prim(:,:,:,3,1)**2+this%fposvec_prim(:,:,:,3,2)**2+this%fposvec_prim(:,:,:,3,3)**2) ! shifted TOP-faces
    END IF

    ! initialize gravitational acceleration and Keplerian angular velocity
    ! loop over ghost cells to ensure that all values are well defined
    ! and remain finite; this is of particular importance for the calculation
    ! of the disk scale height (see below)
!NEC$ COLLAPSE
    DO k=mesh%KGMIN,mesh%KGMAX
      DO j=mesh%JGMIN,mesh%JGMAX
!NEC$ IVDEP
        DO i=mesh%IGMIN,mesh%IGMAX
          ! compute local Keplerian angular velocity
          ! Since omega usually becomes infinite in the limit r -> r_prim (or a, if a > 0 )
          ! (unless softening is enabled, i. e. eps > 0) we avoid devision by zero by first
          ! limiting the inverse of omega to some very small value.
          this%omega2(i,j,k,1) = 1./ max(10*tiny(eps), &
             (( this%r_prim(i,j,k) * (this%r_prim(i,j,k)-r_schwarzschild)**2) + eps**3) / (physics%Constants%GN*this%mass) )
          this%omega(i,j,k) = sqrt(this%omega2(i,j,k,1))
          ! curvilinear components of the gravitational acceleration
!NEC$ UNROLL(3)
          DO l=1,physics%VDIM
            this%accel%data4d(i,j,k,l) = -this%omega2(i,j,k,1) * this%posvec_prim(i,j,k,l)
          END DO
        END DO
      END DO
    END DO

    !! \todo implement computation for pseudo-Newton Paczynski-Wiita potential
    CALL this%CalcPotential(mesh,physics,this%mass,this%r_prim,this%fr_prim,this%pot)


    CALL this%SetOutput(mesh,physics,config,io)

    ! enable mass accretion by setting "outbound" to one of the four boundaries
    ! of the computational domain (depends on mesh geometry)
    SELECT TYPE(geo=>mesh%geometry)
    CLASS IS(geometry_cylindrical)
       CALL getattr(config, "outbound", this%outbound, west)
    CLASS IS(geometry_spherical)
       CALL getattr(config, "outbound", this%outbound, west)
    CLASS DEFAULT
       CALL getattr(config, "outbound", this%outbound, 0)
       CALL this%WARNING("GravitySources_pointmass","geometry does not support accretion")
    END SELECT
    CALL this%InfoGravity()
  END SUBROUTINE initgravity_pointmass

  SUBROUTINE calcpotential(this,Mesh,Physics,mass,dist,dist_faces,pot)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_pointmass)         :: this
    CLASS(mesh_base),    INTENT(IN)  :: Mesh
    CLASS(physics_base), INTENT(IN)  :: Physics
    REAL                             :: mass
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
                         INTENT(IN)  :: dist
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3), &
                         INTENT(IN)  :: dist_faces
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,4), &
                         INTENT(OUT) :: pot
    !------------------------------------------------------------------------!

    !! \todo implement computation for pseudo-Newton Paczynski-Wiita potential
    pot(:,:,:,1) = - physics%constants%GN * mass / dist(:,:,:)
    pot(:,:,:,2) = - physics%constants%GN * mass / dist_faces(:,:,:,1) ! direction EAST !Mesh%radius%faces(:,:,:,EAST)
    pot(:,:,:,3) = - physics%constants%GN * mass / dist_faces(:,:,:,2) ! direction NORTH !Mesh%radius%faces(:,:,:,NORTH)
    pot(:,:,:,4) = - physics%constants%GN * mass / dist_faces(:,:,:,3) ! direction TOP !Mesh%radius%faces(:,:,:,TOP)
  END SUBROUTINE

  SUBROUTINE infogravity(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_pointmass), INTENT(IN) :: this
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32) :: param_str
    !------------------------------------------------------------------------!
    WRITE (param_str,'(ES8.2)') this%mass
    CALL this%Info("            potential:         " // &
         trim(this%potential%GetName()))
    CALL this%Info("            initial mass:      " // trim(param_str))
    IF (this%outbound.GT.0) THEN
      param_str = "enabled"
    ELSE
      param_str = "disabled"
    END IF
    CALL this%Info("            accretion:         " // trim(param_str))
    IF (this%outbound.GT.0.AND.this%acclimit.GT.0.0) THEN
      WRITE (param_str,'(ES8.2)') this%acclimit
      CALL this%Info("            accretion limit:   " // trim(param_str))
    END IF
    IF (this%switchon.GT.0.0) THEN
       WRITE (param_str,'(ES8.2)') this%switchon
       CALL this%Info("            switchon time:     " // trim(param_str))
    END IF
  END SUBROUTINE infogravity


  SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_pointmass), INTENT(INOUT) :: this
    CLASS(mesh_base),         INTENT(IN)    :: Mesh
    CLASS(physics_base),      INTENT(IN)    :: Physics
    CLASS(fluxes_base),       INTENT(IN)    :: Fluxes
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    REAL,                     INTENT(IN)    :: time,dt
    !------------------------------------------------------------------------!
    INTEGER                       :: l
    REAL, DIMENSION(Physics%VNUM) :: bflux
    REAL                          :: scaling,oldmass,massfac,sqrmassfac,dmass,dmasslim
    !------------------------------------------------------------------------!
    ! update accel and omega only in case of accretion
    IF (this%outbound.NE.0) THEN
       ! get boundary flux
#ifdef PARALLEL
       bflux(:)  = fluxes%GetBoundaryFlux(mesh,physics,this%outbound,mpi_comm_world)
#else
       bflux(:)  = fluxes%GetBoundaryFlux(mesh,physics,this%outbound)
#endif
       ! store old mass
       oldmass = this%mass
       ! compute the mass flux difference
       ! REMARK #1: bflux is the accumulated flux across the boundary
       ! REMARK #2: mdot is NOT the accretion rate, but the
       !            mass flux from the previous evalution (see below)
       dmass = (bflux(physics%DENSITY) - this%mdot)
       IF(this%acclimit.GE.0.) THEN
         dmasslim  = min(dmass,this%acclimit*dt*this%mass)
         this%massloss = this%massloss + dmass - dmasslim
         dmass = dmasslim
       END IF
       ! compute new mass and current true, i. e. possibly limited, accretion rate
       this%mass = this%mass + dmass
       IF (dt.GT.0.0) this%accrate = dmass/dt
       ! scale gravitational acceleration and Keplerian angular velocity
       ! with newmass/oldmass to account for accretion
       massfac = this%mass/oldmass
       sqrmassfac = sqrt(massfac)
       this%accel%data4d(:,:,:,:) = this%accel%data4d(:,:,:,:) * massfac
       this%pot(:,:,:,:) = this%pot(:,:,:,:) * massfac
       this%omega(:,:,:) = this%omega(:,:,:) * sqrmassfac
       this%mdot = bflux(physics%DENSITY)
    END IF

    ! modify acceleration during switchon phase
    IF (time.GE.0.0.AND.time.LE.this%switchon) THEN
      ! compute new scaling
      scaling = sin(0.5*pi*time/this%switchon)**2
      ! compute acceleration and scale it;
      ! during the switchon phase accretion is disabled
!NEC$ UNROLL(3)
      DO l=1,physics%VDIM
        this%accel%data4d(:,:,:,l) = -scaling * this%omega2(:,:,:,1) * this%posvec_prim(:,:,:,l)
      END DO
    END IF

  END SUBROUTINE updategravity_single


  PURE SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_pointmass), INTENT(INOUT) :: this
    CLASS(mesh_base),         INTENT(IN)    :: mesh
    CLASS(physics_base),      INTENT(IN)    :: physics
    CLASS(marray_compound),   INTENT(INOUT) :: pvar
    TYPE(marray_base),        INTENT(INOUT) :: bccsound,h_ext,height
    !------------------------------------------------------------------------!
    ! compute disk height
    h_ext%data3d(:,:,:) = getdiskheight(bccsound%data3d(:,:,:),this%omega(:,:,:))
  END SUBROUTINE calcdiskheight_single

  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(gravity_pointmass), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    CHARACTER(LEN=128) :: buffer
    !------------------------------------------------------------------------!
    IF (this%GetRank().EQ.0) THEN
       CALL this%Info("-------------------------------------------------------------------")
       WRITE(buffer, "(A,(ES11.3))") " central mass: ", this%mass
       CALL this%Info(buffer)
       IF (this%acclimit.GT.0.0) THEN
          WRITE(buffer, "(A,(ES11.3))") "    mass loss: ", this%massloss
          CALL this%Info(buffer)
       END IF
    END IF

    DEALLOCATE(this%pot,this%omega,this%omega2,this%r_prim,this%fr_prim, &
               this%posvec_prim,this%posvec_prim_tmp,this%fposvec_prim,this%mass, &
               this%accrate,this%massloss,this%pos)

    DEALLOCATE(this%potential)

    CALL this%Finalize_base()
  END SUBROUTINE finalize

  ELEMENTAL FUNCTION getdiskheight(cs,omega) RESULT(height)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: cs,omega
    REAL             :: height
    !------------------------------------------------------------------------!
    height = cs / omega
  END FUNCTION getdiskheight

END MODULE gravity_pointmass_mod