sources_cooling.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_cooling.f90                                               #
!#                                                                           #
!# Copyright (C) 2009-2024                                                   #
!# Tobias Illenseer <tillense@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_cooling_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 sources_c_accel_mod
  USE logging_base_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER(LEN=32), PARAMETER :: source_name = "optically thin cooling"
  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(sources_base) :: sources_cooling
    TYPE(marray_base), ALLOCATABLE  :: q
    REAL      :: switchon
    REAL      :: tmin
  CONTAINS
    PROCEDURE :: initsources
    PROCEDURE :: setoutput
    PROCEDURE :: externalsources
    PROCEDURE :: calctimestep
    PROCEDURE :: updatecooling
    final :: finalize
  END TYPE
 
  !--------------------------------------------------------------------------!
  PUBLIC :: &
    ! types
    sources_cooling
  !--------------------------------------------------------------------------!
 
CONTAINS
 
  SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
    USE physics_euler_mod, ONLY : physics_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_cooling), 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
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    ! call basic initialization subroutine
    CALL this%InitSources_base(stype,source_name)
 
    ! Courant number, i.e. safety factor for numerical stability
    CALL getattr(config, "cvis", this%cvis, 0.1)
 
    ! minimum temperature
    CALL getattr(config, "Tmin", this%Tmin, 1.0e-30)
 
    ! soft switch on time scale for activating cooling (negative disables)
    CALL getattr(config, "switchon", this%switchon, -1.0)
 
    ! some sanity checks
    ! isothermal modules are excluded
    SELECT TYPE (phys => physics)
    CLASS IS(physics_euler)
      ! do nothing
    CLASS DEFAULT
      ! abort
      CALL this%Error("sources_cooling::InitSources","physics not supported")
    END SELECT
 
    ALLOCATE(this%Q,stat=err)
    IF (err.NE.0) CALL this%Error("sources_cooling::InitSources","memory allocation failed")
    this%Q = marray_base()
 
    ! set initial time < 0
    this%time = -1.0
 
    ! initialize cooling function
    this%Q%data1d(:)  = 0.0
    
    ! register ouput arrays
    CALL this%SetOutput(mesh,config,io)
  END SUBROUTINE initsources
 
 
  SUBROUTINE setoutput(this,Mesh,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_cooling), INTENT(INOUT) :: this
    CLASS(mesh_base), INTENT(IN) :: Mesh
    TYPE(dict_typ), POINTER     :: config,IO
    !------------------------------------------------------------------------!
    INTEGER              :: valwrite
    !------------------------------------------------------------------------!
    ! register cooling term for output
    CALL getattr(config, "output/Qcool", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "Qcool", &
         this%Q%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
  END SUBROUTINE setoutput
 
 
  SUBROUTINE externalsources(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_cooling), 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,time,p)
        END SELECT
      END SELECT
 
      ! energy loss due to radiation processes
      s%energy%data1d(:) = -this%Q%data1d(:)
 
    END SELECT
  END SUBROUTINE externalsources
 
 
  SUBROUTINE calctimestep(this,Mesh,Physics,Fluxes,pvar,cvar,time,dt,dtcause)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_cooling), 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(INOUT) :: dt
    INTEGER,             INTENT(OUT)   :: dtcause
    !------------------------------------------------------------------------!
    REAL              :: invdt
    !------------------------------------------------------------------------!
    ! maximum of inverse cooling timescale t_cool ~ P/Q_cool
    dt = huge(invdt)
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_euler)
      SELECT TYPE(phys => physics)
      TYPE IS (physics_euler)
        SELECT TYPE(p => pvar)
        TYPE IS (statevector_euler)
          CALL this%UpdateCooling(mesh,phys,time,p)
          invdt = maxval(abs(this%Q%data1d(:) / p%pressure%data1d(:)), &
                        mask=mesh%without_ghost_zones%mask1d(:))
          IF (invdt.GT.tiny(invdt)) dt = this%cvis / invdt
        END SELECT
      END SELECT
    END SELECT
  END SUBROUTINE calctimestep
 
 
  SUBROUTINE updatecooling(this,Mesh,Physics,time,pvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_cooling),       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,Namu,scaling
    !------------------------------------------------------------------------!
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_euler)
      ! modify Qcool during switchon phase
      IF (time.GE.0.0.AND.time.LE.this%switchon) THEN
        ! compute new scaling factor < 1 during switchon phase
        scaling = sin(0.5*pi*time/this%switchon)**2
      ELSE
        scaling = 1.0
      END IF
      ! conversion factor cs2 -> T / Kelvin
      murgamma = physics%mu /(physics%Constants%RG*physics%gamma) &
        / physics%Constants%cf_temperature
      ! conversion factor mass density -> particle density / m^3
      namu = physics%Constants%NA/physics%mu &
        / physics%Constants%cf_density * physics%Constants%cf_mass
      ! compute cooling source term ~ n^2 * Λ(T) with
      ! particle density n / m^3 and temperature T / K
      ! return value of lambda is given SI units i.e. W/m^3
      WHERE (mesh%without_ghost_zones%mask1d(:))
        this%Q%data1d(:) = scaling * physics%Constants%cf_energy &
          * (namu * p%density%data1d(:))**2 &
          * lambda(murgamma * physics%bccsound%data1d(:)**2)
      ELSEWHERE
        this%Q%data1d(:) = 0.0
      END WHERE
    END SELECT
  END SUBROUTINE updatecooling
 
 
  ELEMENTAL FUNCTION lambda(T) RESULT(L)
    IMPLICIT NONE
    REAL, INTENT(IN) :: t
    REAL :: l
    IF (t.LT.1.26e+4) THEN
       ! disable cooling for T < 1.26e4 K, i.e. set L to 0
       ! this is a very very rough approximation
       ! dont't trust this function for T < 1.26e4 !!!
       l = 0.0
    ELSE IF ((t.GE.1.26e+4).AND.(t.LT.1.5e+5)) THEN
       l = 3.7d-38 * exp(0.58*log(t))
    ELSE IF ((t.GE.1.5e5).AND.(t.LT.2.2e+7)) THEN
       l = 1.2d-31 * exp(-0.68*log(t))
    ELSE
       l = 1.2d-39 * exp(0.41*log(t))
    END IF
  END FUNCTION lambda
 
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(sources_cooling), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    DEALLOCATE(this%Q)
  END SUBROUTINE finalize
 
END MODULE sources_cooling_mod