gauss2d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: gauss2d.f90                                                       #
!#                                                                           #
!# Copyright (C) 2006-2018                                                   #
!# 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 gauss2d
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameter
  REAL, PARAMETER     :: tsim     = 0.3      ! simulation time
  REAL, PARAMETER     :: gamma    = 1.4      ! ratio of specific heats
  REAL, PARAMETER     :: csiso    = &
                                    0.0      ! non-isothermal simulation
!                                     1.127    ! isothermal simulation
                                             !   with CSISO as sound speed
  ! initial condition (dimensionless units)
  REAL, PARAMETER     :: rho0     = 1.0      ! ambient density
  REAL, PARAMETER     :: rho1     = 100.0      ! peak density above RHO0
  REAL, PARAMETER     :: rwidth   = 0.06     ! half width of the Gaussian
  REAL, PARAMETER     :: p0       = 1.0      ! ambient pressure
  REAL, PARAMETER     :: p1       = 100.0      ! peak pressure above P0
  REAL, PARAMETER     :: pwidth   = 0.06     ! half width of the Gaussian
  REAL, PARAMETER     :: omega0   = 0.0      ! angular velocity
  REAL, PARAMETER     :: eta      = 0.0      ! dynamic viscosity (0.0 disables)
  ! location of the pulse in cylindrical coordinates
  REAL, PARAMETER     :: r0       = 0.0      ! radial position
  REAL, PARAMETER     :: z0       = 0.0      ! vertical position
  ! mesh settings
  INTEGER, PARAMETER  :: mgeo     = cartesian! geometry
  INTEGER, PARAMETER  :: xres     = 100      ! x-resolution
  INTEGER, PARAMETER  :: yres     = 100      ! y-resolution
  INTEGER, PARAMETER  :: zres     = 1        ! z-resolution
  REAL, PARAMETER     :: rmax     = 1.0      ! width of square that fits into
                                             !   computational domain
  REAL, PARAMETER     :: gpar     = 0.8      ! geometry scaling parameter
  ! output parameters
  INTEGER, PARAMETER  :: onum     = 10       ! number of output data sets
  CHARACTER(LEN=256), PARAMETER &            ! output data dir
                      :: odir     = './'
  CHARACTER(LEN=256), PARAMETER &            ! output data file name
                      :: ofname   = 'gauss2D'
  !--------------------------------------------------------------------------!
  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)

  CALL sim%Run()
  ok = .NOT.sim%aborted
  CALL sim%Finalize()
  DEALLOCATE(sim)

  tap_check(ok,"stoptime reached")
  tap_done

CONTAINS

  SUBROUTINE makeconfig(Sim, config)
    USE functions, ONLY : asinh, acosh
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite)           :: Sim
    TYPE(Dict_TYP),POINTER  :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER                 :: bc(6)
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, sources, &
                               timedisc, fluxes
    REAL                    :: x1,x2,y1,y2,z1,z2
    !------------------------------------------------------------------------!
    INTENT(INOUT)           :: sim
    !------------------------------------------------------------------------!
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(cartesian)
       x1         = -0.5
       x2         =  0.5
       y1         = -0.5
       y2         =  0.5
       z1         = -0.0
       z2         =  0.0
       bc(west)   = absorbing
       bc(east)   = absorbing
       bc(south)  = absorbing
       bc(north)  = absorbing
!        bc(WEST)   = NO_GRADIENTS
!        bc(EAST)   = NO_GRADIENTS
!        bc(SOUTH)  = NO_GRADIENTS
!        bc(NORTH)  = NO_GRADIENTS
       bc(bottom) = no_gradients
       bc(top)    = no_gradients
    CASE DEFAULT
       CALL error(sim%Physics,"InitProgram","geometry not supported for this test")
    END SELECT

    !mesh settings
    mesh => dict( &
            "meshtype"  / midpoint, &
            "geometry"  / mgeo, &
            "inum"      / xres, &
            "jnum"      / yres, &
            "knum"      / zres, &
            "xmin"      / x1,   &
            "xmax"      / x2,   &
            "ymin"      / y1,   &
            "ymax"      / y2,   &
            "zmin"      / z1,   &
            "zmax"      / z2,   &
       !     "decomposition"/ (/ 1, 2/), &
            "gparam"    / gpar)

    ! boundary conditions
    boundary => dict( &
            "western"   / bc(west), &
            "eastern"   / bc(east), &
            "southern"  / bc(south), &
            "northern"  / bc(north), &
            "bottomer"  / bc(bottom), &
            "topper"    / bc(top))

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

    ! flux calculation and reconstruction method
    fluxes => dict( &
            "order"     / linear, &
            "fluxtype"  / kt, &
!             "variables" / PRIMITIVE, &
            "variables" / conservative, &
            "limiter"   / vanleer, &
            "output/slopes" / 0)

    NULLIFY(sources)

    ! time discretization settings
    timedisc => dict( &
            "method"    / modified_euler, &
            "order"     / 3, &
            "tol_rel"   / 0.01, &
            "cfl"       / 0.4, &
            "stoptime"  / tsim, &
            "dtlimit"   / 1.0e-8, &
            "maxiter"   / 10000000)

    datafile => dict( &
            "fileformat" / vtk, &
            "filename"   / (trim(odir) // trim(ofname)), &
            "count"      / onum)

    config => dict( &
            "mesh"     / mesh, &
            "physics"  / physics, &
            "boundary" / boundary, &
            "fluxes"   / fluxes, &
            "datafile" / datafile, &
            "timedisc" / timedisc)

    IF (ASSOCIATED(sources)) &
        CALL setattr(config, "sources", sources)

  END SUBROUTINE makeconfig


  SUBROUTINE initdata(Mesh,Physics,Timedisc)
    USE physics_euler_mod, ONLY : statevector_euler
    USE physics_eulerisotherm_mod, ONLY : statevector_eulerisotherm
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(physics_base)         :: Physics
    CLASS(mesh_base)            :: Mesh
    CLASS(timedisc_base)        :: Timedisc
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX)   &
                                :: radius
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) &
                                :: posvec
    !------------------------------------------------------------------------!
    INTENT(IN)                  :: mesh,physics
    INTENT(INOUT)               :: timedisc
    !------------------------------------------------------------------------!
    IF (abs(r0).LE.tiny(r0).AND.abs(z0).LE.tiny(z0)) THEN
       ! no shift of point mass set radius and posvec to Mesh defaults
       radius(:,:,:)   = mesh%radius%bcenter(:,:,:)
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
    ELSE
       ! shifted point mass position:
       ! compute curvilinear components of shift vector
       posvec(:,:,:,1) = r0
       posvec(:,:,:,2) = z0
       CALL convert2curvilinear(mesh%geometry,mesh%bcenter,posvec,posvec)
       ! subtract the result from the position vector:
       ! this gives you the curvilinear components of all vectors pointing
       ! from the point mass to the bary center of any cell on the mesh
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - posvec(:,:,:,:)
       ! compute its absolute value
       radius(:,:,:)   = sqrt(posvec(:,:,:,1)**2+posvec(:,:,:,2)**2)
    END IF

    ! initial condition
    SELECT TYPE(pvar => timedisc%pvar)
    TYPE IS(statevector_eulerisotherm) ! isothermal HD
      ! Gaussian density pulse
      pvar%density%data3d(:,:,:)  = rho0 + rho1*exp(-log(2.0) &
               * (radius(:,:,:)/rwidth)**2)
      ! vanishing velocities
      pvar%velocity%data1d(:) = 0.0
    TYPE IS(statevector_euler) ! non-isothermal HD
      ! constant density
      pvar%density%data1d(:)  = rho0
      ! Gaussian pressure pulse
      pvar%pressure%data3d(:,:,:) = p0 + p1*exp(-log(2.0) &
               * (radius(:,:,:)/rwidth)**2)
      ! vanishing velocities
      pvar%velocity%data1d(:) = 0.0
    END SELECT

    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)

    CALL mesh%Info(" DATA-----> initial condition: 2D Gaussian pulse")

  END SUBROUTINE initdata

END PROGRAM gauss2d