gauss3d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: gauss3d.f90                                                       #
!#                                                                           #
!# Copyright (C) 2006-2014                                                   #
!# 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 gauss3d
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameter
  REAL, PARAMETER     :: tsim     = 0.5      ! 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     = 1.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       = 1.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 initial pulse (cartesian coordinates)
  REAL, PARAMETER     :: x0       = 0.0      ! x-position
  REAL, PARAMETER     :: y0       = 0.0      ! y-position
  REAL, PARAMETER     :: z0       = 0.0      ! z-position
  ! mesh settings
  INTEGER, PARAMETER  :: mgeo     = cartesian! geometry
!   INTEGER, PARAMETER  :: MGEO     = CYLINDRICAL
!   INTEGER, PARAMETER  :: MGEO     = SPHERICAL
  INTEGER, PARAMETER  :: res      = 30       ! x-/y-/z-resolution
  REAL, PARAMETER     :: rmax     = 0.5      ! maximal radial extent
  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   = 'gauss3D'
  !--------------------------------------------------------------------------!
  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)
!     SELECT TYPE(pvar => Sim%Fluxes%prim)
!     TYPE IS(statevector_eulerisotherm) ! isothermal HD
!       PRINT *,pvar%density%data1d(1)
!     CLASS DEFAULT
!       PRINT *, "ERROR: wrong state vector"
!     END SELECT
  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),xres,yres,zres
    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         = -rmax
       x2         =  rmax
       y1         = -rmax
       y2         =  rmax
       z1         = -rmax
       z2         =  rmax
       xres       = res
       yres       = res
       zres       = res
!        bc(WEST)   = NO_GRADIENTS
!        bc(EAST)   = NO_GRADIENTS
!        bc(SOUTH)  = NO_GRADIENTS
!        bc(NORTH)  = NO_GRADIENTS
!        bc(BOTTOM) = NO_GRADIENTS
!        bc(TOP)    = NO_GRADIENTS
       bc(west)   = absorbing
       bc(east)   = absorbing
       bc(south)  = absorbing
       bc(north)  = absorbing
       bc(bottom) = absorbing
       bc(top)    = absorbing
    CASE(cylindrical)
       x1 =  0.01
       x2 =  rmax
       y1 =  0.0
       y2 = 2*pi
       z1 = -rmax
       z2 =  rmax
       xres       = res / 2
       yres       = res / 2
       zres       = res
       bc(west)   = axis
       bc(east)   = absorbing
       bc(south)  = periodic
       bc(north)  = periodic
       bc(bottom) = absorbing
       bc(top)    = absorbing
    CASE(spherical)
       x1 = 0.01
       x2 = rmax
       y1 = 0.0
       y2 = pi
       z1 = 0.0
       z2 = 2*pi
       xres       = res / 2
       yres       = res / 4
       zres       = res / 2
       bc(west)  = axis
       bc(east)  = absorbing
       bc(south) = axis
       bc(north) = axis
       bc(bottom)= periodic
       bc(top)   = periodic
    CASE DEFAULT
       CALL sim%Error("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, 1, 2/), &
            "output/radius" / 1, &
            "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, &
            "limiter"   / vanleer, &
            "output/slopes" / 0)
 
    NULLIFY(sources)
 
    ! time discretization settings
    timedisc => dict( &
            "method"    / modified_euler, &
            "order"     / 3, &
            "tol_rel"   / 10.0, & ! > 1 disables adaptive step size
            "cfl"       / 0.4, &
            "stoptime"  / tsim, &
            "dtlimit"   / 1.0e-8, &
            "maxiter"   / 10000000)
 
    datafile => dict( &
            "fileformat" / xdmf, &
            "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 geometry_base_mod
    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 (any((/abs(x0),abs(y0),abs(z0)/).GT.tiny(z0))) THEN
       ! shifted point mass position:
       ! compute curvilinear components of shift vector
       posvec(:,:,:,1) = x0
       posvec(:,:,:,2) = y0
       posvec(:,:,:,3) = z0
       CALL mesh%geometry%Convert2Curvilinear(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+posvec(:,:,:,3)**2)
    ELSE
       ! no shift of point mass set radius and posvec to Mesh defaults
       radius(:,:,:)   = mesh%radius%bcenter(:,:,:)
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
    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
    CLASS DEFAULT
      CALL physics%Error("gauss3d::InitData","unsupported physics")
    END SELECT
 
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
 
    CALL mesh%Info(" DATA-----> initial condition: 3D Gaussian pulse")
 
  END SUBROUTINE initdata
 
END PROGRAM gauss3d