sblintheo3D.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sblintheo.f90                                                     #
!#                                                                           #
!# Copyright (C) 2015                                                        #
!# 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    :: TSIM       = 10.0/OMEGA     ! simulation time        !
  REAL, PARAMETER    :: gamma      = 1.4            ! 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 ! det. by. Toomre    !
  REAL, PARAMETER    :: height     = soundspeed/omega ! (initial) scale height
!   INTEGER, PARAMETER :: PHYS       = EULER
  INTEGER, PARAMETER :: phys       = euler_isotherm
  ! mesh settings
  INTEGER, PARAMETER :: mgeo       = cartesian
  INTEGER, PARAMETER :: xres       = 64             ! amount of cells in x-  !
  INTEGER, PARAMETER :: yres       = 64             ! y-direction (rho/phi)  !
  INTEGER, PARAMETER :: zres       = 16             ! and z-direction        !
  ! extent of computational domain
  REAL               :: domainx    = 40.0           ! domain size [GEOM]     !
  REAL               :: domainy    = 40.0           ! domain size [GEOM]     !
  REAL               :: domainz    = 10.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 = "sblintheo3D"
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  LOGICAL :: ok
  !--------------------------------------------------------------------------!
 
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()
  ok = .NOT.sim%aborted
 
! search for the amplitude
!maximum = MAXVAL(Sim%Timedisc%pvar(:,:,:,Sim%Physics%DENSITY)) - SIGMA0
!TAP_CHECK_CLOSE(maximum, 0.0316463969505, 0.005, "Last max. < 0.005 deviation")
 
  CALL sim%Finalize()
  DEALLOCATE(sim)
 
tap_check(ok,"stoptime reached")
tap_done
 
CONTAINS
 
  SUBROUTINE makeconfig(Sim,config)
    IMPLICIT NONE
    !--------------------------------------------------------------------------!
    CLASS(fosite)           :: Sim
    TYPE(dict_typ), POINTER :: config
    !--------------------------------------------------------------------------!
    TYPE(dict_typ), POINTER :: mesh,physics,fluxes,boundary,&
                               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*domainz
    zmax       = +0.5*domainz
 
    ! physics settings
    SELECT CASE(phys)
    CASE(euler_isotherm)
      physics => dict("problem" / euler_isotherm, &
                      "cs"      / soundspeed, &      ! isothermal speed of sound
                      "units"   / geometrical)
    CASE(euler)
      physics => dict("problem" / euler, &
                      "gamma"   / gamma, &           ! ratio of specific heats
                      "units"   / geometrical)
    CASE DEFAULT
      CALL sim%Error("sblintheo3D::MakeConfig","unsupported physics")
    END SELECT
 
    ! 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,  &
                "Qsher"       / q &
                )
 
    ! boundary conditions (x-/y-boundaries are set automatically)
    ! z-boundary should be reflecting, otherwise the disk mass
    ! increases and the disk becomes unstable
    boundary => dict(&
               "bottomer" / reflecting,&
               "topper"   / reflecting)
 
    ! 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"   / 0, &
                "self/output/potential" / 0 &
                )
 
    ! sources settings (contains source terms)
    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, &
                  "boundary"    / boundary, &
                  "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, kz
    !------------------------------------------------------------------------!
 
    !---------------------- linear theory test ------------------------------!
    kx = -2*(2*pi/domainx)
    ky = 2*pi/domainy
    kz = 2*pi/domainz
 
    ! initial condition
    SELECT TYPE(p => pvar)
    TYPE IS(statevector_eulerisotherm)
      ! zero all velocity components
      p%velocity%data1d(:) = 0.0
      IF(mesh%shear_dir.EQ.2)THEN
        p%density%data3d(:,:,:) = (sigma0 &
          + delsigma*cos(kx*mesh%bcenter(:,:,:,1) &
                       + ky*mesh%bcenter(:,:,:,2)))/(sqrt(2*pi)*height) &
                    *exp(-0.5*(mesh%bcenter(:,:,:,3)/height)**2)
        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)))/(sqrt(2*pi)*height) &
                    *exp(-0.5*(mesh%bcenter(:,:,:,3)/height)**2)
        p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
      END IF
      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)))/(sqrt(2*pi)*height) &
                    *exp(-0.5*(mesh%bcenter(:,:,:,3)/height)**2)
        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)))/(sqrt(2*pi)*height) &
                    *exp(-0.5*(mesh%bcenter(:,:,:,3)/height)**2)
        p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2)
      END IF
      p%density%data1d(:) = p%density%data1d(:) / mesh%dz
      p%pressure%data1d(:) = soundspeed**2 / gamma * p%density%data1d(:)
    CLASS DEFAULT
      CALL physics%Error("sblintheo3D::InitData","only (non-)isothermal HD supported")
    END SELECT
 
    CALL physics%Convert2Conservative(pvar,cvar)
    CALL mesh%Info(" DATA-----> initial condition: " // &
         "Linear Theory Test - 3D Shearingbox")
  END SUBROUTINE initdata
END PROGRAM sblintheo