sources_diskcooling.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_diskcooling.f90                                           #
!#                                                                           #
!# Copyright (C) 2011-2021                                                   #
!# 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_diskcooling_mod
  USE sources_base_mod
  USE sources_cooling_mod
  USE sources_gravity_mod
  USE physics_base_mod
  USE fluxes_base_mod
  USE mesh_base_mod
  USE marray_compound_mod
  USE marray_base_mod
  USE logging_base_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER(LEN=32), PARAMETER :: source_name = "thin accretion disk cooling"
  INTEGER, PARAMETER :: gray   = 1
  INTEGER, PARAMETER :: gammie = 2
  INTEGER, PARAMETER :: gammie_sb = 3
  CHARACTER(LEN=28), PARAMETER, DIMENSION(3) :: &
                        cooling_name = (/ "Gray disk cooling          ", &
                                          "Gammie disk cooling        ", &
                                          "Gammie shearing box cooling" /)
  !--------------------------------------------------------------------------!
  REAL, PARAMETER :: sqrt_three = 1.73205080757
  REAL, PARAMETER :: sqrt_twopi = 2.50662827463
 
! Rosseland mean opacity constants in SI units;
! taken from Bell & Lin, ApJ, 427, 1994
! kappa_i= kappa_0i * rho**rexp(i) * T**Texp(i)
 
!!$  DOUBLE PRECISION, PARAMETER :: kappa0(8) = (/ &
!!$       2.00D-05, & ! ice grains                     [m^2/kg/K^2]
!!$       2.00D+15, & ! evaporation of ice grains      [m^2/kg*K^7]
!!$       1.00D-02, & ! metal grains                 [m^2/kg/K^0.5]
!!$       2.00D+77, & ! evaporation of metal grains [m^5/kg^2*K^24]
!!$       1.00D-11, & ! molecules                [m^4/kg^(5/3)/K^3]
!!$       1.00D-38, & ! H-scattering            [m^3/kg^(4/3)/K^10]
!!$       1.50D+16, & ! bound-free and free-free [m^5/kg^2*K^(5/2)]
!!$       KE /)       ! electron scattering                [m^2/kg]
  REAL, PARAMETER :: logkappa0(8) = (/ &
       -10.8197782844, & ! ice grains                     [m^2/kg/K^2]
       35.2319235755, &  ! evaporation of ice grains      [m^2/kg*K^7]
       -4.60517018599, & ! metal grains                 [m^2/kg/K^0.5]
       177.992199341, &  ! evaporation of metal grains [m^5/kg^2*K^24]
       -25.3284360229, & ! molecules                [m^4/kg^(5/3)/K^3]
       -87.4982335338, & ! H-scattering            [m^3/kg^(4/3)/K^10]
       37.246826596, &   ! bound-free and free-free [m^5/kg^2*K^(5/2)]
       -3.3581378922 /)  ! electron scattering                [m^2/kg]
  REAL, PARAMETER :: texp(8) = (/ 2.0, -7.0, 0.5, -24.0, 3.0, 10.0, -2.5, 0.0 /)
  REAL, PARAMETER :: rexp(8) = (/ 0.0, 0.0, 0.0, 1.0, 2./3., 1./3., 1.0, 0.0 /)
  REAL, PARAMETER :: t0 = 3000      ! temperature constant (opacity interpolation)
 
  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(sources_cooling) :: sources_diskcooling
    CLASS(logging_base), ALLOCATABLE :: cooling
    TYPE(marray_base), ALLOCATABLE  :: ephir
    TYPE(sources_gravity), POINTER  :: grav
    REAL                            :: b_cool
    REAL                            :: b_start
    REAL                            :: b_final
    REAL                            :: t_start
    REAL                            :: dt_bdec
    REAL                            :: t_0
    REAL                            :: rho_0
  CONTAINS
    PROCEDURE :: initsources
    PROCEDURE :: infosources
    PROCEDURE :: updatecooling
    final :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       sources_diskcooling, &
       ! constants
       gray, gammie, gammie_sb
  !--------------------------------------------------------------------------!
 
CONTAINS
 
  SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
    USE physics_euler_mod, ONLY : physics_euler
    USE geometry_cylindrical_mod, ONLY : geometry_cylindrical
    USE geometry_cartesian_mod, ONLY : geometry_cartesian
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_diskcooling), 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
    !------------------------------------------------------------------------!
    INTEGER :: stype,err
    INTEGER :: cooling_func
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    CALL this%InitSources_base(stype,source_name)
 
    ! some sanity checks
    ! isothermal modules are excluded
    SELECT TYPE (phys => physics)
    CLASS IS(physics_euler)
      ! do nothing
    CLASS DEFAULT
      ! abort
      CALL this%Error("sources_diskcooling::InitSources","physics not supported")
    END SELECT
 
    ! check for dimensions
    IF (physics%VDIM.GT.2) THEN
      CALL this%Error("sources_diskcooling::InitSources","only 2D simulations supported")
    END IF
 
    ! get cooling method
    CALL getattr(config,"method",cooling_func)
    IF (cooling_func.LE.1.OR.cooling_func.GT.3) &
      CALL this%Error("sources_diskcooling::InitSources","unknown cooling method")
 
    ALLOCATE(logging_base::this%cooling,stat=err)
    IF (err.EQ.0) ALLOCATE(this%Q,stat=err)
    IF (err.NE.0) CALL this%Error("sources_diskcooling::InitSources","memory allocation failed")
    CALL this%cooling%InitLogging(cooling_func,cooling_name(cooling_func))
    this%Q = marray_base()
 
    SELECT CASE(this%cooling%GetType())
    CASE(gray)
      ! check units, currently only SI units are supported
      IF (physics%constants%GetType().NE.si) &
         CALL this%Error("sources_diskcooling::InitSources","only SI units supported for gray cooling")        
    CASE(gammie)
      ! check for geometry
      SELECT TYPE (mgeo => mesh%Geometry)
      CLASS IS(geometry_cylindrical)
        ALLOCATE(this%ephir)
        this%ephir = marray_base(2)
        this%ephir%data4d(:,:,:,1) = -mesh%posvec%bcenter(:,:,:,2)/mesh%radius%bcenter(:,:,:)**2
        this%ephir%data4d(:,:,:,2) = mesh%posvec%bcenter(:,:,:,1)/mesh%radius%bcenter(:,:,:)**2
      CLASS DEFAULT
        CALL this%Error("sources_diskcooling::InitSources","geometry not supported")
      END SELECT
    CASE(gammie_sb)
      SELECT TYPE (mgeo => mesh%Geometry)
      TYPE IS(geometry_cartesian)
        ! do nothing
      CLASS DEFAULT
        CALL this%Error("sources_diskcooling::InitSources","geometry not supported")
      END SELECT
    CASE DEFAULT
      ! this should not happen (see above)
      CALL this%Error("sources_diskcooling::InitSources","unknown cooling method")
    END SELECT
 
    ! Courant number, i.e. safety factor for numerical stability
    CALL getattr(config, "cvis", this%cvis, 0.1)
 
    ! minimum temperature
    CALL getattr(config, "Tmin", this%T_0, 1.0e-30)
 
    ! minimum density
    CALL getattr(config, "Rhomin", this%rho_0, 1.0e-30)
 
    ! cooling time
    CALL getattr(config, "b_cool", this%b_cool, 1.0e+00)
 
    ! initial and final cooling time
    CALL getattr(config, "b_start", this%b_start, 0.0)
    CALL getattr(config, "b_final", this%b_final, this%b_cool)
 
    ! starting point for beta
    CALL getattr(config, "t_start", this%t_start, 0.0)
 
    ! timescale for changing beta
    CALL getattr(config, "dt_bdec", this%dt_bdec, -1.0)
 
    ! set initial time < 0
    this%time = -1.0
 
    ! initialize cooling function
    this%Q%data1d(:)  = 0.0
 
    !initialise output
    CALL this%SetOutput(mesh,config,io)
    CALL this%InfoSources(mesh)
  END SUBROUTINE initsources
 
  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_diskcooling), INTENT(IN) :: this
    CLASS(mesh_base),           INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32) :: param_str
    !------------------------------------------------------------------------!
    CALL this%Info("            cooling function:  " // trim(this%cooling%GetName()))
    SELECT CASE(this%cooling%GetType())
    CASE(gray)
       WRITE (param_str,'(ES10.2)') this%T_0
       CALL this%Info("            min. temperature:  " // trim(param_str))
       WRITE (param_str,'(ES10.2)') this%rho_0
       CALL this%Info("            minimum density:   " // trim(param_str))
    CASE(gammie)
       WRITE (param_str,'(ES10.2)') this%b_cool
       CALL this%Info("            cooling parameter: " // trim(param_str))
    CASE(gammie_sb)
       WRITE (param_str,'(ES10.2)') this%b_cool
       CALL this%Info("            cooling parameter: " // trim(param_str))
       IF(this%dt_bdec.GE.0.0) THEN
        WRITE (param_str,'(ES10.2)') this%b_start
        CALL this%Info("            initial cooling parameter: " // trim(param_str))
        WRITE (param_str,'(ES10.2)') this%b_final
        CALL this%Info("            final cooling parameter: " // trim(param_str))
        WRITE (param_str,'(ES10.2)') this%t_start
        CALL this%Info("            starting b_dec time: " // trim(param_str))
        WRITE (param_str,'(ES10.2)') this%dt_bdec
        CALL this%Info("            operating time: " // trim(param_str))
       END IF
    END SELECT
  END SUBROUTINE infosources
 
 
  SUBROUTINE updatecooling(this,Mesh,Physics,time,pvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_diskcooling),   INTENT(INOUT) :: this
    CLASS(mesh_base),             INTENT(IN)    :: Mesh
    CLASS(physics_euler),         INTENT(IN)    :: Physics
    REAL,                         INTENT(IN)    :: time
    CLASS(statevector_euler),     INTENT(IN)    :: pvar
    !------------------------------------------------------------------------!
    REAL              :: muRgamma,Qfactor
    !------------------------------------------------------------------------!
 
    ! energy loss due to radiation processes
    SELECT CASE(this%cooling%GetType())
       CASE(gray)
          ! make sure the
          IF (ASSOCIATED(this%grav)) THEN
            IF (ALLOCATED(this%grav%height)) THEN
              IF (ASSOCIATED(this%grav%height%data1d)) THEN
                ! some constants
                murgamma = physics%mu/(physics%Constants%RG*physics%gamma)
                qfactor  = 8./3.*physics%Constants%SB
                ! compute gray cooling term
                this%Q%data1d(:) = lambda_gray(pvar%density%data1d(:),this%grav%height%data1d(:), &
                    murgamma*physics%bccsound%data1d(:)*physics%bccsound%data1d(:), &
                    this%rho_0,this%T_0,qfactor)
              ELSE
                CALL this%Error("sources_diskcooling::UpdateCooling","disk height seems not properly assigned in gray cooling")
              END IF
            END IF
          END IF
       CASE(gammie)
          ! compute Gammie cooling term based on cooling time scale
          ! t_cool = b_cool / Omega with Omega = ephi*v / r and
          ! dimensionless scaling parameter beta which is usually
          ! held constant except for the initial phase
          IF (this%dt_bdec.GE.0.0) THEN
            IF (time .LT. this%t_start) THEN
              this%b_cool = this%b_start
            ELSE IF ((time .GE. this%t_start) .AND. ((time-this%t_start) .LT. &
            this%dt_bdec)) THEN
              this%b_cool = this%b_start +  &
                  (this%b_final-this%b_start)/this%dt_bdec*(time - this%t_start)
            ELSE IF ((time-this%t_start) .GE. this%dt_bdec) THEN
              this%b_cool = this%b_final
            END IF
          END IF
          ! set cooling function
          this%Q%data3d(:,:,:) = lambda_gammie(pvar%pressure%data3d(:,:,:) / (physics%gamma-1.), &
              abs((this%ephir%data4d(:,:,:,1)*pvar%velocity%data4d(:,:,:,1) &
                  +this%ephir%data4d(:,:,:,2)*(pvar%velocity%data4d(:,:,:,2)&
                    +mesh%OMEGA*mesh%radius%bcenter(:,:,:)))) / this%b_cool)
       CASE(gammie_sb)
          ! in sb t_cool = b_cool
          this%Q%data1d(:) = lambda_gammie(pvar%pressure%data1d(:) / (physics%gamma-1.), &
              mesh%OMEGA/this%b_cool)
    END SELECT
 
  END SUBROUTINE updatecooling
 
  ELEMENTAL FUNCTION rosselandmeanopacity_new(logrho,logT) RESULT(kappa_R)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    REAL, INTENT(IN) :: logrho,logt
    REAL             :: kappa_r
    !------------------------------------------------------------------------!
    REAL :: kappa4(8)
    REAL :: tfactor
    !------------------------------------------------------------------------!
    ! compute kappa_i^4 for all cooling mechanisms
!NEC$ unroll(8)
    kappa4(:) = exp(max(-40.0,min(40.0, &
         4.*(logkappa0(:)+rexp(:)*logrho+texp(:)*logt))))
    ! compute (T/T0)**10
    tfactor = exp(max(-40.0,min(40.0,10.*(logt-log(t0)))))
    ! compute Rosseland mean using Gails interpolation formula
    kappa_r = 1. /(tiny(kappa_r) + &
           (1./kappa4(1) + 1./(1.+tfactor)/(kappa4(2)+kappa4(3)))**0.25 &
         + (1./(kappa4(4)+kappa4(5)+kappa4(6)) + 1./(kappa4(7)+kappa4(8)))**0.25)
  END FUNCTION rosselandmeanopacity_new
 
  ELEMENTAL FUNCTION lambda_gray(Sigma,h,Tc,rho0,T0,Qf) RESULT(Qcool)
    IMPLICIT NONE
    !---------------------------------------------------------------------!
    REAL, INTENT(IN) :: sigma,h,tc,rho0,t0,qf
    REAL :: qcool
    REAL :: logrho,logt,kappa,tau,tau_eff,t0tc
    !---------------------------------------------------------------------!
    ! logarithm of midplane density
    ! log(SQRT(2*Pi)^-1 * Sigma / H ) = -log(SQRT(2*Pi) * H / Sigma)
    logrho = -log(max(rho0,sqrt_twopi * h / sigma))
    ! logarithm of midplane temperature
    logt = log(max(t0,tc))
    ! compute Rosseland mean absorption coefficient using Gails formula
!CDIR IEXPAND
    kappa = rosselandmeanopacity_new(logrho,logt)
    ! optical depth
    tau = 0.5*kappa*sigma
    ! effective optical depth
    tau_eff = 0.5*tau + 1./(3.*tau) + 1./sqrt_three
    ! temperature ratio
    t0tc = t0/tc
    ! cooling term
    qcool = qf/tau_eff * exp(4.0*logt) * (1.0-t0tc*t0tc*t0tc*t0tc) ! = Tc**4-T0**4
  END FUNCTION lambda_gray
 
  ELEMENTAL FUNCTION lambda_gammie(Eint,t_cool_inv) RESULT(Qcool)
    IMPLICIT NONE
    !---------------------------------------------------------------------!
    REAL, INTENT(IN) :: eint,t_cool_inv
    REAL :: qcool
    !---------------------------------------------------------------------!
    qcool = eint * t_cool_inv
  END FUNCTION lambda_gammie
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(sources_diskcooling), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    IF (ALLOCATED(this%ephir))  DEALLOCATE(this%ephir)
  END SUBROUTINE finalize
 
END MODULE sources_diskcooling_mod