pringle3d.f90

Jim Pringles famous viscous spreading ring solution in 3D.

Jim Pringles famous viscous spreading ring solution in 3D

Author

Björn Sperling

Tobias Illenseer

This variant simulates the ring with vertical extent with and without rotational symmetry. The equations are solved in non-dimensional units, using the Keplerian velocity at the location of the initial ring at R=1 as the velocity scale.

Remarks

(1) This test is very sensitive to mesh settings. If there are not enough mesh points at small radii instabilities grow starting at the inner boundary.

(2) There is an unsolved stability issue with the viscous source term. The stability condition doesn't seem to be sufficient. Thus in case of high Reynolds numbers one has to adjust the viscous CFL number cvis to enforce smaller time steps and preserve stability.

References:

• [32] Lynden-Bell, D. & Pringle, J. E. Evolution of Viscous Disks and Origin of Nebular Variables, M.N.R.A.S vol. 168(3) (1974) pp. 603-637 http://adsabs.harvard.edu/abs/1974MNRAS.168..603L
• [42] Pringle, J. E. Accretion discs in astrophysics Annual review of astronomy and astrophysics. vol. 19 (1981) pp. 137-162 DOI: 10.1146/annurev.aa.19.090181.001033
• [50] Speith, R. and Kley, W. Stability of the viscously spreading ring Astron. Astrophys. vol. 399 (2003), pp. 395-407 DOI: 10.1051/0004-6361:20021783

!#############################################################################
!#                                                                           #
!# 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 sources_base_mod, ONLY : sources_base
    USE sources_gravity_mod, ONLY : sources_gravity
    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 => null()
    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%GetSourcesPointer(gravity)
    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,"(ES11.3)") tvis
    CALL mesh%Info("                               " // "viscous timescale:  " //trim(value))
    WRITE(value,"(ES11.3)") re
    CALL mesh%Info("                               " // "Reynolds number:    " //trim(value))
    WRITE(value,"(ES11.3)") ma
    CALL mesh%Info("                               " // "Mach number:        " //trim(value))
    WRITE(value,"(ES11.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