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
  ! 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 :: res_xy     = 128            ! resolution in x/y-direction
  INTEGER, PARAMETER :: res_z      = 1              ! resolution in z-direction
  REAL, PARAMETER    :: box_size   = 320*gn*sigma0/(omega*omega)
  REAL, PARAMETER    :: box_height = 0.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=>null(),pvar_init=>null()
  REAL               :: sigma
  LOGICAL            :: ok
  !--------------------------------------------------------------------------!
  tap_plan(1)
 
  ! with south-north 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)
#ifndef PARALLEL
  ! store transposed initial data
  CALL sim%Physics%new_statevector(pvar_init,primitive)
  CALL rotatedata(sim%Mesh,sim%Timedisc%pvar,pvar_init,"xy")
#endif
  CALL sim%Run()
  ok = .NOT.sim%aborted
#ifndef PARALLEL
  ! store transposed result of the first run
  CALL sim%Physics%new_statevector(pvar,primitive)
  CALL rotatedata(sim%Mesh,sim%Timedisc%pvar,pvar,"xy")
 
  ! simulate west-east shear
  CALL sim%InitFosite()
  CALL makeconfig(sim, sim%config)
  CALL setattr(sim%config, "/datafile/filename", (trim(odir) // trim(ofname) // "_rotate"))
  CALL setattr(sim%config, "mesh/shearingbox", 2)
  CALL sim%Setup()
  sim%Timedisc%pvar%data1d(:) = pvar_init%data1d(:)
  SELECT TYPE(phys => sim%Physics)
  CLASS IS(physics_eulerisotherm)
    CALL phys%Convert2Conservative(sim%Timedisc%pvar,sim%Timedisc%cvar)
  CLASS DEFAULT
    CALL phys%Error("shear","only (non-)isothermal HD supported")
  END SELECT
  CALL sim%Run()
 
  ! compare results
  sigma = sqrt(sum((sim%Timedisc%pvar%data4d(:,:,:,:)-pvar%data4d(:,:,:,:))**2)/ &
                    SIZE(pvar%data4d(:,:,:,:)))
  DEALLOCATE(pvar,pvar_init)
#endif
  CALL sim%Finalize()
  DEALLOCATE(sim)
 
#ifndef PARALLEL
  ! x-y symmetry test does not work in parallel, because west-east shearing is not implemented
  tap_check_small(sigma,1e-13,"Shear x-y symmetry test")
#else
  tap_check(ok,"stoptime reached")
#endif
  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
    !------------------------------------------------------------------------!
    xmin       = -0.5*box_size
    xmax       = +0.5*box_size
    ymin       = -0.5*box_size
    ymax       = +0.5*box_size
    zmin       = -0.5*box_height
    zmax       = +0.5*box_height
 
    ! mesh settings
    mesh =>     dict(&
                "meshtype"    / midpoint, &
                "geometry"    / mgeo, &
                "shearingbox" / 1, &
                "decomposition" / (/1,-1,1/), & ! decompose along y-direction
                "inum"        / res_xy, &
                "jnum"        / res_xy, &
                "knum"        / res_z, &
                "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 &
              )
 
    ! 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, &
              "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)
    CLASS IS(statevector_eulerisotherm) ! isothermal HD
      WHERE ((abs(mesh%bcenter(:,:,:,1)).LT.0.15*box_size.AND. &
              abs(mesh%bcenter(:,:,:,2)).LT.0.15*box_size))
        p%density%data3d(:,:,:) = sigma0
      ELSEWHERE
        p%density%data3d(:,:,:) = sigma0*1e-2
      END WHERE
      SELECT CASE(mesh%shear_dir)
      CASE(1)
        p%velocity%data2d(:,2) = 0.0
        p%velocity%data4d(:,:,:,1) = mesh%Q*mesh%bcenter(:,:,:,2)*mesh%Omega
      CASE(2)
        p%velocity%data2d(:,1) = 0.0
        p%velocity%data4d(:,:,:,2) = -mesh%Q*mesh%bcenter(:,:,:,1)*mesh%Omega
      CASE DEFAULT
        CALL physics%Error("shear::InitData","shearing disabled or unsupported direction")
      END SELECT
    CLASS DEFAULT
      CALL physics%Error("shear::InitData","only (non-)isothermal HD supported")
    END SELECT
 
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_euler) ! non-isothermal HD
      p%pressure%data1d(:) = 1e-1
    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 CASE(dir)
    CASE("xy")
      SELECT TYPE(pin => pvar_in)
      CLASS IS(statevector_eulerisotherm) ! (non-)isothermal HD
        SELECT TYPE(pout => pvar_out)
        CLASS IS(statevector_eulerisotherm) ! (non-)isothermal HD
          ! rotate at middle of the field, because x' = -y in shearingsheet.
          FORALL(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
            pout%density%data3d(i,j,k) = pin%density%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 FORALL
        END SELECT
      END SELECT
      SELECT TYPE(pin => pvar_in)
      CLASS IS(statevector_euler) ! non-isothermal HD
        SELECT TYPE(pout => pvar_out)
        CLASS IS(statevector_euler) ! non-isothermal HD
          ! rotate at middle of the field, because x' = -y in shearingsheet.
          FORALL(i=mesh%IGMIN:mesh%IGMAX,j=mesh%JGMIN:mesh%JGMAX,k=mesh%KGMIN:mesh%KGMAX)
            pout%pressure%data3d(i,j,k) = pin%pressure%data3d(mesh%IGMAX-mesh%IGMIN-j-2,i,k)
          END FORALL
        END SELECT
      END SELECT
    CASE DEFAULT
      CALL mesh%Error("shear::RotateData","only 'xy'-plane is currently supported for rotating data")
    END SELECT
  END SUBROUTINE rotatedata
 
END PROGRAM shear