mestel.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: mestel.f90                                                        #
!#                                                                           #
!# Copyright (C) 2012-2019                                                   #
!# Manuel Jung    <mjung@astrophysik.uni-kiel.de>                            #
!# Björn 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 mestel
  USE fosite_mod
#ifdef NECSXAURORA
    USE asl_unified
#endif
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER :: gn = 6.67384e-11       ! Newtons grav. constant [m^3/kg/s^2]
  REAL, PARAMETER :: year = 3.15576e+7      ! year [s]
  REAL, PARAMETER :: parsec = 3.0857e+16    ! parsec [m]
  REAL, PARAMETER :: au = 1.49598e+11       ! astronomical unit [m]
  REAL, PARAMETER :: msun = 1.989e+30       ! solar mass [kg]
  REAL, PARAMETER :: rg = 8.31              ! gas constant
  ! simulation parameters
  REAL, PARAMETER :: simtime = 1.0e+4*year  ! simulation time [s]
  REAL, PARAMETER :: mbh = 1.0e+7*msun      ! initial black hole mass [kg]
  REAL, PARAMETER :: mratio = 0.1           ! initial mdisk/mbh mass ratio
  REAL, PARAMETER :: mdisk = mratio*mbh     ! initial disk mass [kg]
  REAL, PARAMETER :: temp = 100.0           ! initial temperature [K]
  REAL, PARAMETER :: noise = 0.3            ! initial noise level
  ! physics settings
  INTEGER, PARAMETER :: phys = euler        ! transport model
  REAL, PARAMETER :: mu = 6.02e-4           ! mean molecular weight [kg/mol]
  REAL, PARAMETER :: gamma = 1.4            ! ratio of specific heats
  REAL, PARAMETER :: beta_vis = 1.0e-3      ! beta viscosity parameter
  ! mesh settings
  REAL, PARAMETER :: rmin = 5.0e-2 * parsec ! inner radius [m]
  REAL, PARAMETER :: rmax = 1.e+0 * parsec  ! outer radius [m]
  REAL, PARAMETER :: rgeo = 1.0 * parsec    ! geometry scaling constant
  INTEGER, PARAMETER :: mgeo = logcylindrical ! mesh geometry
  INTEGER, PARAMETER :: xres = 64         ! mesh resolution (radial)
  INTEGER, PARAMETER :: yres = 128         ! mesh resolution (azimuthal)
  INTEGER, PARAMETER :: zres = 1            ! mesh resolution (z)
  ! output settings
  INTEGER, PARAMETER :: onum = 100         ! number of output time steps
  CHARACTER(LEN=256), PARAMETER :: &
  ofname = 'markward', &                    ! data file name
  odir   = "./"                             ! output directory
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  !--------------------------------------------------------------------------!

ALLOCATE(sim)

#ifdef NECSXAURORA
CALL asl_library_initialize()
#endif

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

#ifdef NECSXAURORA
CALL asl_library_finalize()
#endif

DEALLOCATE(sim)

CONTAINS
  SUBROUTINE makeconfig(config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(Dict_TYP),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(Dict_TYP),POINTER :: mesh,boundary,timedisc,datafile,&
                              sources,fluxes,grav,physics,&
                              vis,cooling
    !------------------------------------------------------------------------!
    physics => dict( &
        "problem"         / phys, &
        "output/bccsound" / 1, &
        "mu"              / mu, &
        "gamma"           / gamma, &
        "units"           / si)

    fluxes => dict( &
        "fluxtype"        / kt, &
        "order"           / linear, &
        "variables"       / primitive, &
        "limiter"         / vanleer, &
        "theta"           / 1.2)

    mesh => dict( &
        "meshtype"        / midpoint, &
        "geometry"        / mgeo, &
        "inum"            / xres, &
        "jnum"            / yres, &
        "knum"            / zres, &
        "xmin"            / log(rmin/parsec), &
        "xmax"            / log(rmax/parsec), &
        "ymin"            / (-pi), &
        "ymax"            / ( pi), &
        "zmin"            / 0.0, &
        "zmax"            / 0.0, &
        "gparam"          / rgeo, &
        "use_fargo"       / 1, &
        "fargo"           / 1, &
        "decomposition"   / (/-1,1,1/), &
        "output/volume"   / 1 )

    boundary => dict( &
        "western"         / custom, &
        "eastern"         / custom, &
        "southern"        / periodic, &
        "northern"        / periodic, &
        "bottomer"        / reflecting, &
        "topper"          / reflecting)

    grav => dict( &
        "stype"           / gravity, &
        "cvis"            / 0.9, &
        "energy"          / 0, &
        "output/height"   / 1, &
        "self/gtype"      / spectral, &
        "self/green"      / 1, &
        "pmass/gtype"     / pointmass, &
        "pmass/potential" / newton, &
        "pmass/mass"      / mbh, &
        "pmass/outbound"  / 1)

    ! cooling model
    cooling => dict( &
        "stype"           / disk_cooling, &
!         "method"          / GRAY, &
        "method"          / gammie, &
        "b_cool"          / 10.0, &
        "Tmin"            / 3.0, &
        "rhomin"          / 1.0e-30, &
        "cvis"            / 0.1)

    ! viscosity model
    vis => dict( &
        "stype"           / viscosity, &
        "vismodel"        / beta, &
        "dynconst"        / beta_vis, &
        "output/stress"   / 1, &
        "output/dynvis"   / 1, &
        "output/kinvis"   / 1, &
        "cvis"            / 0.5)

    ! collect source terms
    sources => dict( &
        "diskcooling"     / cooling, &
!         "viscosity"       / vis, &
        "grav"            / grav)

    ! time discretization settings
    timedisc => dict( &
        "method"          / ssprk, &
        "tol_rel"         / 1.0e-3, &
        "cfl"             / 0.3, &
        "stoptime"        / simtime, &
        "dtlimit"         / 1.0e-4, &
        "tmin"            / 3.0, &
        "rhstype"         / 1, &
        "maxiter"         / 2000000000)

    ! add absolute error bounds and output fields depending on physics
    SELECT CASE(phys)
    CASE(euler)
       CALL setattr(timedisc, "tol_abs", (/1.0e-3, 1.0e-1, 1.0e-01, 1.0e+10/))
       CALL setattr(timedisc, "output/xmomentum", 0)
       CALL setattr(timedisc, "output/ymomentum", 0)
       CALL setattr(timedisc, "output/energy", 0)
       CALL setattr(timedisc, "output/rhs", 0)
       CALL setattr(timedisc, "output/external_sources", 0)
    CASE DEFAULT
       CALL error(sim,"MakeConfig","Physics model not supported.")
    END SELECT

    ! data file settings
    datafile => dict( &
        "fileformat"      / vtk, &
        "filename"        / "mestel", &
        "count"           / onum)

    ! create configuration
    config => dict( &
        "physics"         / physics, &
        "fluxes"          / fluxes, &
        "mesh"            / mesh, &
        "boundary"        / boundary, &
        "sources"         / sources, &
        "timedisc"        / timedisc, &
        "datafile"        / datafile)

  END SUBROUTINE makeconfig


  SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes,Sources)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(INOUT) :: Physics
    CLASS(fluxes_base),   INTENT(INOUT) :: Fluxes
    CLASS(sources_base),  POINTER       :: Sources
    !------------------------------------------------------------------------!
    ! Local variable declaration
    CLASS(sources_base), POINTER :: sp
    CLASS(sources_gravity), POINTER :: gp
    CLASS(marray_base), ALLOCATABLE :: rands,Sigma
#ifdef PARALLEL
    INTEGER :: ierror
#endif
    REAL    :: mass
    CHARACTER(LEN=20) :: mdisk_str
#ifdef NECSXAURORA
    INTEGER :: rng, n
#endif
    !------------------------------------------------------------------------!
    ALLOCATE(rands,sigma)
    ! get gravitational acceleration
    sp => 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("mestel::InitData","no gravity term initialized")

    ! get random numbers for density noise
    rands = marray_base()
#ifndef NECSXAURORA
    CALL initrandseed(physics)
    CALL random_number(rands%data1d)
#else
    CALL asl_random_create(rng, asl_randommethod_mt19937_64)
    CALL asl_random_distribute_uniform(rng)
    n = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)
    CALL asl_random_generate_d(rng, n, rands%data1d)
    CALL asl_random_destroy(rng)
#endif
    rands%data1d(:) = rands%data1d(:) * noise * 2.0 + (1.0 - noise)

    ! set surface density using radial power law (1/r) with a little noise
    sigma = marray_base()
    sigma%data1d(:) = rands%data1d(:) * rmin/mesh%radius%data2d(:,2)

    ! determine disk mass
    mass = sum(mesh%volume%data1d(:)*sigma%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

    ! setting for custom boundary conditions (western boundary)
    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

    ! setting for custom boundary conditions (eastern boundary)
    SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
    CLASS IS (boundary_custom)
      CALL beast%SetCustomBoundaries(mesh,physics, &
        (/custom_reflect,custom_reflect,custom_logexpol,custom_reflect/))
    END SELECT

    SELECT TYPE (pvar => timedisc%pvar)
    CLASS IS(statevector_euler)
      ! 1. set surface density by rescaling Sigma
      pvar%density%data1d(:) = sigma%data1d(:) * mdisk / mass
      ! 2. initialize gravitational acceleration
      CALL gp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
      ! 3. set azimuthal velocity: balance initial radial gravitational
      !    acceleration with centrifugal acceleration
      pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
        timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),gp%accel%data4d)
      ! 4. transform velocities to rotating frame
      ! ATTENTION: this works only if the second velocity is the azimuthal velocity
      IF (mesh%OMEGA.GT.0.0) &
        pvar%velocity%data2d(:,2) = pvar%velocity%data2d(:,2) - mesh%OMEGA*mesh%radius%data2d(:,2)
      ! 5. set pressure using surface density and initial temperature
      pvar%pressure%data1d(:) = physics%constants%RG/physics%MU * temp * pvar%density%data1d(:)
    CLASS DEFAULT
      CALL physics%Error("mestel::InitData","unsupported state vector")
    END SELECT

    IF (mesh%FARGO.EQ.2) &
       timedisc%w(:,:) = sqrt(physics%constants%GN*(mbh/mesh%radius%bcenter(:,mesh%JMIN,:)))

    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)

    ! print some information
    WRITE (mdisk_str, '(ES8.2)') mdisk/msun
    CALL mesh%Info(" DATA-----> initial condition: Mestel's disk")
    CALL mesh%Info("            disk mass:         " // trim(mdisk_str) // " M_sun")

    CALL rands%Destroy()
    CALL sigma%Destroy()
    DEALLOCATE(rands,sigma)
  END SUBROUTINE initdata

  SUBROUTINE initrandseed(Physics)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(physics_base), INTENT(IN) :: Physics
    INTEGER :: i, n, clock
    INTEGER, DIMENSION(:), ALLOCATABLE :: seed
    !------------------------------------------------------------------------!
    ! Initialize random number generator with a seed based on the systems time
    ! source: http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
    CALL random_seed(size = n)
    ALLOCATE(seed(n))
    CALL system_clock(count=clock)
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
#ifdef PARALLEL
    seed = seed + physics%GetRank()
#endif
    CALL random_seed(put = seed)
    DEALLOCATE(seed)
  END SUBROUTINE initrandseed
END PROGRAM mestel