mmsn.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: bindisk.f90                                                       #
!#                                                                           #
!# Copyright (C) 2010-2017                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Anna Feiler      <afeiler@astrophysik.uni-kiel.de>                        #
!# Roman Avramenko  <ravramenko@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 init
  USE fosite_mod
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER    :: gn      = 1.0             ! Newtons grav. constant    !
  REAL, PARAMETER    :: au      = 1.0             ! astronomical unit         !
  REAL, PARAMETER    :: msun    = 1.0             ! solar mass [kg]              !
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 10              ! simulation time [binary orbits]
  REAL, PARAMETER    :: valpha  = 1e-3            ! alpha viscosity parameter !
  ! 1. binary system
!  REAL, PARAMETER    :: MBH1    = 0.690*MSUN
!  REAL, PARAMETER    :: MBH2    = 0.203*MSUN
  REAL, PARAMETER    :: mbh1    = 0.4465*msun
  REAL, PARAMETER    :: mbh2    = 0.4465*msun
!  REAL, PARAMETER    :: EXCENT  = 0.159          ! excentricity              !
  REAL, PARAMETER    :: excent  = 0.0            ! excentricity              !
  REAL, PARAMETER    :: semma   = 0.224*au       ! semi mayor axis           !
  ! 2. disk
  REAL, PARAMETER    :: rgap    = 2.5*semma      ! estimated gap size        !
  REAL               :: xscale  = 1.0            ! scaling factor
  REAL               :: hratio  = 0.05           ! scaling factor
  ! gas paramter
  REAL, PARAMETER    :: mu      = 2.35e-3        ! mean molecular mass [kg/mol] !
  REAL, PARAMETER    :: rg      = 8.31447        ! molar gas constant        !
  ! mesh settings
  REAL, PARAMETER    :: gpar = au                ! geometry scaling paramete !
  REAL, PARAMETER    :: rmin = 1.5*semma         ! inner radius of the disk  !
  REAL, PARAMETER    :: rmax = 5.0*au            ! outer radius of the grid  !
  INTEGER, PARAMETER :: xres = 128               ! x-resolution              !
  INTEGER, PARAMETER :: yres = 128               ! y-resolution              !
  INTEGER, PARAMETER :: zres = 1                 ! y-resolution              !
  ! output file parameter
  INTEGER, PARAMETER :: onum = 100               ! number of output time st  !
  CHARACTER(LEN=256), PARAMETER &
                     :: odir  = "./"
  CHARACTER(LEN=256), PARAMETER &                ! output data file name     !
                     :: ofname = 'mmsn'
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  REAL               :: omega,period
  !--------------------------------------------------------------------------!

ALLOCATE(sim)
CALL sim%InitFosite()
CALL makeconfig(sim, sim%config)
CALL sim%Setup()
CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, sim%Fluxes, sim%Sources)
CALL sim%Run()
CALL sim%Finalize()
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, binary, vis, timedisc, datafile, &
                                    rotframe
    !------------------------------------------------------------------------!
    ! some derived simulation parameters
    omega  = sqrt(gn*(mbh1+mbh2)/semma)/semma
    period = 2.*pi / omega
!     OMEGA  = 0.0

    ! mesh settings
    ! stellar orbits must be inside the central hole of the mesh
    mesh => dict( &
              "meshtype"        / midpoint, &
              "geometry"        / logcylindrical, &
!              "geometry"        / CYLINDRICAL, &
              "inum"            / xres, &
              "jnum"            / yres, &
              "knum"            / zres, &
              "xmin"            / log(rmin/gpar), &
              "xmax"            / log(rmax/gpar), &
!              "xmin"            / (RMIN/GPAR), &
!              "xmax"            / (RMAX/GPAR), &
              "ymin"            / 0.0, &
              "ymax"            / (2.0*pi), &
              "zmin"            / 0.0, &
              "zmax"            / 0.0, &
              "use_fargo"       / 1, &
              "fargo"           / 1, &
              "decomposition"   / (/ -1, 1, 1/), &
              "gparam"          / gpar)
    IF (omega.GT.0.0) CALL setattr(mesh,"omega",omega)

    ! physics settings
    physics => dict( &
              "problem"         / euler_isotherm, &
              "units"           / geometrical,&
              "output/fcsound"  / 1, &
              "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
    vis => dict( &
              "stype"           / viscosity, &
              "vismodel"        / alpha, &
              "cvis"            / 0.1,&
              "dynconst"        / valpha, &
              "output/dynvis"   / 1)

    rotframe => dict( &
              "stype"           / rotating_frame)

    ! gravitational acceleration due to binary system
    binary => dict( &
              "gtype"           / pointmass_binary, &
              "mass1"           / mbh1, &
              "mass2"           / mbh2, &
              "mesh_rot"        / 1, &
!              "excentricity"    / EXCENT, &
              "output/binpos"   / 1, &
              "output/omega"    / 1, &
              "semimayoraxis"   / semma)

    ! source term due to all gravity terms
    grav => dict( &
              "stype"           / gravity, &
              "binary"          / binary,&
!              "self/gtype"      / SPECTRAL, &
!              "self/green"      / 1, &
              "output/height"   / 1, &
              "output/potential"/ 1, &
              "output/accel"    / 1)


    sources => dict( &
              "vis"             / vis, &
              "grav"            / grav)
!     IF (OMEGA.GT.0.0) CALL SetAttr(sources,"rotframe",rotframe)

    ! time discretization settings
    timedisc => dict( &
              "method"          / ssprk,&
              "cfl"             / 0.3, &
              "stoptime"        / (period*tsim), &
              "dtlimit"         / 1.0e-40,&
              "tol_rel"         / 1.0e-3, &
              "tol_abs"         / (/ 1.0e-08, 1., 1.0e-08 /), &
              "rhstype"         / 1, &
              "maxiter"         / 100000000)

    ! initialize data input/output
    datafile => dict( &
              "fileformat"      / vtk, &
              "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)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    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_base),  POINTER        :: Sources
    !------------------------------------------------------------------------!
    ! Local variable declaration
    CLASS(sources_base), POINTER :: sp
    CLASS(sources_gravity), POINTER :: gp
    INTEGER           :: i
#ifdef PARALLEL
    INTEGER           :: ierror
#endif
    REAL              :: Sigma0
    REAL, DIMENSION(:,:,:), POINTER :: r, Sigma
    REAL, DIMENSION(:,:,:,:), POINTER :: r_faces
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX)    :: bccsound
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,6)  :: fcsound
    CHARACTER(LEN=32) :: info_str
    !------------------------------------------------------------------------!
    ! set some pointers for convenience
    r => mesh%RemapBounds(mesh%radius%bcenter)
    r_faces => mesh%RemapBounds(mesh%radius%faces)

    ! Set sound speed
    bccsound = hratio*sqrt((mbh1+mbh2)*physics%constants%GN/r(:,:,:))
    DO i =1,6
      fcsound(:,:,:,i) = hratio*sqrt((mbh1+mbh2)*physics%constants%GN/r_faces(:,:,:,i))
    END DO

    ! set isothermal sound speeds
    SELECT TYPE (phys => physics)
    CLASS IS(physics_eulerisotherm)
      CALL phys%SetSoundSpeeds(mesh,bccsound)
      CALL phys%SetSoundSpeeds(mesh,fcsound)
    CLASS DEFAULT
      ! abort
      CALL phys%Error("InitData","physics not supported")
    END SELECT

    ! boundary conditions
    ! custom boundary conditions at western boundary 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/))
    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/))
    END SELECT

    ! 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("mmsn::InitData","no gravity term initialized")

    ! choose Sigma0 in a way that 2% of the mass of the binary is contained in 30 au
!    Sigma0     = 0.0008*(MBH1+MBH2)*SQRT(30.0/AU)/(4.*PI)
!    Sigma0     = 0.02*(MBH1+MBH2)/(4.*PI*(-(30.0/GPAR)**(-0.5)+(RMIN/GPAR)**(-0.5)))
!    Sigma0     = 0.02*(MBH1+MBH2)*SQRT(30.0/AU)/(4.*PI)
    sigma0 = 0.00311
    SELECT TYPE (pvar => timedisc%pvar)
    CLASS IS(statevector_eulerisotherm)
      ! 1. set surface density and initialize gravitational acceleration
      pvar%density%data3d(:,:,:) = fgap(r(:,:,:),rgap)*sigma0*0.1*xscale*r(:,:,:)**(-1.5)
!       pvar%velocity%data4d(:,:,:,:) = 0.0
!       IF (Mesh%FARGO.GT.0) Timedisc%w(:,:) = 0.0
!       ! get conservative variables
!       CALL Physics%Convert2Conservative(Timedisc%pvar,Timedisc%cvar)
      ! 2. azimuthal velocity: balance initial radial gravitational acceleration
      !    with centrifugal acceleration
      CALL gp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
      pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
        timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),gp%accel%data4d)
      ! transform velocities to rotating frame
      ! ATTENTION: this works only if the second velocity is the azimuthal velocity
      pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) - mesh%OMEGA*r(:,:,:)
    CLASS DEFAULT
      CALL physics%Error("mmsn::InitData","unsupported state vector")
    END SELECT

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

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

    ! print some information on stdout
    CALL mesh%Info(" DATA-----> initial condition: " // "MMSN - Setup")
    WRITE (info_str, '(ES8.2)') xscale
    CALL mesh%Info("                    disk mass: " // trim(info_str) // " MMSN")

  END SUBROUTINE initdata

  ELEMENTAL REAL FUNCTION fgap(R, Rgap)
    IMPLICIT NONE
    REAL, INTENT(IN) :: r, rgap

    fgap = 1./(1. + exp(-(r-rgap)/(0.1*rgap)))

  END FUNCTION

END PROGRAM init