sblintheo.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sblintheo.f90                                                     #
!#                                                                           #
!# Copyright (C) 2015-2018                                                   #
!# 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 sblintheo
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER    :: gn         = 1.0            ! grav. constant [GEOM]  !
  ! make a shearingsheet simulation and denote the direction of shearing
  ! (applies boundaries, ficitious forces and fargo automatically)
  INTEGER, PARAMETER :: shearsheet_direction = 1
  ! simulation parameter
  REAL, PARAMETER    :: omega      = 1.0            ! rotation at fid. point !
  REAL, PARAMETER    :: sigma0     = 1.0            ! mean surf.dens.        !
  REAL, PARAMETER    :: delsigma   = sigma0*5.0e-4  ! disturbance            !
  REAL, PARAMETER    :: tsim       = 4.54845485/omega ! simulation time (first maximum) !
  REAL, PARAMETER    :: gamma      = 2.0            ! dep. on vert. struct.  !
  REAL, PARAMETER    :: q          = 1.5            ! shearing parameter     !
  ! set (initial) speed of sound using the Toomre criterion and
  ! assuming marginal stabily; in non-isothermal simulations this is
  ! used to set the initial pressure (see InitData) whereas in isothermal
  ! simulations this is constant throughout the simulation
  REAL, PARAMETER    :: soundspeed = pi*gn*sigma0/omega
  ! mesh settings
  INTEGER, PARAMETER :: mgeo       = cartesian
  INTEGER, PARAMETER :: xres       = 128            ! amount of cells in x-  !
  INTEGER, PARAMETER :: yres       = 128            ! y-direction (rho/phi)  !
  INTEGER, PARAMETER :: zres       = 1
  REAL, PARAMETER :: thickness     = 0e-5           ! small value or zero for!
                                                    !   razor thin disks     !
  REAL               :: domainx    = 40.0           ! domain size [GEOM]     !
  REAL               :: domainy    = 40.0           ! domain size [GEOM]     !
  ! number of output time steps
  INTEGER, PARAMETER :: onum       = 10
  ! output directory and output name
  CHARACTER(LEN=256), PARAMETER :: odir   = "./"
  CHARACTER(LEN=256), PARAMETER :: ofname = "sblintheo"
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  REAL               :: maximum
  !--------------------------------------------------------------------------!
 
  tap_plan(1)
 
  ALLOCATE(sim)
 
  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()
 
  ! check amplitude
  IF (sim%Mesh%dz.GT.0.0) THEN
    maximum = maxval(sim%Timedisc%pvar%data4d(:,:,:,sim%Physics%DENSITY))*sim%Mesh%dz - sigma0
  ELSE
    maximum = maxval(sim%Timedisc%pvar%data4d(:,:,:,sim%Physics%DENSITY)) - sigma0
  END IF
 
  CALL sim%Finalize()
  DEALLOCATE(sim)
 
  tap_check_close(maximum, 0.02395931, 0.0003, "First max. < 3e-4 deviation")
  tap_done
 
CONTAINS
 
  SUBROUTINE makeconfig(Sim,config)
    IMPLICIT NONE
    !--------------------------------------------------------------------------!
    CLASS(fosite)           :: Sim
    TYPE(dict_typ), POINTER :: config
    !--------------------------------------------------------------------------!
    TYPE(dict_typ), POINTER :: mesh,physics,fluxes,&
                               grav,sources,timedisc,datafile
    REAL                    :: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
    !--------------------------------------------------------------------------!
    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.5*thickness*domainx
    zmax       = +0.5*thickness*domainx
 
    ! physics settings
    physics =>  dict(&
                "problem"     / euler_isotherm, &
                "cs"          / soundspeed, &
                "units"       / geometrical &
                )
 
    ! mesh settings
    mesh =>     dict(&
                "meshtype"    / midpoint, &
                "geometry"    / mgeo, &
                "shearingbox" / shearsheet_direction, &
                "inum"        / xres, &
                "jnum"        / yres, &
                "knum"        / zres, &
                "xmin"        / xmin, &
                "xmax"        / xmax, &
                "ymin"        / ymin, &
                "ymax"        / ymax, &
                "zmin"        / zmin, &
                "zmax"        / zmax, &
                "omega"       / omega, &
                "decomposition"/ (/ 1, -1, 1/), &
                "output/rotation" / 0, &
                "output/volume"   / 0, &
                "output/bh"   / 0, &
                "output/dl"   / 0,  &
                "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/accel"     / 1, &
                "self/output/potential" / 1 &
                )
 
    ! add gravitational sources
    sources =>  dict(&
                "grav"        / grav &
                )
 
    ! time discretization settings
    timedisc => dict( &
              "method"        / dormand_prince, &
              "cfl"           / 0.4, &
              "stoptime"      / tsim, &
              "dtlimit"       / 1.0e-40, &
              "tol_rel"       / 0.1, &
              "maxiter"       / 100000000)
 
    ! initialize data input/output
    datafile => dict( &
              "fileformat"    / vtk, &
              "filename"      / (trim(odir) // 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
    !------------------------------------------------------------------------!
    ! local variable declaration
    REAL              :: kx, ky
    !------------------------------------------------------------------------!
 
    !---------------------- linear theory test ------------------------------!
    kx = -2*(2*pi/domainx)
    ky = 2*pi/domainy
 
    ! initial condition
    SELECT TYPE(p => pvar)
    TYPE IS(statevector_eulerisotherm)
      IF(mesh%shear_dir.EQ.2)THEN
        p%density%data3d(:,:,:) = sigma0 &
          + delsigma*cos(kx*mesh%bcenter(:,:,:,1) + ky*mesh%bcenter(:,:,:,2))
        p%velocity%data2d(:,1) = 0.0
        p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1)
      ELSE IF(mesh%shear_dir.EQ.1)THEN
        p%density%data3d(:,:,:) = sigma0 &
          + delsigma*cos(kx*mesh%bcenter(:,:,:,2) - ky*mesh%bcenter(:,:,:,1))
        p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
        p%velocity%data2d(:,2) = 0.0
      END IF
      ! compute volume density from surface density
      ! for disks with vertical extent
      IF (mesh%dz.GT.0.0) p%density%data1d(:) = p%density%data1d(:) / mesh%dz
    TYPE IS(statevector_euler) ! non-isothermal HD
      IF(mesh%shear_dir.EQ.2)THEN
        p%density%data3d(:,:,:) = sigma0 &
          + delsigma*cos(kx*mesh%bcenter(:,:,:,1) + ky*mesh%bcenter(:,:,:,2))
        p%velocity%data2d(:,1) = 0.0
        p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1)
      ELSE IF(mesh%shear_dir.EQ.1)THEN
        p%density%data3d(:,:,:) = sigma0 &
          + delsigma*cos(kx*mesh%bcenter(:,:,:,2) - ky*mesh%bcenter(:,:,:,1))
        p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
        p%velocity%data2d(:,2) = 0.0
      END IF
      p%pressure%data3d(:,:,:) = soundspeed**2 * sigma0 / gamma
      ! compute volume density/pressure from surface density/pressure
      ! for disks with vertical extent
      IF (mesh%dz.GT.0.0) THEN
        p%density%data1d(:) = p%density%data1d(:) / mesh%dz
        p%pressure%data1d(:) = p%pressure%data1d(:) / mesh%dz
      END IF
    CLASS DEFAULT
      CALL physics%Error("sblintheo::InitData","only (non-)isothermal HD supported")
    END SELECT
 
    CALL physics%Convert2Conservative(pvar,cvar)
 
    !------------------------------------------------------------------------!
 
    CALL mesh%Info(" DATA-----> initial condition: " // &
         "Linear Theory Test - Shearingsheet")
  END SUBROUTINE initdata
END PROGRAM sblintheo