sedov3d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sedov3d.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 sedov3d
  USE fosite_mod
  USE solutions
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 0.05        ! simulation stop time
  REAL, PARAMETER    :: gamma   = 1.4         ! ratio of specific heats
  REAL, PARAMETER    :: r       = 1.033       ! scaling parameter: 1.0 for 2D, 1.033 for 3D
  ! initial condition (dimensionless units)
  REAL, PARAMETER    :: rho0    = 1.0         ! ambient density
  REAL, PARAMETER    :: p0      = 1.0e-05     ! ambient pressure
  REAL, PARAMETER    :: e1      = 1.0         ! initial energy input
  ! Spatial with of the initial pulse should be at least 5 cells;
  ! if you wish to compare the results on different grids
  ! R0 should be of the same order
  REAL, PARAMETER    :: r0      = 5.0e-2
  ! mesh settings
!  INTEGER, PARAMETER :: MGEO    = LOGSPHERICAL   ! geometry
  INTEGER, PARAMETER :: mgeo    = spherical   ! geometry
!  INTEGER, PARAMETER :: MGEO    = CARTESIAN   ! geometry
  INTEGER, PARAMETER :: xres    = 50          ! x-resolution
  INTEGER, PARAMETER :: yres    = 2           ! y-resolution
  INTEGER, PARAMETER :: zres    = 2           ! z-resolution
  REAL, PARAMETER    :: rmax    = 0.4         ! outer radius of comput. domain
  REAL, PARAMETER    :: gpar    = 0.2         ! 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  = 'sedov3d'
  !-------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE  :: sim
  !-------------------------------------------------------------------------!
  REAL, DIMENSION(:),ALLOCATABLE  :: pvar_diff
#ifdef PARALLEL
  REAL, DIMENSION(:), POINTER     :: pvar,pvar_all,radius,radius_all
  INTEGER                         :: err
#endif
  INTEGER :: i
  REAL    :: rt,rshock
  !-------------------------------------------------------------------------!
  ALLOCATE(sim,pvar_diff(1:xres))
  CALL sim%InitFosite()
#ifdef PARALLEL
  IF(sim%GetRank().EQ.0) &
#endif
  tap_plan(1)

  CALL makeconfig(sim, sim%config)
  CALL sim%Setup()
  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)

  CALL sim%Run()
  
#ifdef PARALLEL
  ALLOCATE(pvar(sim%Mesh%IMIN:sim%Mesh%IMAX),pvar_all(sim%GetNumProcs()*(sim%Mesh%IMAX-sim%Mesh%IMIN+1)))
  ALLOCATE(radius(sim%Mesh%IMIN:sim%Mesh%IMAX),radius_all(sim%GetNumProcs()*(sim%Mesh%IMAX-sim%Mesh%IMIN+1)))
  pvar(:) = sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,1,1,1)
  radius(:)= sim%Mesh%radius%center(sim%Mesh%IMIN:sim%Mesh%IMAX,1,1)
 
  ! Compare results with analytical solution. Check only for correct shock velocity
  ! even if the full analytical solution is implemented

!  CALL sedov(GAMMA,E1,RHO0,P0,TSIM,3, SIM%Mesh%radius%bcenter(1:XRES,1,1),pvar)
  
  !The shock location is where the density gradient is greatest
  CALL mpi_gather(pvar,int(sim%Mesh%IMAX-sim%Mesh%IMIN+1),mpi_double, &
    pvar_all, int(sim%Mesh%IMAX-sim%Mesh%IMIN+1),mpi_double,0,mpi_comm_world,err)
  CALL mpi_gather(radius,int(sim%Mesh%IMAX-sim%Mesh%IMIN+1),mpi_double, &
    radius_all, int(sim%Mesh%IMAX-sim%Mesh%IMIN+1),mpi_double,0,mpi_comm_world,err)

 IF(sim%GetRank().EQ.0) THEN 
  DO i=1,xres-1
    pvar_diff(i) = pvar_all(i)-pvar_all(i+1)
  END DO
  rshock = radius_all(maxloc(pvar_diff,dim=1))
  
  !analytical shock position
  rt = r*(e1*tsim**2/rho0)**0.2

  !Check whether analytical solution is within +/- one cell of simulated shock
  tap_check((rt.LT.rshock+sim%Mesh%dx).AND.(rt.GT.rshock-sim%Mesh%dx),"Shock velocity correct")
END IF
#else
  DO i=1,xres
      pvar_diff(i) = sim%Timedisc%pvar%data4d(i,1,1,1)-sim%Timedisc%pvar%data4d(i+1,1,1,1)
  END DO
  rshock = sim%Mesh%radius%center(maxloc(pvar_diff,dim=1),1,1)
  
  !analytical shock position
  rt = r*(e1*tsim**2/rho0)**0.2

  tap_check((rt.LT.rshock+sim%Mesh%dx).AND.(rt.GT.rshock-sim%Mesh%dx),"Shock velocity correct")
#endif
  CALL sim%Finalize()
  tap_done
  DEALLOCATE(sim,pvar_diff)

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, &
                               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.5
       z2 =  0.5
       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(spherical)
       x1 = 0.0
       x2 = rmax
       y1 = 0.0
       y2 = pi
       z1 = 0.0
       z2 = 2*pi
       bc(west)  = reflecting
       bc(east)  = no_gradients
       bc(south) = axis
       bc(north) = axis
       bc(bottom)= periodic
       bc(top)   = periodic
    CASE(logspherical)
       x1 = log(1e-3/gpar)
       x2 = log(rmax/gpar)
       y1 = 0.0
       y2 = pi
       z1 = 0.0
       z2 = 2*pi
       bc(west)  = reflecting
       bc(east)  = no_gradients
       bc(south) = axis
       bc(north) = axis
       bc(bottom)= periodic
       bc(top)   = periodic
    CASE(cylindrical)
       x1 = -rmax
       x2 = rmax
       y1 = 0.0
       y2 = 2.0*pi
       z1 = 0.0
       z2 = rmax
       bc(west)  = no_gradients    !default: ABSORBING
       bc(east)  = no_gradients    !default: ABSORBING
       bc(south) = no_gradients    !default: AXIS
       bc(north) = no_gradients    !default: ABSORBING
       bc(bottom)= no_gradients
       bc(top)   = no_gradients
!     CASE(OBLATE_SPHEROIDAL)       !TODO anderen Meshs zu 3D konvertieren
!        x1 = 0.0
!        x2 = Acosh(RMAX/GPAR)
!        y1 = -0.5*PI
!        y2 = 0.5*PI
!        bc(WEST)  = FOLDED
!        bc(EAST)  = NO_GRADIENTS
!        bc(SOUTH) = AXIS
!        bc(NORTH) = AXIS
!     CASE(TANCYLINDRICAL)
!        x1 = ATAN(-RMAX/GPAR)
!        x2 = ATAN(RMAX/GPAR)
!        y1 = 0.0
!        y2 = RMAX
!        bc(WEST)  = NO_GRADIENTS
!        bc(EAST)  = NO_GRADIENTS
!        bc(SOUTH) = AXIS
!        bc(NORTH) = NO_GRADIENTS
!     CASE(SINHSPHERICAL)
!        x1 = 0.0
!        x2 = Asinh(RMAX/GPAR)
!        y1 = 0.0
!        y2 = PI
!        bc(WEST)  = REFLECTING
!        bc(EAST)  = NO_GRADIENTS
!        bc(SOUTH) = AXIS
!        bc(NORTH) = AXIS
    CASE DEFAULT
       CALL sim%Physics%Error("InitProgram","geometry not supported for 3D Sedov explosion")
    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, &
              "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
    physics => dict( &
              "problem"     / euler, &
              "gamma"       / gamma)

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

    ! time discretization settings
    timedisc => dict( &
      "method"      / dormand_prince,               &
              "order"       / 5,                            &
              "cfl"         / 0.4,                          &
              "stoptime"    / tsim,                         &
              "dtlimit"     / 1.0e-13,                      &
              "tol_rel"     / 0.01,                         &
              "tol_abs"     / (/1e-5,1e-5,1e-5,1e-5,1e-5/), &
              "maxiter"     / 1000000                       )

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

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

  SUBROUTINE initdata(Mesh,Physics,Timedisc)
    USE physics_euler_mod, ONLY : physics_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(physics_base),  INTENT(IN)    :: Physics
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER                             :: n
    REAL                                :: P1
    !------------------------------------------------------------------------!
    ! isothermal modules are excluded
    SELECT TYPE (phys => physics)
    CLASS IS(physics_euler)
      ! peak pressure
      n  = 3 ! 3 for 3D
      p1 = 3.*(phys%gamma - 1.0)*e1 / ((n + 1)*pi*r0**n)
    CLASS DEFAULT
      ! abort
      CALL phys%Error("InitData","physics not supported")
    END SELECT

    ! uniform density
    timedisc%pvar%data4d(:,:,:,physics%DENSITY)   = rho0
    ! vanishing initial velocities
    timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.
    timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.
    timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = 0.
    ! pressure
    WHERE (mesh%radius%bcenter(:,:,:).LE.r0)
       ! behind the shock front
       timedisc%pvar%data4d(:,:,:,physics%PRESSURE)  = p1
    ELSEWHERE
       ! in front of the shock front (ambient medium)
       timedisc%pvar%data4d(:,:,:,physics%PRESSURE)  = p0
    END WHERE

    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
    CALL mesh%Info(" DATA-----> initial condition: 3D Sedov explosion")

  END SUBROUTINE initdata

END PROGRAM sedov3d