shear.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: shear.f90                                                         #
!#                                                                           #
!# Copyright (C) 2008-2018                                                   #
!# 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 shear
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER    :: gn         = 1.0            ! grav. constant [GEOM]  !
  INTEGER, PARAMETER :: units      = geometrical
  ! make a shearingsheet simulation and denote the direction of shearing
  ! (applies boundaries, ficitious forces and fargo automatically)
  INTEGER, PARAMETER :: shearsheet_direction = 2
  ! simulation parameter
  REAL, PARAMETER    :: csiso    = &
                                     0.0      ! non-isothermal simulation
!                                     1.0      ! isothermal simulation
                                               !   with CSISO as sound speed
  REAL, PARAMETER    :: omega      = 1.0
  REAL, PARAMETER    :: sigma0     = 1.0
  REAL, PARAMETER    :: tsim       = 10./omega
  REAL, PARAMETER    :: gamma      = 2.0
  REAL, PARAMETER    :: q          = 1.5            ! shearing parameter     !
  ! mesh settings
  INTEGER, PARAMETER :: mgeo       = cartesian
  INTEGER, PARAMETER :: xres       = 128
  INTEGER, PARAMETER :: yres       = 128
  INTEGER, PARAMETER :: zres       = 1
  REAL               :: domainx    = 320.0
  REAL               :: domainy    = 320.0
  ! number of output time steps
  INTEGER, PARAMETER :: onum       = 10
  ! output directory and output name
  CHARACTER(LEN=256), PARAMETER :: odir   = "./"
  CHARACTER(LEN=256), PARAMETER :: ofname = "sbox"
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  CLASS(marray_compound), POINTER :: pvar,pvar_init
  REAL               :: sigma
  !--------------------------------------------------------------------------!
  tap_plan(1)

  ! with west-east shear
  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)
  ! store transposed initial data
  CALL sim%Physics%new_statevector(pvar_init,primitive)
  CALL rotatedata(sim%Mesh,sim%Timedisc%pvar,pvar_init,"xy")
  CALL sim%Run()
  ! store transposed result of the first run
  CALL sim%Physics%new_statevector(pvar,primitive)
  CALL rotatedata(sim%Mesh,sim%Timedisc%pvar,pvar,"xy")
  ! finish the simulation
  CALL sim%Finalize()
  DEALLOCATE(sim)

  ! simulate south-north shear
  ALLOCATE(sim)
  CALL sim%InitFosite()
  CALL makeconfig(sim, sim%config)
  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_rotate"))
  CALL setattr(sim%config, "mesh/shear_dir", 1)
  CALL sim%Setup()
  sim%Timedisc%pvar%data1d(:) = pvar_init%data1d(:)
  CALL sim%Physics%Convert2Conservative(sim%Timedisc%pvar,sim%Timedisc%cvar)
  CALL sim%Run()

  ! compare results
  sigma = sqrt(sum((sim%Timedisc%pvar%data4d(:,:,:,:)-pvar%data4d(:,:,:,:))**2)/ &
                    SIZE(pvar%data4d(:,:,:,:)))

  CALL pvar%Destroy()
  CALL sim%Finalize()
  DEALLOCATE(pvar,sim)

  tap_check_small(sigma,1e-13,"Shear x-y symmetry test")
  tap_done

  CONTAINS

  SUBROUTINE makeconfig(Sim, config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite)            :: Sim
    TYPE(Dict_TYP), POINTER  :: config
    !------------------------------------------------------------------------!
    TYPE(Dict_TYP), POINTER  :: mesh, physics, datafile, &
                                sources, timedisc, fluxes, shear
    !------------------------------------------------------------------------!
    INTENT(INOUT)            :: sim
    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.0
    zmax       = 0.0

    ! mesh settings
    mesh =>     dict(&
                "meshtype"    / midpoint, &
                "geometry"    / mgeo, &
                "shear_dir"   / shearsheet_direction, &
                "inum"        / xres, &
                "jnum"        / yres, &
                "knum"        / zres, &
                "xmin"        / xmin, &
                "xmax"        / xmax, &
                "ymin"        / ymin, &
                "ymax"        / ymax, &
                "zmin"        / zmin, &
                "zmax"        / zmax, &
                "omega"       / omega, &
                "output/rotation" / 0, &
                "output/volume"   / 0, &
                "output/bh"   / 0, &
                "output/dl"   / 0  &
                )

    ! physics settings
    IF (csiso.GT.tiny(csiso)) THEN
      physics => dict("problem" / euler_isotherm, &
                        "units" / units, &
                           "cs" / csiso)             ! isothermal speed of sound
    ELSE
      physics => dict("problem" / euler, &
                        "units" / units, &
                        "gamma" / gamma)             ! ratio of specific heats
    END IF

    ! flux calculation and reconstruction method
    fluxes => dict( &
              "fluxtype"      / kt, &
              "order"         / linear, &
              "variables"     / primitive, &
              "limiter"       / vanleer &
              )

    ! fluxes settings
    fluxes =>   dict(&
                "order"       / linear, &
                "fluxtype"    / kt, &
                "variables"   / primitive, &
                "limiter"     / vanleer &
                )

    ! shearing box fictious forces
    shear => dict( &
                "stype"           / shearbox, &
                "output/accel_x"  / 0, &
                "output/accel_y"  / 0 &
                )

    ! collect sources in dictionary
    sources => dict("shear"   / shear)

    ! time discretization settings
    timedisc => dict( &
              "method"        / modified_euler, &
              "order"         / 3, &
              "cfl"           / 0.4, &
              "stoptime"      / tsim, &
              "dtlimit"       / 1.0e-40, &
              "tol_rel"       / 0.1, &
              "output/external_sources" / 1, &
              "maxiter"       / 100000000)

    ! initialize data input/output
    datafile => dict( &
              "fileformat"    / vtk, &
              "filename"      / (trim(odir) // trim(ofname)), &
              "count"         / onum)

    ! collect all above dicts in the configuration dict
    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
    !------------------------------------------------------------------------!
    ! Test for shearing and boundary module
    ! initial condition
    SELECT TYPE(p => pvar)
    TYPE IS(statevector_eulerisotherm) ! isothermal HD
      WHERE ((abs(mesh%bcenter(:,:,:,1)).LT.0.15*domainx.AND. &
              abs(mesh%bcenter(:,:,:,2)).LT.0.15*domainy))
        p%density%data3d(:,:,:) = sigma0
      ELSEWHERE
        p%density%data3d(:,:,:) = sigma0*1e-2
      END WHERE
      p%velocity%data2d(:,1) = 0.0
      p%velocity%data4d(:,:,:,2) = -mesh%Q*mesh%bcenter(:,:,:,1)*mesh%Omega
    TYPE IS(statevector_euler) ! non-isothermal HD
      WHERE ((abs(mesh%bcenter(:,:,:,1)).LT.0.15*domainx.AND. &
              abs(mesh%bcenter(:,:,:,2)).LT.0.15*domainy))
        p%density%data3d(:,:,:) = sigma0
      ELSEWHERE
        p%density%data3d(:,:,:) = sigma0*1e-2
      END WHERE
      p%velocity%data2d(:,1) = 0.0
      p%velocity%data4d(:,:,:,2) = -mesh%Q*mesh%bcenter(:,:,:,1)*mesh%Omega
      p%pressure%data1d(:) = 1e-1
    CLASS DEFAULT
      CALL physics%Error("shear::InitData","only non-isothermal HD supported")
    END SELECT

    CALL physics%Convert2Conservative(pvar,cvar)
    CALL mesh%Info(" DATA-----> initial condition: " // &
         "Shearing patch")
  END SUBROUTINE initdata

  ! only for symmetric matrix
  SUBROUTINE rotatedata(Mesh,pvar_in,pvar_out,dir)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base), INTENT(IN)            :: Mesh
    CLASS(marray_compound), INTENT(INOUT)   :: pvar_in,pvar_out
    CHARACTER(LEN=2), INTENT(IN)            :: dir
    !------------------------------------------------------------------------!
    INTEGER :: i,j,k
    !------------------------------------------------------------------------!
    SELECT TYPE(pin => pvar_in)
    TYPE IS(statevector_euler) ! non-isothermal HD
      SELECT TYPE(pout => pvar_out)
      TYPE IS(statevector_euler) ! non-isothermal HD
        SELECT CASE(dir)
        CASE("xy")
          ! rotate at middle of the field, because x' = -y in shearingsheet.
          DO k = mesh%KGMIN,mesh%KGMAX
            DO j = mesh%JGMIN,mesh%JGMAX
              DO i = mesh%IGMIN,mesh%IGMAX
                pout%density%data3d(i,j,k) = pin%density%data3d(mesh%IGMAX-mesh%IGMIN-j-2,i,k)
                pout%pressure%data3d(i,j,k) = pin%pressure%data3d(mesh%IGMAX-mesh%IGMIN-j-2,i,k)
                pout%velocity%data4d(i,j,k,1) = pin%velocity%data4d(mesh%IGMAX-mesh%IGMIN-j-2,i,k,2)
                pout%velocity%data4d(i,j,k,2) = pin%velocity%data4d(mesh%IGMAX-mesh%IGMIN-j-2,i,k,1)
              END DO
            END DO
          END DO
        CASE DEFAULT
          CALL mesh%Error("shear::RotateData","directions must be one of 'xy','xz' or 'yz'")
        END SELECT
      END SELECT
    END SELECT
  END SUBROUTINE rotatedata

END PROGRAM shear