pringle.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: pringle.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 pringle_test
  USE fosite_mod
  USE common_dict
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 2.0e-2     ! 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-20    ! 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 = 400        ! x-resolution
  INTEGER, PARAMETER :: yres = 1           ! y-resolution
  INTEGER, PARAMETER :: zres = 1           ! z-resolution
  REAL, PARAMETER    :: rmin = 0.05        ! 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 = 10          ! number of output data sets
  CHARACTER(LEN=256), PARAMETER &          ! output data dir
                     :: odir = './'
  CHARACTER(LEN=256), PARAMETER &          ! output data file name
                     :: ofname = 'pringle'
  ! derived parameters
  REAL, PARAMETER    :: csiso = 1./ma      ! isothermal speed of sound
  REAL               :: x0 = 0.0           ! cartesian position of point
  REAL               :: y0 = 0.0           !   mass (default: origin)
  REAL               :: z0 = 0.0           !
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE   :: sim
  REAL, DIMENSION(:), ALLOCATABLE :: sigma
  REAL :: sum_numer, sum_denom
  INTEGER :: n,den,vx,vy
  REAL    :: next_output_time
  INTEGER :: i,j,k
  LOGICAL :: ok
  !--------------------------------------------------------------------------!

  ALLOCATE(sim)
  CALL sim%InitFosite()

#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_plan(4)
#ifdef PARALLEL
  END IF
#endif

  CALL makeconfig(sim, sim%config)
  CALL sim%Setup()
  CALL initdata(sim,sim%Mesh, sim%Physics, sim%Timedisc)
  next_output_time = sim%Datafile%time
  IF (ASSOCIATED(sim%Timedisc%solution)) THEN
#ifdef HAVE_FGSL
    DO WHILE((sim%Timedisc%maxiter.LE.0).OR.(sim%iter.LE.sim%Timedisc%maxiter))
      ! advance numerical solution
      IF(sim%Step()) EXIT
      IF (next_output_time.LT.sim%Datafile%time) THEN
        ! compute exact solution for next output time step
        next_output_time = sim%Datafile%time
        SELECT TYPE(pvar => sim%Timedisc%solution)
        TYPE IS (statevector_eulerisotherm)
          k=sim%Mesh%KMIN
          DO j=sim%Mesh%JMIN,sim%Mesh%JMAX
            DO i=sim%Mesh%IMIN,sim%Mesh%IMAX
              IF (vistype.EQ.beta) THEN
                CALL viscousring_analytic(1.0,next_output_time/tvis+tinit, &
                  sim%Mesh%radius%bcenter(i,j,k), &
                  pvar%density%data3d(i,j,k),pvar%velocity%data4d(i,j,k,1))
              ELSE
                CALL viscousring_analytic(pl_exp,next_output_time/tvis+tinit, &
                  sim%Mesh%radius%bcenter(i,j,k), &
                  pvar%density%data3d(i,j,k),pvar%velocity%data4d(i,j,k,1))
              END IF
            END DO
          END DO
        END SELECT
      END IF
    END DO
#else
    CALL sim%Error("pringle","computation of analytical solution requires GSL support")
#endif
  ELSE
    CALL sim%Run()
  END IF
  ok = .NOT.sim%aborted
  ALLOCATE(sigma(sim%Physics%VNUM))
  ! compare with exact solution if requested
  IF (ASSOCIATED(sim%Timedisc%solution)) THEN
    DO n=1,sim%Physics%VNUM
      ! use L1 norm to estimate the deviation from the exact solution:
      !   Σ |pvar - pvar_exact| / Σ |pvar_exact|
      sum_numer = sum(abs(sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
                           sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n) &
                        -sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
                           sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
      sum_denom = sum(abs(sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
                           sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
#ifdef PARALLEL
      IF (sim%GetRank().GT.0) THEN
        CALL mpi_reduce(sum_numer,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
        CALL mpi_reduce(sum_denom,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
      ELSE
        CALL mpi_reduce(mpi_in_place,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
        CALL mpi_reduce(mpi_in_place,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
#endif
      sigma(n) = sum_numer / sum_denom
#ifdef PARALLEL
      END IF
#endif
    END DO
  ELSE
    sigma(:) = 0.0
  END IF

  den = sim%Physics%DENSITY
  vx  = sim%Physics%XVELOCITY
  vy  = sim%Physics%YVELOCITY

#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_check(ok,"stoptime reached")
! These lines are very long if expanded. So we can't indent it or it will be cropped.
tap_check_small(sigma(den),5.0e-02,"density deviation < 5%")
tap_check_small(sigma(vx),8.0e-02,"radial velocity deviation < 8%")
! skip azimuthal velocity deviation, because exact value is 0
tap_check_small(sigma(vy),2.0e-04,"azimuthal velocity deviation < 0.02%")
tap_done
#ifdef PARALLEL
  END IF
#endif

  IF (ALLOCATED(sigma)) DEALLOCATE(sigma)
  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
!       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) = no_gradients
      bc(top)    = no_gradients
      !!! 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.05*PI
      y2 = pi/2 ! + 0.05*PI
      z1 = -pi
      z2 = pi
!       bc(WEST)  = NO_GRADIENTS
!       bc(EAST)  = NO_GRADIENTS
      bc(west)  = custom
      bc(east)  = custom
      bc(south) = no_gradients
      bc(north) = no_gradients
      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/x"   / x0, &            ! cartesian position of point mass
            "pmass/y"   / y0, &
            "pmass/outbound" / 0)          ! disable accretion

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

    ! 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"   / (/0.0,1e-5,0.0/))

    ! 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)
    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
    INTEGER           :: i,j,k
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) &
                      :: posvec,ephi
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: radius
    REAL              :: r,vr,vphi,sden,mu
    CHARACTER(LEN=64) :: value
    !------------------------------------------------------------------------!
    IF (abs(x0).LE.tiny(x0).AND.abs(y0).LE.tiny(y0)) THEN
       ! no shift of point mass set radius and posvec to Mesh defaults
       radius(:,:,:) = mesh%radius%bcenter(:,:,:)
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
    ELSE
       ! shifted point mass position:
       ! compute curvilinear components of shift vector
       posvec(:,:,:,1) = x0
       posvec(:,:,:,2) = y0
       posvec(:,:,:,3) = z0
       CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,posvec,posvec)
       ! subtract the result from the position vector:
       ! this gives you the curvilinear components of all vectors pointing
       ! from the point mass to the bary center of any cell on the mesh
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - posvec(:,:,:,:)
       ! compute its absolute value
       radius(:,:,:) = sqrt(posvec(:,:,:,1)**2+posvec(:,:,:,2)**2+posvec(:,:,:,3)**2)
    END IF

    ! curvilinear components of azimuthal unit vector
    ! (maybe with respect to shifted origin)
    ! from ephi = ez x er = ez x posvec/radius = ez x (rxi*exi + reta*eeta)/r
    !             = rxi/r*(ez x exi) + reta/r*(ez x eeta) = rxi/r*eeta - reta/r*exi
    ! because (ez,exi,eeta) is right handed orthonormal set of basis vectors
    ephi(:,:,:,1) = -posvec(:,:,:,2)/radius(:,:,:)
    ephi(:,:,:,2) = posvec(:,:,:,1)/radius(:,:,:)
    ephi(:,:,:,3) = 0.0

    ! 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_extrapol,custom_kepler/))
    END SELECT
    SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
    CLASS IS (boundary_custom)
      CALL beast%SetCustomBoundaries(mesh,physics, &
        (/custom_logexpol,custom_outflow,custom_kepler/))
    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 center of mass
            r = radius(i,j,k)
            ! Keplerian velocity
            vphi = sqrt(1.0/r)
#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
            ! curvilinear velocity components
            pvar%velocity%data4d(i,j,k,1:2) = &
                vr*posvec(i,j,k,1:2)/r + vphi*ephi(i,j,k,1:2)
!               pvar%velocity%data4d(i,j,k,1:2) = vr*posvec(i,j,k,1:2)/r
!               p%velocity%data4d(i,j,k,2) = vphi
          END DO
        END DO
      END DO
!         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

    ! check fargo and set background velocity field
    IF (mesh%FARGO.EQ.2) &
       timedisc%w(:,:) = sqrt(1.0/radius(:,mesh%JMIN,:))

    ! 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))

  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 pringle_test