ttauridisk.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 2D hydrodynamical simulation program                             #
!# module: ttauridisk.f90                                                    #
!#                                                                           #
!# Copyright (C) 2006-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.                 #
!#                                                                           #
!#############################################################################
 
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
PROGRAM ttauri
  USE fosite_mod
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPIF_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER :: gn      = 6.6742d-11     ! Newtons grav. constant       !
  REAL, PARAMETER :: cc      = 2.99792458d+8  ! speed of light               !
  REAL, PARAMETER :: msun    = 1.989d+30      ! solar mass [kg]              !
  REAL, PARAMETER :: au      = 1.49597870691e+11 ! astronomical unit [m]     !
  REAL, PARAMETER :: year    = 3.15576e+7           ! Julian year [sec]      !
  REAL, PARAMETER :: parsec  = 0.5/pi*1.296e+6 * au ! parsec [m]             !
  ! simulation parameters
  REAL, PARAMETER :: tsim    = 1.0e+4 * year
  REAL, PARAMETER :: m_bh    = 0.5e+0 * msun  ! mass of central point mass   !
  REAL, PARAMETER :: m_disk  = 1.0e-1 * m_bh  ! accretion disk mass          !
  REAL, PARAMETER :: rs = 2 * gn * m_bh / cc**2 ! Schwarzschild radius       !
  REAL, PARAMETER :: lscale  = au             ! typical length scale         !
  REAL, PARAMETER :: etasgs = 1.0e-6
  REAL, PARAMETER :: gamma = 1.4
  ! computational domain
  REAL, PARAMETER :: rmin = 2.0d+0 * lscale
  REAL, PARAMETER :: rmax = 3.0d+2 * lscale
  REAL, PARAMETER :: zmin = -1.3d+2 * lscale
  REAL, PARAMETER :: zmax = 1.3d+2 * lscale
  REAL, PARAMETER :: ztan = 1.0e+2 * lscale   ! scale for tancyl. coordinates!
  ! initial condition:
  !   hydrostatic equilibrium for slightly non-keplerian flow
  !   deviation from Keplerian: v_phi(s,z) = v_kep(s,z) / lambda(z)
  !                             lambda(z)  = 1 + A0/(1+(z/Z0)**2)
  !   pressure in equatorial plane: Pc(s) = P0 * EXP(-s/S0)
  REAL, PARAMETER :: z0      = 3.0e+02 * lscale ! vertical scale             !
  REAL, PARAMETER :: a0      = 1.0e-01        ! velocity deviation           !
  REAL, PARAMETER :: s0      = 1.5e+02 * lscale ! radial scale               !
  REAL, PARAMETER :: b0      = -1.0 ! <1 !    ! slope for density power law  !
  REAL, PARAMETER :: rho_inf = 1.0e-15        ! density at infinity          !
  REAL, PARAMETER :: p_inf   = 1.4e-10        ! pressure at infinity         !
  INTEGER, PARAMETER :: problem = euler
  ! resolution
  INTEGER, PARAMETER :: rres = 100
  INTEGER, PARAMETER :: phires = 1
  INTEGER, PARAMETER :: zres = 128
  ! output file parameter
  INTEGER, PARAMETER :: onum = 1000                  ! number of outputs      !
  CHARACTER(LEN=256), PARAMETER :: odir &! output directory
                                 = "./"
  CHARACTER(LEN=256), PARAMETER &                   ! output data file name  !
                     :: ofname = 'Ttauri'
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE  :: sim
  !--------------------------------------------------------------------------!
 
  ALLOCATE(sim)
 
  CALL sim%InitFosite()
 
  CALL makeconfig(sim%config)
 
  CALL sim%Setup()
 
  ! set initial condition
  CALL initdata(sim%Mesh, sim%Physics, sim%Fluxes, sim%Timedisc, sim%Sources)
 
  CALL sim%Run()
 
  CALL sim%Finalize()
 
CONTAINS
 
   SUBROUTINE makeconfig(config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(dict_typ),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, logfile, sources, &
                               timedisc,fluxes,stemp
    !------------------------------------------------------------------------!
 
   ! mesh settings
    mesh => dict(&
             "meshtype" / midpoint, &
             "geometry" / cylindrical, &
             "inum" / rres, &        ! resolution in x and            !
             "jnum" / phires, &        ! resolution in y and            !
             "knum" / zres, &        !   z direction                  !             
             "xmin" / rmin, &
             "xmax" / rmax, &
             "ymin" / (0.0), &
             "ymax" / (2.*pi), &
             "zmin" / zmin, &
             "zmax" / zmax)
 
    ! physics settings
    physics => dict(&
         "problem" / problem, &
         "gamma"   / gamma, &          ! ratio of specific heats        !
         "mu"      / 6.02e-4, &      ! mean molecular weight          !
         "dpmax"   / 1.0e+3)         ! for advanced time step control !
 
    ! numerical scheme for flux calculation
    fluxes => dict(&
         "fluxtype"  / kt, &
         "order"     / linear, &
         "variables" / primitive, &   ! vars. to use for reconstruction!
         "limiter"   / minmod, &    ! one of: minmod, monocent,...   !
         "theta"     / 1.2)           ! optional parameter for limiter !
 
     ! boundary conditions
    boundary => dict( &
!!$         "western"  / REFLECTING, &
!!$         +"eastern"  / REFLECTING, &
!!$         +"western"  / CUSTOM, &
!!$         +"eastern"  / CUSTOM, &
!!$         +"western"  / NO_GRADIENTS, &
!!$         +"eastern"  / NO_GRADIENTS, &
         "western"  / custom, &
         "eastern"  / farfield, &
!!$         +"southern" / AXIS, &
!!$         +"southern" / REFLECTING, &
!!$         +"northern" / REFLECTING)
         "southern" / periodic, &
!!$         +"northern" / CUSTOM)
!!$         +"southern" / NO_GRADIENTS, &
!!$         +"northern" / NO_GRADIENTS)
!!$         +"southern" / FARFIELD, &
         "northern" / periodic, &
     !    "bottomer" / NO_GRADIENTS, &
         "bottomer" / farfield, &
         "topper" / farfield)
!!$         +"northern" / REFLECTING)
 
    ! source term due to a point mass
    stemp => dict("stype" / gravity, &
           "pmass/gtype"  / pointmass, &            ! grav. accel. of a point mass     !
           "pmass/mass" / m_bh)               ! mass of the accreting object[kg] !
 
    NULLIFY(sources)
    CALL setattr(sources, "grav", stemp)
 
 
    ! account for energy losses due to radiative cooling
!!$    stemp => Dict(, &
!!$         "stype" / COOLING,, &
!!$         "cvis"  / 0.9)                     ! CFL number for cooling time    !
!!$    CALL SetAttr(sources, "cooling", stemp)
 
!    ! viscosity source term
!    stemp => Dict("stype"   / VISCOSITY, &
!         "cvis"     / 0.6, &
!!!$         +"vismodel" / ALPHA, &
!!!$         +"dynconst" / 0.1)
!         "vismodel" / BETA, &
!         "dynconst" / 0.001, &
!           "output/stress" / 1, &
!!           "output/kinvis" / 1, &
!           "output/dynvis" / 1)
 
    ! time discretization settings
    timedisc => dict(&
         "method"   / modified_euler, &
         "order"    / 3, &
         "cfl"      / 0.4, &
         "stoptime" / tsim, &
         "dtlimit"  / 1.0e-3, &
         "maxiter"  / 1000000000, &
         "output/external_sources" / 1, &
         "output/geometrical_sources" / 1, &
         "output/fluxes" / 0)
 
!    ! initialize log input/output
!    logfile => Dict(, &
!         "fileformat" / BINARY, &
!         +"filename"   / TRIM(lfname), &
!         +"dtwall"     / 3600, &
!         +"filecycles" / 1)
 
    ! initialize data input/output
    datafile => dict("fileformat" / vtk, &
          "unit" / 5555, &
          "filename"   / (trim(odir) // trim(ofname)), &
          "stoptime"   / tsim, &
          "count"      / onum, &
          "filecycles" / (onum+1))
 
    config => dict("mesh" / mesh, &
             "physics"  / physics, &
             "boundary" / boundary, &
             "fluxes"   / fluxes, &
             "sources"  / sources, &
             "timedisc" / timedisc, &
!             "logfile"  / logfile, &
             "datafile" / datafile)
 
  END SUBROUTINE makeconfig
 
  SUBROUTINE initdata(Mesh,Physics,Fluxes,Timedisc,Sources)
    USE sources_base_mod, ONLY : sources_base
    USE sources_gravity_mod, ONLY : sources_gravity
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base)    ,INTENT(IN) :: Mesh
    CLASS(physics_base) ,INTENT(INOUT) :: Physics
    CLASS(fluxes_base)  ,INTENT(INOUT) :: Fluxes
    CLASS(timedisc_base),INTENT(INOUT):: Timedisc
    CLASS(sources_list),ALLOCATABLE,INTENT(INOUT) :: Sources
    !------------------------------------------------------------------------!
    ! Local variable declaration
    CLASS(sources_base), POINTER :: sp
    CLASS(sources_gravity), POINTER :: gp => null()
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,2) :: dv
    REAL              :: s,z,tau,lambda,tau1,tau2,dPs,dPz,s1,s2,z1,z2,P0,rho0,cs_inf
    REAL              :: mdisk,am
#ifdef PARALLEL
    INTEGER           :: ierror
    REAL              :: mdisk_local,am_local
#endif
    INTEGER           :: i,j,k,n
    CHARACTER(LEN=20) :: mdisk_str,am_str
    !------------------------------------------------------------------------!
    ! get gravitational acceleration
    IF (ALLOCATED(sources)) &
      sp => sources%GetSourcesPointer(gravity)
    IF (.NOT.ASSOCIATED(sp)) &
      CALL physics%Error("ttauridisk::InitData","no gravity term initialized")
 
    ! initial condition (computational domain)
    ! the total disk mass is proportional to the central pressure;
    ! to determine the necessary central pressure for a given
    ! disk mass, the initial condition is computed twice
    p0 = 1.0e+10
    rho0 = 1.0e-09
    cs_inf = sqrt(gamma*p_inf/rho_inf)
    DO n=1,2
      DO k=mesh%KGMIN,mesh%KGMAX
       DO j=mesh%JGMIN,mesh%JGMAX
          DO i=mesh%IGMIN,mesh%IGMAX
             s = abs(mesh%bccart(i,j,k,1))
             z = mesh%bccart(i,j,k,3)
if (s .eq. 0.0 .or. z .eq. 0.0) print *,s,z,i,j
 
             tau = func_tau(s,z)
             lambda = func_lambda(z)
             timedisc%pvar%data4d(i,j,k,physics%XVELOCITY) = 0.0
             timedisc%pvar%data4d(i,j,k,physics%ZVELOCITY) = 0.0
             timedisc%pvar%data4d(i,j,k,physics%PRESSURE) = pc(tau,rho0,cs_inf)
             z1 = mesh%cart%faces(i,j,k,5,3) ! western cell faces z coordinate
             z2 = mesh%cart%faces(i,j,k,6,3) ! eastern cell faces z coordinate
             tau1 = func_tau(s,z1)
             tau2 = func_tau(s,z2)
             dpz = (pc(tau2,rho0,cs_inf)-pc(tau1,rho0,cs_inf))/mesh%dlz%data3d(i,j,k)
             IF (abs(z).GT.0.0) THEN
!!$                Timedisc%pvar(i,j,Physics%DENSITY) = rho0 * rho(s,z,tau,lambda)
                timedisc%pvar%data4d(i,j,k,physics%DENSITY)   = dpz / gz(s,z)
             ELSE
                timedisc%pvar%data4d(i,j,k,physics%DENSITY)   = rho0 * (s/s0)**b0
             END IF
             s1 = abs(mesh%cart%faces(i,j,k,1,1)) !vorher (3,1) ! southern cell faces s coordinate
             s2 = abs(mesh%cart%faces(i,j,k,2,1)) !vorher (4,1) ! northern cell faces s coordinate
             tau1 = func_tau(s1,z)
             tau2 = func_tau(s2,z)
             dps = (pc(tau2,rho0,cs_inf)-pc(tau1,rho0,cs_inf))/mesh%dlx%data3d(i,j,k)
!             Timedisc%pvar%data4d(i,j,k,Physics%YVELOCITY) = SQRT(MAX(0.0, &
!                  s*(dPs/Timedisc%pvar%data4d(i,j,k,Physics%DENSITY) - gs(s,z))))
             timedisc%pvar%data4d(i,j,k,physics%YVELOCITY) = sqrt(max(0.0,-s*gs(s,z) / lambda))
          END DO
        END DO
       END DO
       ! compute disk mass
       mdisk = sum(timedisc%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%DENSITY) &
            * mesh%volume%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
       ! compute angular momentum
       am = sum(mesh%hy%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) &
            * timedisc%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%YVELOCITY) &
            * timedisc%pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%DENSITY) &
            * mesh%volume%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
 
#ifdef PARALLEL
       mdisk_local = mdisk 
       CALL mpi_allreduce(mdisk_local,mdisk,1,default_mpi_real,mpi_sum, &
            mesh%comm_cart,ierror)
       am_local = am 
       CALL mpi_allreduce(am_local,am,1,default_mpi_real,mpi_sum, &
            mesh%comm_cart,ierror)
#endif
!!$       ! adjust central pressure
!!$       P0 = M_DISK / mdisk * P0
       ! adjust density at (S0,Z=0)
       rho0 = m_disk / mdisk * rho0
    END DO
      ! 2. initialize gravitational acceleration
!      CALL gp%UpdateGravity(Mesh,Physics,Fluxes,Timedisc%pvar,0.0,0.0)
!       Timedisc%pvar%data4d(:,:,:,2:2+Physics%VDIM) = &
!                Timedisc%GetCentrifugalVelocity(Mesh,Physics,Fluxes,Sources,(/0.,0.,1./),gp%accel%data4d)
 
 
    ! add velocity perturbations
!!$    CALL RANDOM_SEED
!!$    CALL RANDOM_NUMBER(dv)
!!$    Timedisc%pvar(:,:,Physics%XVELOCITY) = Timedisc%pvar(:,:,Physics%XVELOCITY) &
!!$         + (dv(:,:,1)-0.5)*2.0E+2
!!$    Timedisc%pvar(:,:,Physics%YVELOCITY) = Timedisc%pvar(:,:,Physics%YVELOCITY) &
!!$         + (dv(:,:,2)-0.5)*2.0E+2
 
      SELECT TYPE (bbottom => timedisc%Boundary%Boundary(bottom)%p)
      CLASS IS (boundary_farfield)
       DO k=1,mesh%GKNUM
         DO j=mesh%JMIN,mesh%JMAX
            DO i=mesh%IMIN,mesh%IMAX
               bbottom%data(i,j,k,:) = timedisc%pvar%data4d(i,j,mesh%KMIN-k,:)
            END DO
         END DO
       END DO
     END SELECT
 
      SELECT TYPE (btop => timedisc%Boundary%Boundary(top)%p)
      CLASS IS (boundary_farfield)
       DO k=1,mesh%GKNUM
         DO j=mesh%JMIN,mesh%JMAX
            DO i=mesh%IMIN,mesh%IMAX
               btop%data(i,j,k,:) = timedisc%pvar%data4d(i,j,mesh%KMAX+k,:)
            END DO
         END DO
       END DO
     END SELECT
 
      SELECT TYPE (beast => timedisc%Boundary%Boundary(east)%p)
      CLASS IS (boundary_farfield)
       DO k=mesh%KMIN,mesh%KMAX
         DO j=mesh%JMIN,mesh%JMAX
            DO i=1,mesh%GINUM
               beast%data(i,j,k,:) = timedisc%pvar%data4d(mesh%IMAX+i,j,k,:)
            END DO
         END DO
       END DO
     END SELECT
 
      ! western
      SELECT TYPE (bwest => timedisc%Boundary%Boundary(west)%p)
      CLASS IS(boundary_custom)
        CALL bwest%SetCustomBoundaries(mesh,physics, &
          (/custom_nograd,custom_outflow,custom_kepler,custom_nograd,custom_nograd/))
      END SELECT
 
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
    CALL timedisc%Boundary%CenterBoundary(mesh,physics,0.0,timedisc%pvar, &
         timedisc%cvar)
 
    ! print some information
    WRITE (mdisk_str, '(ES10.2)') mdisk/msun
    WRITE (am_str, '(ES10.2)') am
    CALL mesh%Info(" DATA-----> initial condition: non Keplerian flow")
    CALL mesh%Info("            disk mass:         " // trim(mdisk_str) // " M_sun")
    CALL mesh%Info("            angular momentum:  " // trim(am_str) // " kg/m^2/s")
 
  END SUBROUTINE initdata
    
    ELEMENTAL FUNCTION vkepler(sc,zc) RESULT(res)
      IMPLICIT NONE
      REAL,INTENT(IN) :: sc,zc
      REAL :: res
      res = sqrt(gn*m_bh/sc) * (1.0d0+(zc/sc)**2)**(-0.75)
    END FUNCTION vkepler
 
    ELEMENTAL FUNCTION gs(sc,zc) RESULT(res)
      IMPLICIT NONE
      REAL,INTENT(IN) :: sc,zc
      REAL :: res
      res = -(gn*m_bh/sc**2) / sqrt(1.0d0+(zc/sc)**2)**3
    END FUNCTION gs
    
    ELEMENTAL FUNCTION gz(sc,zc) RESULT(res)
      IMPLICIT NONE
      REAL,INTENT(IN) :: sc,zc
      REAL :: res
      res = -(gn*m_bh/sc**2) * (zc/sc) / sqrt(1.0d0+(zc/sc)**2)**3
    END FUNCTION gz
    
    ELEMENTAL FUNCTION func_lambda(zc) RESULT(res)
      IMPLICIT NONE
      REAL,INTENT(IN) :: zc
      REAL :: res
      res = 1.0d0 + a0/(1.0d0+(zc/z0)**2)
    END FUNCTION func_lambda
 
    ELEMENTAL FUNCTION func_tau(sc,zc) RESULT(res)
      IMPLICIT NONE
      REAL,INTENT(IN) :: sc,zc
      REAL :: res
      res = sc*sqrt(1.0d0 + (zc/sc)**2*(1.0d0 + (1.0d0+0.5d0*(zc/z0)**2)/a0))
    END FUNCTION func_tau
 
    ELEMENTAL FUNCTION pc(tau,rho_s0,cs_inf) RESULT(res)
      IMPLICIT NONE
      REAL,INTENT(IN) :: tau,rho_s0,cs_inf
      REAL :: res
!      res = P0 * EXP(-tau/S0)
      res = p_inf * (1.0 + a0/(1.0+a0) * 1.0/(1.0-b0) * 0.5 * gamma * &
           (cc/cs_inf)**2 * rho_s0/rho_inf * (rs/s0) * (tau/s0)**(b0-1.0))
    END FUNCTION pc
 
!!$    ELEMENTAL FUNCTION dPc(tau) RESULT(res)
!!$      IMPLICIT NONE
!!$      REAL,INTENT(IN) :: tau
!!$      REAL :: res
!!$      res = -1.0D0/S0 * Pc(tau)
!!$    END FUNCTION dPc
 
!    ELEMENTAL FUNCTION rho(sc,zc,tau,lambda) RESULT(res)
!      IMPLICIT NONE
!      REAL,INTENT(IN) :: sc,zc,tau,lambda
!      REAL :: res
!      res = A0/(1.0+A0) * lambda/(lambda-1.0) * sqrt(1.0+(zc/sc)**2)**3 &
!           * (tau/S0)**B0 * (s/tau)**3
!    END FUNCTION rho
 
 
END PROGRAM ttauri