shearingsheet.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: shearingsheet.f90                                                 #
!#                                                                           #
!# Copyright (C) 2015-2022                                                   #
!# Jannes Klee      <jklee@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 shearingsheet
  USE fosite_mod
#ifdef NECSXAURORA
  USE asl_unified
#endif
#ifdef PARALLEL
#ifdef HAVE_MPI_MOD
  USE mpi
#endif
#endif
  IMPLICIT NONE
#ifdef PARALLEL
#ifdef HAVE_MPI_H
  include 'mpif.h'
#endif
#endif
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER    :: gn         = 1.0            ! grav. constant [GEOM]  !
  INTEGER, PARAMETER :: units      = geometrical    ! needs to be consistent !
  ! simulation parameter
  REAL, PARAMETER    :: omega      = 1.0            ! rotation at fid. point !
  REAL, PARAMETER    :: sigma0     = 1.0            ! mean surf.dens.        !
  REAL, PARAMETER    :: tsim       = 20./omega     ! simulation time        !
  REAL, PARAMETER    :: gamma      = 2.0            ! dep. on vert. struct.  !
  REAL, PARAMETER    :: beta_c     = 10.0           ! cooling parameter      !
!  REAL, PARAMETER    :: BETA_C     = 2.0           ! 2 -> collapse          !
  REAL, PARAMETER    :: q          = 1.5            ! shearing parameter     !
  ! mesh settings
  INTEGER, PARAMETER :: mgeo       = cartesian
  INTEGER, PARAMETER :: shear_direction = 1         ! enables shearingsheet with direcion !
  INTEGER, PARAMETER :: xres       = 256            ! cells in x-direction   !
  INTEGER, PARAMETER :: yres       = 256            ! cells in x-direction   !
  INTEGER, PARAMETER :: zres       = 1              ! cells in z-direction   !
  REAL               :: domainx    = 200.0          ! domain size [GEOM]     !
  REAL               :: domainy    = 200.0          ! domain size [GEOM]     !
  ! number of output time steps
  INTEGER, PARAMETER :: onum       = 20
  ! output directory and output name
  CHARACTER(LEN=256), PARAMETER :: odir   = "./"
  CHARACTER(LEN=256), PARAMETER :: ofname = "shearingsheet"
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  !--------------------------------------------------------------------------!
 
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%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)           :: Sim
    TYPE(dict_typ), POINTER :: config
    !--------------------------------------------------------------------------!
    ! local variable declaration
    TYPE(dict_typ), POINTER :: mesh,physics,fluxes,grav,cooling, &
                               sources,timedisc,datafile
    REAL :: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX, SOUNDSPEED
    !--------------------------------------------------------------------------!
    domainx    = domainx*gn*sigma0/(omega*omega)
    domainy    = domainy*gn*sigma0/(omega*omega)
    xmin       = -0.5*domainx
    xmax       = +0.5*domainx
    ymin       = -0.5*domainy
    ymax       = +0.5*domainy
    zmin       = 0.0
    zmax       = 0.0
    soundspeed = pi*gn*sigma0/omega ! Toomre-criterion
 
    ! physics settings
    physics =>  dict(&
                "problem"     / euler, &
                "gamma"       / gamma, &
                "units"       / geometrical &
                )
 
    ! mesh settings
    mesh =>     dict(&
                "meshtype"    / midpoint, &
                "geometry"    / mgeo, &
                "omega"       / omega, &
                "shearingbox" / shear_direction, &
                "decomposition" / (/ 1, -1, 1/), &
                "inum"        / xres, &
                "jnum"        / yres, &
                "knum"        / zres, &
                "xmin"        / xmin, &
                "xmax"        / xmax, &
                "ymin"        / ymin, &
                "ymax"        / ymax, &
                "zmin"        / zmin, &
                "zmax"        / zmax, &
                "Qshear"      / q &
                )
 
    ! fluxes settings
    fluxes =>   dict(&
                "order"       / linear, &
                "fluxtype"    / kt, &
                "variables"   / primitive, &
                "limiter"     / vanleer &
                )
 
    ! gravity settings (source term)
    grav =>     dict(&
                "stype"               / gravity, &
                "self/gtype"          / sboxspectral, &
                "self/output/phi"     / 0, &
                "self/output/accel_x" / 0, &
                "self/output/Fmass2D" / 0, &
                "self/output/accel_y" / 0 &
                )
 
    ! parametrized cooling from Gammie (2001)
    cooling =>  dict(&
                "stype"        / disk_cooling, &
                "method"       / gammie_sb, &
                "output/Qcool" / 1, &
                "b_cool"       / beta_c &
                )
 
    ! sources settings (contains source terms)
    sources =>  dict(&
                "cooling"     / cooling, &
                "grav"        / grav &
                )
 
    ! time discretization settings
    timedisc => dict(&
                "method"      / dormand_prince, &
                "cfl"         / 0.4, &
                "stoptime"    / tsim, &
                "dtlimit"     / 1e-40, &
                "output/external_sources" / 0, &
                "output/density_cvar" / 0, &
                "output/energy" / 0, &
                "output/xmomentum" / 0, &
                "output/ymomentum" / 0, &
                "output/rhs" / 0, &
                "output/geometrical_sources" / 0, &
                "maxiter"     / 100000000, &
                "tol_rel"     / 1.0e-3 &
                )
 
    ! data i/o settings
    datafile => dict(&
                "fileformat"  / vtk, &
                "filepath"    / trim(odir), &
                "filename"    / trim(ofname), &
                "count"       / onum &
                )
 
    ! overall config settings
    config =>   dict(&
                "mesh"        / mesh, &
                "physics"     / physics, &
                "fluxes"      / fluxes, &
                "sources"     / sources, &
                "timedisc"    / timedisc, &
                "datafile"    / datafile &
                )
  END SUBROUTINE makeconfig
 
  SUBROUTINE initdata(Mesh,Physics, pvar, cvar)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base),    INTENT(IN)  :: Mesh
    CLASS(physics_base), INTENT(IN)  :: Physics
    CLASS(marray_compound), POINTER, INTENT(INOUT) :: pvar,cvar
    !------------------------------------------------------------------------!
    INTEGER           :: i,j,k,err
    REAL              :: SOUNDSPEED
    REAL, ALLOCATABLE :: rands(:,:,:)
#ifdef NECSXAURORA
    INTEGER :: rng, n
#endif
    !------------------------ standard run ----------------------------------!
    SELECT TYPE(p => pvar)
    TYPE IS(statevector_euler)
      ! constant initial pressure determined by Q = 1
      soundspeed = pi*physics%Constants%GN*sigma0/omega ! Toomre-criterion
      ! constant initial density
      DO concurrent(i=1:SIZE(p%density%data1d))
        p%density%data1d(i) = sigma0
        p%pressure%data1d(i) = 2.5**2*pi*pi* &
            physics%Constants%GN**2.*sigma0**3./(gamma*omega**2)
      END DO
 
 
      ! create random numbers for setup of initial velocities
#ifndef NECSXAURORA
      CALL initrandseed(physics)
#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)
#endif
      ALLOCATE(rands(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),stat=err)
      IF (err.NE.0) CALL physics%Error("shearingsheet::InitData","Memory allocation failed")
 
      IF (mesh%shear_dir.EQ.2) THEN
#ifndef NECSXAURORA
        CALL random_number(rands)
#else
        CALL asl_random_generate_d(rng, n, rands)
#endif
        DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
          p%velocity%data4d(i,j,k,1) = (rands(i,j,k)-0.5)*0.1*soundspeed
        END DO
 
#ifndef NECSXAURORA
        CALL random_number(rands)
#else
        CALL asl_random_generate_d(rng, n, rands)
#endif
        DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
          p%velocity%data4d(i,j,k,2)  = -q*omega*mesh%bcenter(i,j,k,1) &
            + (rands(i,j,k)-0.5)*0.1*soundspeed
        END DO
      ELSE IF (mesh%shear_dir.EQ.1) THEN
#ifndef NECSXAURORA
        CALL random_number(rands)
#else
        CALL asl_random_generate_d(rng, n, rands)
#endif
        DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
          p%velocity%data4d(i,j,k,1)  = q*omega*mesh%bcenter(i,j,k,2) &
            + (rands(i,j,k)-0.5)*0.1*soundspeed
        END DO
 
#ifndef NECSXAURORA
        CALL random_number(rands)
#else
        CALL asl_random_generate_d(rng, n, rands)
#endif
        DO concurrent(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
          p%velocity%data4d(i,j,k,2) = (rands(i,j,k)-0.5)*0.1*soundspeed
        END DO
      END IF
      DEALLOCATE(rands)
    CLASS DEFAULT
      CALL physics%Error("shear::InitData","only non-isothermal HD supported")
    END SELECT
#ifdef NECSXAURORA
    CALL asl_random_destroy(rng)
#endif
    !------------------------------------------------------------------------!
    CALL physics%Convert2Conservative(pvar,cvar)
    CALL mesh%Info(" DATA-----> initial condition: " // &
         "Standard run shearingsheet")
  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 shearingsheet