sources_diskcooling.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_diskcooling.f90                                           #
!#                                                                           #
!# Copyright (C) 2011-2019                                                   #
!# 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 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
  INTEGER, PARAMETER :: gray   = 1
  INTEGER, PARAMETER :: gammie = 2
  INTEGER, PARAMETER :: gammie_sb = 3
  !--------------------------------------------------------------------------!
  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_base) :: sources_diskcooling
    CHARACTER(LEN=32) :: source_name = "thin accretion disk cooling"
    CLASS(logging_base), ALLOCATABLE :: cooling
    REAL                            :: b_cool
    REAL                            :: b_start
    REAL                            :: b_final
    REAL                            :: t_start
    REAL                            :: dt_bdec
    TYPE(marray_base) :: qcool
    TYPE(marray_base) :: ephir
    REAL                            :: t_0
    REAL                            :: rho_0
  CONTAINS
    PROCEDURE :: initsources_diskcooling
    PROCEDURE :: infosources
    PROCEDURE :: setoutput
    PROCEDURE :: externalsources_single
    PROCEDURE :: calctimestep_single
    PROCEDURE :: updatecooling
    PROCEDURE :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       sources_diskcooling, &
       ! constants
       gray, gammie, gammie_sb
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initsources_diskcooling(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) :: 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
    INTEGER :: cooling_func
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    CALL this%InitLogging(stype,this%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("InitSources_diskcooling","physics not supported")
    END SELECT

    SELECT CASE(cooling_func)
    CASE(gammie)
      ! check for geometry
      SELECT TYPE (mgeo => mesh%Geometry)
      CLASS IS(geometry_cylindrical)
        ! do nothing
      CLASS DEFAULT
        CALL this%Error("InitSources_diskcooling","geometry not supported")
      END SELECT
    CASE(gammie_sb)
      SELECT TYPE (mgeo => mesh%Geometry)
      TYPE IS(geometry_cartesian)
        ! do nothing
      CLASS DEFAULT
        CALL this%Error("InitSources_diskcooling","geometry not supported")
      END SELECT
    END SELECT

    ! check for dimensions
    IF (physics%VDIM.GT.2) THEN
      CALL this%Error("InitSources_diskcooling","only 2D simulations supported")
    END IF

    ! get cooling method
    CALL getattr(config,"method",cooling_func)
    ALLOCATE(logging_base::this%cooling)
    SELECT CASE(cooling_func)
    CASE(gray)
      IF (physics%constants%GetType().NE.si) &
         CALL this%Error("InitSources_diskcooling","only SI units supported for gray cooling")
      CALL this%cooling%InitLogging(gray,"gray cooling")

      this%Qcool = marray_base()

    CASE(gammie)
      CALL this%cooling%InitLogging(gammie,"Gammie cooling")

      this%Qcool = marray_base()
      this%ephir = marray_base(2)

    CASE(gammie_sb)
      CALL this%cooling%InitLogging(gammie_sb,"Gammie cooling (SB)")

      this%Qcool = marray_base()

    CASE DEFAULT
       CALL this%Error("UpdateCooling","Cooling function not supported!")
    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 arrays
    this%Qcool%data1d(:)  = 0.0
    SELECT CASE(cooling_func)
    CASE(gray)
      ! do nothing
    CASE(gammie)
      this%ephir%data4d(:,:,:,1) = -mesh%posvec%bcenter(:,:,:,2)/mesh%radius%bcenter(:,:,:)**2
      this%ephir%data4d(:,:,:,2) = mesh%posvec%bcenter(:,:,:,1)/mesh%radius%bcenter(:,:,:)**2
    CASE(gammie_sb)
      ! do nothing
    END SELECT

    !initialise output
    CALL this%SetOutput(mesh,config,io)
    CALL this%InitSources(mesh,fluxes,physics,config,io)

  END SUBROUTINE initsources_diskcooling


  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,'(ES8.2)') this%T_0
       CALL this%Info("            min. temperature:  " // trim(param_str))
       WRITE (param_str,'(ES8.2)') this%rho_0
       CALL this%Info("            minimum density:   " // trim(param_str))
    CASE(gammie)
       WRITE (param_str,'(ES8.2)') this%b_cool
       CALL this%Info("            cooling parameter: " // trim(param_str))
    CASE(gammie_sb)
       WRITE (param_str,'(ES8.2)') this%b_cool
       CALL this%Info("            cooling parameter: " // trim(param_str))
       IF(this%dt_bdec.GE.0.0) THEN
        WRITE (param_str,'(ES8.2)') this%b_start
        CALL this%Info("            initial cooling parameter: " // trim(param_str))
        WRITE (param_str,'(ES8.2)') this%b_final
        CALL this%Info("            final cooling parameter: " // trim(param_str))
        WRITE (param_str,'(ES8.2)') this%t_start
        CALL this%Info("            starting b_dec time: " // trim(param_str))
        WRITE (param_str,'(ES8.2)') this%dt_bdec
        CALL this%Info("            operating time: " // trim(param_str))
       END IF
    END SELECT
  END SUBROUTINE infosources


  SUBROUTINE setoutput(this,Mesh,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_diskcooling)  :: this
    CLASS(mesh_base), INTENT(IN) :: Mesh
    TYPE(Dict_TYP), POINTER     :: config,IO
    !------------------------------------------------------------------------!
    INTEGER              :: valwrite
    !------------------------------------------------------------------------!

    !cooling energy source term
    CALL getattr(config, "output/Qcool", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "Qcool", &
         this%Qcool%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))

  END SUBROUTINE setoutput


  SUBROUTINE externalsources_single(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    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_base),INTENT(INOUT)   :: Physics
    CLASS(sources_base),INTENT(INOUT)   :: Sources
    CLASS(fluxes_base),INTENT(IN)       :: Fluxes
    REAL,INTENT(IN)                     :: time, dt
    CLASS(marray_compound),INTENT(INOUT):: pvar,cvar,sterm
    !------------------------------------------------------------------------!
    SELECT TYPE(s => sterm)
    TYPE IS (statevector_euler)
      s%density%data1d(:) = 0.0
      s%momentum%data1d(:) = 0.0

      SELECT TYPE(phys => physics)
      TYPE IS (physics_euler)
        SELECT TYPE(p => pvar)
        TYPE IS (statevector_euler)
          CALL this%UpdateCooling(mesh,phys,sources,time,p)
        END SELECT
      END SELECT

      ! energy loss due to radiation processes
      s%energy%data1d(:) = -this%Qcool%data1d(:)

    END SELECT
  END SUBROUTINE externalsources_single


  SUBROUTINE calctimestep_single(this,Mesh,Physics,Fluxes,pvar,cvar,time,dt)
    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_base), INTENT(INOUT) :: Physics
    CLASS(fluxes_base),  INTENT(IN)    :: Fluxes
    CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
    REAL,                INTENT(IN)       :: time
    REAL,                INTENT(OUT)      :: dt
    !------------------------------------------------------------------------!
    REAL              :: invdt
    !------------------------------------------------------------------------!
    ! maximum of inverse cooling timescale t_cool ~ P/Q_cool
    dt = huge(invdt)
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_euler)
      invdt = maxval(abs(this%Qcool%data1d(:) / p%pressure%data1d(:)), &
                     mask=mesh%without_ghost_zones%mask1d(:))
      IF (invdt.GT.tiny(invdt)) dt = this%cvis / invdt
    END SELECT
  END SUBROUTINE calctimestep_single

  SUBROUTINE updatecooling(this,Mesh,Physics,Sources,time,pvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_diskcooling)         :: this
    CLASS(mesh_base),    INTENT(IN)    :: Mesh
    CLASS(physics_euler),INTENT(INOUT) :: Physics
    CLASS(sources_base), INTENT(INOUT) :: Sources
    CLASS(statevector_euler),INTENT(INOUT):: pvar
    REAL,                INTENT(IN)    :: time
    !------------------------------------------------------------------------!
    REAL              :: muRgamma,Qfactor
    !------------------------------------------------------------------------!
    ! calculate value for beta
    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

    ! energy loss due to radiation processes
    SELECT CASE(this%cooling%GetType())
       CASE(gray)
          ! some constants
          murgamma = physics%mu/(physics%Constants%RG*physics%gamma)
          qfactor  = 8./3.*physics%Constants%SB
          ! compute gray cooling term
          this%Qcool%data1d(:) = lambda_gray(pvar%data2d(:,physics%DENSITY),sources%height%data1d(:), &
               murgamma*physics%bccsound%data1d(:)*physics%bccsound%data1d(:), &
               this%rho_0,this%T_0,qfactor)
       CASE(gammie)
          ! compute Gammie cooling term with
          ! t_cool = b_cool / Omega
          ! and Omega = ephi*v / r
          this%Qcool%data3d(:,:,:) = lambda_gammie(pvar%data4d(:,:,:,physics%PRESSURE) / (physics%gamma-1.), &
              abs((this%ephir%data4d(:,:,:,1)*pvar%data4d(:,:,:,physics%XVELOCITY) &
                  +this%ephir%data4d(:,:,:,2)*(pvar%data4d(:,:,:,physics%YVELOCITY)&
                    +mesh%OMEGA*mesh%radius%bcenter(:,:,:)))) / this%b_cool)
       CASE(gammie_sb)
          ! in sb t_cool = b_cool
          this%Qcool%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
    !------------------------------------------------------------------------!
    CLASS(sources_diskcooling), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    SELECT CASE(this%cooling%GetType())
    CASE(gray)
      CALL this%Qcool%Destroy()
    CASE(gammie)
      CALL this%Qcool%Destroy()
      CALL this%ephir%Destroy()
    CASE(gammie_sb)
      CALL this%Qcool%Destroy()
    END SELECT
  END SUBROUTINE finalize

END MODULE sources_diskcooling_mod