pringle3d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: pringle3d.f90                                                     #
!#                                                                           #
!# Copyright (C) 2008-2019                                                   #
!# Bjoern Sperling  <sperling@astrophysik.uni-kiel.de>                       #
!# 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 pringle3d
  USE fosite_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 5.0e-3     ! simulation time [TVIS] see below
                                             ! should be < 1 ~ viscous time scale
  ! these are the basic parameters; for stable solutions one requires
  ! MA >> 1 and RE > MA
  REAL, PARAMETER    :: re      = 1.0e+4     ! Reynolds number (at R=1)
  REAL, PARAMETER    :: ma      = 1.0e+2     ! Mach number (at R=1)
  REAL, PARAMETER    :: mdisk   = 1.0e-2     ! disk mass << 1 = central mass
  ! lower limit for nitial density
  REAL, PARAMETER    :: rhomin  = 1.0e-8    ! minimal initial density
  ! viscosity prescription
  INTEGER, PARAMETER :: vistype = beta
!   INTEGER, PARAMETER :: VISTYPE = POWERLAW
  REAL, PARAMETER    :: pl_exp  = 0.0        ! exponent for power law viscosity
  REAL, PARAMETER    :: tvis    = 4./3.*re   ! viscous time scale
  REAL, PARAMETER    :: tinit   = 1.0e-3     ! time for initial condition [TVIS]
  ! mesh settings
!   INTEGER, PARAMETER :: MGEO = CYLINDRICAL
!   INTEGER, PARAMETER :: MGEO = LOGCYLINDRICAL
  INTEGER, PARAMETER :: mgeo = spherical !!! ATTENTION: not applicable in 1D
  INTEGER, PARAMETER :: xres = 200        ! x-resolution
  INTEGER, PARAMETER :: yres = 16           ! y-resolution
  INTEGER, PARAMETER :: zres = 6           ! z-resolution
  REAL, PARAMETER    :: rmin = 0.1        ! min radius of comp. domain
  REAL, PARAMETER    :: rmax = 2.0         ! max radius of comp. domain
  REAL, PARAMETER    :: gpar = 1.0         ! geometry scaling parameter
  ! output parameters
  INTEGER, PARAMETER :: onum = 100          ! number of output data sets
  CHARACTER(LEN=256), PARAMETER &          ! output data dir
                     :: odir = './'
  CHARACTER(LEN=256), PARAMETER &          ! output data file name
                     :: ofname = 'pringle3d'
  ! derived parameters
  REAL, PARAMETER    :: csiso = 1./ma      ! isothermal speed of sound
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE   :: sim
  !--------------------------------------------------------------------------!

  ALLOCATE(sim)
  CALL sim%InitFosite()
  CALL makeconfig(sim, sim%config)
  CALL sim%Setup()
  CALL initdata(sim,sim%Mesh, sim%Physics, sim%Timedisc)
  CALL sim%Run()
  CALL sim%Finalize()
  DEALLOCATE(sim)

CONTAINS

  SUBROUTINE makeconfig(Sim, config)
    USE functions, ONLY : asinh
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite)          :: Sim
    TYPE(Dict_TYP),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: bc(6)
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
                               grav, vis, timedisc, fluxes, sources
    REAL              :: x1,x2,y1,y2,z1,z2,cvis
    !------------------------------------------------------------------------!
    INTENT(INOUT)     :: sim
    !------------------------------------------------------------------------!
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(cylindrical)
      x1 = rmin
      x2 = rmax
      y1 = -pi
      y2 = pi
!       z1 = 0.0
!       z2 = 0.0
      z1 = -0.1
      z2 = 0.1
!       bc(WEST)  = AXIS
      bc(west)  = no_gradients
      bc(east)  = no_gradients
!       bc(WEST)  = ABSORBING
!       bc(EAST)  = ABSORBING
!       bc(WEST)  = CUSTOM
!       bc(EAST)  = CUSTOM
      bc(south) = periodic
      bc(north) = periodic
!       bc(BOTTOM) = ABSORBING
!       bc(TOP)    = ABSORBING
      bc(bottom) = reflecting
      bc(top)    = reflecting
      !!! ATTENTION:  empirical formula found for standard parameters
      !!!             seems to work for RMIN=0.05 and XRES <= 1000
      cvis = (xres/1207.0)**1.9
    CASE(logcylindrical)
      x1 = log(rmin/gpar)
      x2 = log(rmax/gpar)
      y1 = -pi
      y2 = pi
      z1 = 0.0
      z2 = 0.0
      bc(west)  = no_gradients
      bc(east)  = no_gradients
      bc(south) = periodic
      bc(north) = periodic
      bc(bottom) = no_gradients
      bc(top)    = no_gradients
      cvis = 0.1
    CASE(spherical)
      x1 = rmin
      x2 = rmax
      y1 = pi/2 - 0.02*pi
      y2 = pi/2 + 0.02*pi
      z1 = -pi
      z2 = pi
!       bc(WEST)  = ABSORBING
      bc(east)  = absorbing
!       bc(WEST)  = NO_GRADIENTS
!       bc(EAST)  = NO_GRADIENTS
      bc(west)  = custom
!       bc(EAST)  = CUSTOM
!       bc(SOUTH) = NO_GRADIENTS
!       bc(NORTH) = NO_GRADIENTS
      bc(south) = absorbing
      bc(north) = absorbing
      bc(bottom) = periodic
      bc(top)    = periodic
      !!! ATTENTION:  empirical formula found for standard parameters
      !!!             seems to work for RMIN=0.05 and XRES <= 1000
      cvis = (xres/1207.0)**1.9
    CASE DEFAULT
      CALL sim%Error("pringle::MakeConfig","mesh geometry not supported for Pringle disk")
    END SELECT

    ! mesh settings
    mesh => dict( &
           "meshtype"   / midpoint, &
           "geometry"   / mgeo, &
           "inum"       / xres, &
           "jnum"       / yres, &
           "knum"       / zres, &
           "xmin"       / x1, &
           "xmax"       / x2, &
           "ymin"       / y1, &
           "ymax"       / y2, &
           "zmin"       / z1, &
           "zmax"       / z2, &
!            "use_fargo"  / 1, &
!            "fargo"      / 2, &
           "gparam"     / gpar)

    ! boundary conditions
    boundary => dict( &
            "western"   / bc(west), &
            "eastern"   / bc(east), &
            "southern"  / bc(south), &
            "northern"  / bc(north), &
            "bottomer"  / bc(bottom), &
            "topper"    / bc(top))


    ! physics settings
    physics => dict( &
            "problem"   / euler_isotherm, &
            "units"     / geometrical, &
            "cs"        / csiso)

    ! flux calculation and reconstruction method
    fluxes => dict( &
            "order"     / linear, &
            "fluxtype"  / kt, &
            "variables" / primitive, &
            "limiter"   / vanleer, &
            "theta"     / 1.2)

    ! source term due to viscosity
    vis => dict( &
            "stype"     / viscosity, &
            "vismodel"  / vistype, &
            "dynconst"  / (1./re), &
            "exponent"  / pl_exp, &
            "cvis"      / cvis)

    ! source term due to a point mass
    grav => dict( &
            "stype"     / gravity, &
            "pmass/gtype"/ pointmass, &    ! grav. accel. of a point mass
            "pmass/mass" / 1.0, &          ! non-dim. mass of the accreting object
            "pmass/outbound" / 0)          ! disable accretion

    ! combine all source terms
    sources => dict( &
            "viscosity" / vis, &
            "grav"      / grav)

    ! time discretization settings
    timedisc => dict( &
            "method"    / ssprk, &
            "stoptime"  / (tsim * tvis), &
            "cfl"       / 0.4, &
            "dtlimit"   / (1e-10 * tvis), &
!             "rhstype"   / 1, &
            "maxiter"   / 100000000, &
            "tol_rel"   / 0.01, &
            "tol_abs"   / (/1e-40,1e-5,1e-8,1e-5/))

!     ! enable output of analytical solution
! #ifdef HAVE_FGSL
!     IF (MGEO.EQ.CYLINDRICAL) &
!       CALL SetAttr(timedisc,"output/solution",1)
! #endif

    ! initialize data input/output
    datafile => dict( &
            "fileformat"/ vtk, &
            "filename"  / (trim(odir) // trim(ofname)), &
            "count"     / onum)

    config => dict( &
            "mesh"      / mesh, &
            "physics"   / physics, &
            "boundary"  / boundary, &
            "fluxes"    / fluxes, &
            "sources"   / sources, &
            "timedisc"  / timedisc, &
            "datafile"  / datafile)
  END SUBROUTINE makeconfig

  SUBROUTINE initdata(Sim,Mesh,Physics,Timedisc)
    USE geometry_generic_mod
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite),   INTENT(INOUT) :: Sim
    CLASS(physics_base),  INTENT(IN)    :: Physics
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    !------------------------------------------------------------------------!
    ! Local variable declaration
    CLASS(sources_base), POINTER :: sp
    CLASS(sources_gravity), POINTER :: gp
    CLASS(geometry_base), ALLOCATABLE :: geo_cyl
    TYPE(Dict_TYP), POINTER :: geo_config
    INTEGER           :: i,j,k
#ifdef PARALLEL
    INTEGER           :: ierror
#endif
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) &
                      :: bccyl
    REAL              :: r,z,vr,vphi,sden,mu,height,mass
    CHARACTER(LEN=64) :: value
    !------------------------------------------------------------------------!
    geo_config => dict("geometry" / cylindrical)
    CALL new_geometry(geo_cyl,geo_config)
    CALL deletedict(geo_config)

    ! compute cylindrical coordinates for each cell bary center
    CALL geo_cyl%Convert2Curvilinear(mesh%bccart,bccyl)

    ! custom boundary conditions if requested
    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/))
    END SELECT
    SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
    CLASS IS (boundary_custom)
      CALL beast%SetCustomBoundaries(mesh,physics, &
        (/custom_nograd,custom_outflow,custom_kepler,custom_nograd/))
    END SELECT

    ! get gravitational acceleration
    sp => sim%Sources
    DO
      IF (ASSOCIATED(sp).EQV..false.) RETURN
      SELECT TYPE(sp)
      CLASS IS(sources_gravity)
        gp => sp
        EXIT
      END SELECT
      sp => sp%next
    END DO
    IF (.NOT.ASSOCIATED(sp)) CALL physics%Error("pringle::InitData","no gravity term initialized")

    SELECT CASE(vistype)
    CASE(beta)
      mu=1.0
    CASE(powerlaw)
      mu=pl_exp
    CASE DEFAULT
      CALL timedisc%Error("pringle::InitData","only pringle and beta viscosity possible")
    END SELECT

    ! initial condition
    SELECT TYPE(pvar => timedisc%pvar)
    TYPE IS(statevector_eulerisotherm)
      DO k=mesh%KGMIN,mesh%KGMAX
        DO j=mesh%JGMIN,mesh%JGMAX
          DO i=mesh%IMIN,mesh%IGMAX
            ! distance to axis (r-coordinate in cylindrical coordinates)
            r = bccyl(i,j,k,1)
            ! z-coordinate
            z = bccyl(i,j,k,3)
            ! Keplerian velocity
            vphi = r/(r**2+z**2)**0.75 
            ! compute pressure scale height
            height = csiso * sqrt(r**3) ! h = cs / Omega = cs / SQRT(GM/r**3), with GM=1
#ifdef HAVE_FGSL
            ! use exact solution at time TINIT
            CALL viscousring_analytic(mu,tinit,r,sden,vr)
#else
            ! use approximation for TINIT << 1
            ! no need for modified Bessel functions
            CALL viscousring_approx(mu,tinit,r,sden,vr)
#endif
            pvar%density%data3d(i,j,k) = rhomin + sden / (sqrt(2*pi) * height) &
              * exp(-0.5*(z/height)**2)
!             pvar%velocity%data4d(i,j,k,1) = vr
!             pvar%velocity%data4d(i,j,k,2) = vphi
!             pvar%velocity%data4d(i,j,k,3) = 0.0
          END DO
        END DO
      END DO
      ! determine disk mass
      mass = sum(mesh%volume%data1d(:)*pvar%density%data1d(:),mask=mesh%without_ghost_zones%mask1d)
#ifdef PARALLEL
      CALL mpi_allreduce(mpi_in_place,mass,1,default_mpi_real,mpi_sum, &
                         mesh%comm_cart,ierror)
#endif
      ! rescale disk mass
      pvar%density%data1d(:) = pvar%density%data1d(:) * mdisk / mass
      pvar%velocity%data1d(:) = 0.0 ! initialize all velocities with zero
      CALL gp%UpdateGravity(mesh,sim%Physics,sim%Fluxes,pvar,0.0,0.0)
      pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
        timedisc%GetCentrifugalVelocity(mesh,sim%Physics,sim%Fluxes,sim%Sources,(/0.,0.,1./), &
          gp%accel%data4d)
    CLASS DEFAULT
      CALL timedisc%Error("pringle::InitData","only isothermal physics possible")
    END SELECT

    IF (ASSOCIATED(timedisc%solution)) THEN
      ! store initial condition in exact solution array
      timedisc%solution = timedisc%pvar
    END IF

    ! transform to conservative variables
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)

    CALL mesh%Info(" DATA-----> initial condition: " // "pringle disk")
    WRITE(value,"(ES9.3)") tvis
    CALL mesh%Info("                               " // "viscous timescale:  " //trim(value))
    WRITE(value,"(ES9.3)") re
    CALL mesh%Info("                               " // "Reynolds number:    " //trim(value))
    WRITE(value,"(ES9.3)") ma
    CALL mesh%Info("                               " // "Mach number:        " //trim(value))
    WRITE(value,"(ES9.3)") mdisk
    CALL mesh%Info("                               " // "Disk mass:          " //trim(value))

    CALL geo_cyl%Finalize()
    DEALLOCATE(geo_cyl)
  END SUBROUTINE initdata


  SUBROUTINE viscousring_approx(mu,tau,r,Sigma,vr)
    IMPLICIT NONE
    !--------------------------------------------------------------------!
    REAL, INTENT(IN)  :: mu,tau,r
    REAL, INTENT(OUT) :: Sigma,vr
    !--------------------------------------------------------------------!
    REAL :: lambda,x,r4l,invr4l
    !--------------------------------------------------------------------!
    lambda = 1./(4.0-mu)
    r4l    = r**(0.25/lambda)
    x      = 2*lambda**2 / tau * r**(0.25/lambda)
    invr4l = 1./r4l
    sigma  = mdisk/sqrt((4*pi)**3 * tau) * r4l**(1.5-9*lambda) &
               * exp(-lambda**2/tau*(1.0-r4l)**2)
    vr     = -1.5/re * ((x*tau)/(2*lambda*lambda))**(4*lambda-2.) &
               * (1.0 - 0.5*x/lambda * (x*tau/(2*lambda*lambda) - 1.0))
  END SUBROUTINE viscousring_approx

#ifdef HAVE_FGSL
  SUBROUTINE viscousring_analytic(mu,tau,r,Sigma,vr)
    IMPLICIT NONE
    !--------------------------------------------------------------------!
    REAL, INTENT(IN)  :: mu,tau,r
    REAL, INTENT(OUT) :: Sigma,vr
    !--------------------------------------------------------------------!
    REAL :: lambda,r4l,x
    !--------------------------------------------------------------------!
    lambda = 1./(4.0-mu)
    r4l    = r**(0.25/lambda)
    x      = 2*lambda**2 / tau * r4l
    sigma  = mdisk/(4*pi) * (lambda/tau) * r4l**2 * r**(-2.25) &
               * exp(-lambda**2/tau*(1.0+r4l**2)) &
               * modified_bessel_inu(lambda,x)
    vr     = -1.5/re * ((x*tau)/(2*lambda*lambda))**(4*lambda-2.) &
               * (1.0 - 0.5*x/lambda * (x*tau/(2*lambda*lambda) &
               - modified_bessel_inu(1+lambda,x)/modified_bessel_inu(lambda,x)))
  END SUBROUTINE viscousring_analytic

  FUNCTION modified_bessel_inu(nu,x) RESULT(Inu)
    USE fgsl
    IMPLICIT NONE
    !--------------------------------------------------------------------!
    REAL, INTENT(IN) :: nu,x
    REAL             :: inu
    !--------------------------------------------------------------------!
    INTEGER          :: n
    !--------------------------------------------------------------------!
    IF (aint(nu).EQ.nu) THEN
      ! nu is an integer -> use modified bessel functions of integer order
      n = int(nu)
      inu = fgsl_sf_bessel_icn(n,x)
    ELSE
      ! nu is not an integer -> use representation of modified bessel functions
      ! through hypergeometric functions 0F1
      inu = (0.5*x)**nu / fgsl_sf_gamma(1.0+nu) * fgsl_sf_hyperg_0f1(1.0+nu,(0.5*x)**2)
    END IF
  END FUNCTION modified_bessel_inu
#endif

END PROGRAM pringle3d