sgdisk.f90

Program and data initialization a self-gravitating disk (sgdisk)

Program and data initialization a self-gravitating disk (sgdisk)

Author

Jannes Klee

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sgdisk.f90                                                        #
!#                                                                           #
!# Copyright (C) 2010-2024                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Jannes Klee      <jklee@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 sgdisk
  USE fosite_mod
#ifdef NECSXAURORA
    USE asl_unified
#endif
#ifdef HAVE_KROME
  USE krome_main
  USE krome_user
#endif
#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-8      ! grav. constant (cgs)     !
  REAL,     PARAMETER :: RG      = 8.31447d+7     ! molar gas constant (cgs) !
  REAL,     PARAMETER :: MSUN    = 1.989d+33      ! solar mass [g]           !
  REAL,     PARAMETER :: AU      = 1.49597870691e+13 ! ast. unit [cm]        !
  REAL,     PARAMETER :: YEAR    = 3.15576e+7     ! Julian year [s]          !
  REAL,     PARAMETER :: MSCALE  = msun           ! mass scaling param.      !
  ! simulation parameters
  REAL,     PARAMETER :: TSIM    = 6.0            ! simulation time [ORP]    !
  REAL,     PARAMETER :: MU      = 2.35d+0        ! mean mol. mass [g/mol]  !
  REAL,     PARAMETER :: VALPHA  = 0.05           ! alpha visc. parameter    !
  ! disk & central object
  REAL,     PARAMETER :: MBH     = 1.0*mscale     ! central mass             !
  REAL,     PARAMETER :: MDISK   = 0.1*mscale     ! initial disk mass        !
  REAL,     PARAMETER :: T_INIT  = 120.0          ! initial temperatur       !
  REAL,     PARAMETER :: T_MIN   = 10.0            ! minimal temperatur       !
  ! cooling parameters
  REAL,     PARAMETER :: BETA_C  = 10.0           ! cooling parameter        !
  REAL,     PARAMETER :: GAMMA   = 1.6            ! adiabatic index          !
  ! mesh settings
  INTEGER,  PARAMETER :: MGEO    = logcylindrical ! geometry of the mesh     !
  REAL,     PARAMETER :: GPAR    = au             ! geom. scal. parameter    !
  REAL,     PARAMETER :: RMIN    = 1.0*gpar       ! inner radius of the disk !
  REAL,     PARAMETER :: RMAX    = 25.0*gpar      ! outer radius of the grid !
  INTEGER,  PARAMETER :: XRES    = 256            ! x-resolution             !
  INTEGER,  PARAMETER :: YRES    = xres*3          ! y-resolution             !
  INTEGER,  PARAMETER :: ZRES    = 1              ! z-resolution             !
  ! mestel
  LOGICAL             :: DAMP                     ! wave damping             !
  REAL, PARAMETER     :: NOISE   = 0.3            ! initial noise level      !
  ! output file parameter
  INTEGER, PARAMETER  :: ONUM    = 100            ! number of outputs        !
  CHARACTER(LEN=256), PARAMETER &
                      :: ODIR    = "./"           ! output directory         !
  CHARACTER(LEN=256)  :: OFNAME  = 'sgdisk'       ! output name              !
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: Sim
  REAL                :: CSISO,OMEGA,ORP
  !--------------------------------------------------------------------------!
 
ALLOCATE(sim)
 
#ifdef NECSXAURORA
CALL asl_library_initialize()
#endif
 
CALL sim%InitFosite()
CALL makeconfig(sim, sim%config)
CALL sim%Setup()
CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, sim%Fluxes, sim%Sources, &
              sim%Timedisc%pvar, sim%Timedisc%cvar)
CALL sim%Run()
CALL sim%Finalize()
 
#ifdef NECSXAURORA
    CALL asl_library_finalize()
#endif
 
DEALLOCATE(sim)
 
CONTAINS
 
   SUBROUTINE makeconfig(Sim, config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite), INTENT(INOUT) :: Sim
    TYPE(Dict_TYP), POINTER      :: config
    TYPE(Dict_TYP), POINTER      :: mesh, physics, fluxes, boundary, grav,&
                                    sources, timedisc, datafile, cooling, &
                                    rotframe, pmass, self, damping
    !------------------------------------------------------------------------!
    ! some derived simulation parameters
    csiso   = sqrt(rg/mu*t_init)
    omega   = 0.0
    orp     = 2.*pi/(sqrt(gn*mbh/rmax**3.))
    damp    = .false.
 
    ! mesh settings
    ! stellar orbits must be inside the central hole of the mesh
    mesh => dict(&
                "meshtype"        / midpoint, &
                "geometry"        / mgeo, &
                "inum"            / xres, &
                "jnum"            / yres, &
                "knum"            / zres, &
                "xmin"            / log(rmin/gpar), &
                "xmax"            / log(rmax/gpar), &
                "ymin"            / 0.0, &
                "ymax"            / (2.0*pi), &
                "zmin"            / 0.0, &
                "zmax"            / 0.0, &
                "use_fargo"       / 1, &
                "fargo"           / 2, &
                "decomposition"   / (/ -1, 1, 1/), &
!                "output/radius"   / 1, &
!                "output/dl"       / 1, &
                "output/position_vector" / 1, &
                "gparam"          / gpar)
 
 
    ! physics settings
    physics => dict(&
                "problem"         / euler, &
                "mu"              / mu, &
                "gamma"           / gamma, &
                "units"           / cgs)
!                "output/bccsound" / 1)
 
    ! boundary conditions
    boundary => dict(&
                "western"         / custom,&
                "eastern"         / custom,&
                "southern"        / periodic, &
                "northern"        / periodic, &
                "bottomer"        / reflecting, &
                "topper"          / reflecting)
 
 
    ! numerical fluxes and reconstruction method
    fluxes => dict(&
                "order"           / linear, &
                "fluxtype"        / kt, &
                "variables"       / primitive, &
                "limiter"         / vanleer, &
                "theta"           / 1.2)
 
 
    ! viscosity source term
    rotframe => dict(&
                "stype"           / rotating_frame)
 
 
    ! gravitational acceleration due to binary system
    pmass => dict(&
                "gtype"           / pointmass, &
                "mass"            / mbh, &
                "output/position" / 1)
 
    ! cooling function
    cooling =>  dict(&
                "stype"           / disk_cooling, &
                "output/Qcool"    / 1, &
                "rhomin"          / 1.0e-30, &
                "method"          / gammie, &
                "b_cool"          / beta_c)
 
    ! self-gravity
    self => dict(&
                "gtype"           / spectral, &
                "output/potential" / 1, &
                "self/green"      / 1)
 
    ! collect all gravity-source terms
    grav => dict(&
                "stype"             / gravity, &
                "pmass"             / pmass, &
                "self"              / self, &
                "output/potential"  / 1, &
                "energy"            / 0, &
                "output/height"     / 1, &
                "output/accel"      / 1)
 
    ! collect all source terms
    sources => dict(&
                "cooling"           / cooling, &
                "grav"              / grav)
 
    ! add wave damping if requested
    IF (damp) &
       CALL setattr(sources,"damping",damping)
 
    ! time discretization settings
    timedisc => dict(&
                "method"            / modified_euler, &
                "cfl"               / 0.3, &
                "dtlimit"           / 1.0e-50, &
                "tol_rel"           / 1.0e-3, &
                "rhstype"           / 1, &
                "tol_abs"           / (/ 1.0e-16, 1.0, 1.0e-16, 1.0 /), &
                "output/energy"     / 1, &
                "output/rhs"        / 1, &
                "stoptime"          / (orp*tsim), &
                "pmin"              /  1e-14, &
                "tmin"              / t_min, &
                "checkdata"         / ior(check_nothing,check_tmin), &
                "maxiter"           / 100000000)
 
    ! initialize data input/output
    datafile => dict(&
               "fileformat"         / xdmf, &
               "filename"           / (trim(odir) // trim(ofname)), &
               "count"              / onum)
 
    config => dict(&
                "physics"           / physics, &
                "fluxes"            / fluxes, &
                "mesh"              / mesh, &
                "boundary"          / boundary, &
                "sources"           / sources, &
                "timedisc"          / timedisc, &
                "datafile"          / datafile)
 
  END SUBROUTINE makeconfig
 
 
  SUBROUTINE initdata(Mesh,Physics,Timedisc,Fluxes,Sources,pvar,cvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    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(timedisc_base), INTENT(INOUT) :: Timedisc
    CLASS(fluxes_base),   INTENT(INOUT) :: Fluxes
    CLASS(sources_list),ALLOCATABLE,INTENT(INOUT) :: Sources
    CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
    !------------------------------------------------------------------------!
    ! Local variable declaration
    CLASS(sources_base), POINTER :: sp => null()
    INTEGER           :: i,j,k
#ifdef PARALLEL
    INTEGER           :: ierror
#endif
    REAL              :: mass
    REAL, DIMENSION(:,:,:), POINTER :: r, Sigma
    REAL, DIMENSION(:,:,:,:), POINTER :: r_vec
    CHARACTER(LEN=32) :: mdisk_str
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
                      :: rands
#ifdef NECSXAURORA
    INTEGER :: rng, n
#endif
    !------------------------------------------------------------------------!
    ! MESTEL'S DISK
    ! distance from origin to cell bary centers and position vector
    r => mesh%RemapBounds(mesh%radius%bcenter(:,:,:))
    r_vec => mesh%RemapBounds(mesh%posvec%bcenter(:,:,:,:))
    ! pointer to density array
    sigma => mesh%RemapBounds(pvar%data4d(:,:,:,physics%DENSITY))
 
    ! set surface density using radial power law (1/r) with a little noise
#ifndef NECSXAURORA
    CALL initrandseed(timedisc)
    CALL random_number(rands)
#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)
#endif
    rands = rands * noise * 2.0 + (1.0 - noise)
    sigma = rands*(rmin/r(:,:,:))
 
#ifdef NECSXAURORA
    CALL asl_random_destroy(rng)
#endif
 
    ! determine disk mass
    mass = sum(mesh%volume%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) * &
               sigma(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
#ifdef PARALLEL
       CALL mpi_allreduce(mpi_in_place,mass,1,default_mpi_real,mpi_sum, &
            mesh%comm_cart,ierror)
#endif
    ! rescale disk mass
    sigma(:,:,:) = sigma(:,:,:) * mdisk / mass
 
    ! 2. azimuthal velocity: balance initial radial acceleration with centrifugal acceleration
    SELECT TYPE(phys => physics)
    TYPE IS(physics_euler)
      pvar%data4d(:,:,:,physics%PRESSURE) = &
      csiso*csiso*pvar%data4d(:,:,:,physics%DENSITY)/gamma
      ! necessary in order to calculte height below in UpdateGravity
      SELECT TYPE(pvarious => pvar)
      CLASS IS(statevector_euler)
        CALL phys%UpdateSoundSpeed(pvarious)
      END SELECT
    END SELECT
 
    ! reset velocities
    pvar%data2d(:,physics%XVELOCITY:physics%YVELOCITY) = 0.0
 
    ! check and update gravity
    IF (ALLOCATED(sources)) &
      sp => sources%GetSourcesPointer(gravity)
    IF (ASSOCIATED(sp)) THEN
      SELECT TYPE(sp)
      CLASS IS(sources_gravity)
        CALL sp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
        ! ATTENTION: Don't use GetCentrifugalVelocity without the optional acceleration array!
        ! This would yield undefined data, because GetCentrifugalVelocity calls ComputeRHS
        ! which calls CenterBoundary. Since the FARFIELD boundary conditions are not
        ! initialized at this stage (see below), the result is undefined.
        pvar%data4d(:,:,:,physics%XVELOCITY:physics%XVELOCITY+physics%VDIM-1) = &
            timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),sp%accel%data4d)
      CLASS DEFAULT
        CALL physics%Error("sgdisk::InitData","this should not happen -> gravity pointer is not of class sources_gravity")
      END SELECT
    ELSE
      CALL physics%Error("sgdisk::InitData","gravity source term not initialized")
    END IF
 
    IF (mesh%fargo%GetType().EQ.2) &
       timedisc%w(:,:) = sqrt(physics%constants%GN*mbh/r(:,mesh%JMIN,:))-mesh%OMEGA*r(:,mesh%JMIN,:)
 
    ! transform velocities to rotating frame
    pvar%data4d(:,:,:,physics%YVELOCITY) = pvar%data4d(:,:,:,physics%YVELOCITY) &
        - mesh%OMEGA*r(:,:,:)
 
    ! get conservative variables
    CALL physics%Convert2Conservative(pvar,cvar)
 
    ! 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_nograd,custom_outflow,custom_kepler,custom_nograd/))
    END SELECT
 
    CALL physics%Convert2Conservative(pvar,cvar)
    ! print some information
    WRITE (mdisk_str, '(ES10.2)') mdisk
    CALL mesh%Info(" DATA-----> initial condition: Mestel's disk")
    CALL mesh%Info("            disk mass:         " // trim(mdisk_str))
 
  END SUBROUTINE initdata
 
 
  SUBROUTINE initrandseed(Timedisc)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_base),INTENT(IN) :: Timedisc
    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 + timedisc%GetRank()
#endif
    CALL random_seed(put = seed)
    DEALLOCATE(seed)
  END SUBROUTINE initrandseed
 
 
END PROGRAM sgdisk