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    = 1.0         ! simulation stop time
  REAL, PARAMETER    :: gamma   = 1.4         ! ratio of specific heats
  ! 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    = 200          ! x-resolution
  INTEGER, PARAMETER :: yres    = 6           ! y-resolution
  INTEGER, PARAMETER :: zres    = 1           ! z-resolution
  REAL, PARAMETER    :: rmax    = 2.0         ! 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 :: sigma, sum_numer, sum_denom
  INTEGER :: n,den,vel,pre
  LOGICAL :: ok
  !-------------------------------------------------------------------------!
  ALLOCATE(sim)
  CALL sim%InitFosite()
 
#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_plan(4)
#ifdef PARALLEL
  END IF
#endif
 
  CALL makeconfig(sim, sim%config)
  CALL sim%Setup()
  ALLOCATE(sum_numer(sim%Physics%VNUM),sum_denom(sim%Physics%VNUM),sigma(sim%Physics%VNUM))
  sigma(:) = 0.0
 
  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
  CALL run(sim%Mesh, sim%Physics, sim%Timedisc)
  ok = .NOT.sim%aborted
 
  ! compare with exact solution if requested
  IF (ASSOCIATED(sim%Timedisc%solution)) THEN
    DO n=1,sim%Physics%VNUM
      ! use L1 norm to estimate the deviation from the exact solution:
      !   Σ |pvar - pvar_exact| / Σ |pvar_exact|
      sum_numer(n) = sum(abs(sim%Timedisc%pvar%data2d(:,n)-sim%Timedisc%solution%data2d(:,n)), &
                         mask=sim%Mesh%without_ghost_zones%mask1d)
      sum_denom(n) = sum(abs(sim%Timedisc%solution%data2d(:,n)),mask=sim%Mesh%without_ghost_zones%mask1d)
    END DO
#ifdef PARALLEL
    IF (sim%GetRank().GT.0) THEN
      CALL mpi_reduce(sum_numer,sum_numer,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
      CALL mpi_reduce(sum_denom,sum_denom,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
    ELSE
      CALL mpi_reduce(mpi_in_place,sum_numer,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
      CALL mpi_reduce(mpi_in_place,sum_denom,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
    END IF
#endif
    WHERE (abs(sum_denom(:)).GT.2*tiny(sum_denom))
      sigma(:) = sum_numer(:) / sum_denom(:)
    END WHERE
  END IF
 
  den = sim%Physics%DENSITY
  vel = sim%Physics%XVELOCITY
  pre = sim%Physics%PRESSURE
 
#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_check(ok,"stoptime reached")
! These lines are very long if expanded. So we can't indent it or it will be cropped.
tap_check_small(sigma(den),5.0e-02,"density deviation < 5%")
tap_check_small(sigma(vel),4.0e-02,"radial velocity deviation < 4%")
tap_check_small(sigma(pre),5.0e-02,"pressure deviation < 5%")
tap_done
!   PRINT *,sigma(:)
#ifdef PARALLEL
  END IF
#endif
 
  CALL sim%Finalize()
  DEALLOCATE(sim,sigma,sum_numer,sum_denom)
 
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 = 0.0 !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"     / vanleer)
 
    ! 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/), &
              "output/solution" / 1, &
              "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)
      SELECT TYPE (pvar => timedisc%pvar)
      CLASS IS(statevector_euler)
        ! uniform density everywhere
        pvar%density%data1d(:)  = rho0
        ! vanishing initial velocities
        pvar%velocity%data1d(:) = 0.0
        ! set initial peak pressure P1 inside sphere with radius R0 centered on the origin
        n  = 3 ! 3 for 3D
        p1 = 3.*(phys%gamma - 1.0)*e1 / ((n + 1)*pi*r0**n)
        WHERE (mesh%radius%bcenter(:,:,:).LE.r0)
          ! behind the shock front
          pvar%pressure%data3d(:,:,:)  = p1
        ELSEWHERE
          ! in front of the shock front (ambient medium)
          pvar%pressure%data3d(:,:,:)  = p0
        END WHERE
      CLASS DEFAULT
        ! abort
        CALL phys%Error("sedov3d:InitData","statevector must be of class euler")
      END SELECT
    CLASS DEFAULT
      ! abort
      CALL phys%Error("sedov3d:InitData","physics not supported")
    END SELECT
 
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
    CALL mesh%Info(" DATA-----> initial condition: 3D Sedov explosion")
 
  END SUBROUTINE initdata
 
  SUBROUTINE run(Mesh,Physics,Timedisc)
    USE sedov_mod
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(IN)    :: Physics
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    !------------------------------------------------------------------------!
    TYPE(sedov_typ), ALLOCATABLE :: sedov
    REAL :: T0,vr
    TYPE(marray_base) :: er
    INTEGER :: n,m,i,j,k
    !------------------------------------------------------------------------!
    IF (ASSOCIATED(timedisc%solution)) THEN
      ! initialize sedov class
      ! parameters: ndim, gamma, w, rho0, E1, P0
      sedov = sedov_typ(3,gamma,0.0,rho0,e1,p0)
#ifdef PARALLEL
      IF (sim%GetRank().EQ.0) &
#endif
      CALL sedov%PrintConfiguration()
      ! store initial condition in exact solution array
      timedisc%solution = timedisc%pvar
 
      ! compute initial time using R0 as the initial position of the shock
      t0 = sqrt((e1/rho0)*(r0/sedov%alpha)**5)
 
      ! set radial unit vector to compute full 3D radial velocities in any geometry
      er = marray_base(3)
      DO n=1,3
        er%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,n) = &
            mesh%posvec%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,n) &
          / mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)
      END DO
      DO WHILE((timedisc%maxiter.LE.0).OR.(sim%iter.LE.timedisc%maxiter))
        SELECT TYPE(solution => timedisc%solution)
        TYPE IS (statevector_euler)
          ! compute exact solution of the spherical Sedov explosion problem
          DO k=mesh%KMIN,mesh%KMAX
            DO j=mesh%JMIN,mesh%JMAX
              DO i=mesh%IMIN,mesh%IMAX
                CALL sedov%ComputeSolution(timedisc%time+t0,mesh%radius%bcenter(i,j,k), &
                  solution%density%data3d(i,j,k),vr,solution%pressure%data3d(i,j,k))
                  ! set the velocities
                  m = 1
                  DO n=1,3
                    IF (btest(mesh%VECTOR_COMPONENTS,n-1)) THEN
                      solution%velocity%data4d(i,j,k,m) = vr * er%data4d(i,j,k,n)
                      m = m + 1
                    END IF
                  END DO
              END DO
            END DO
          END DO
        END SELECT
        IF(sim%Step()) EXIT
      END DO
      CALL sedov%Destroy()
    ELSE
      CALL sim%Run()
    END IF
  END SUBROUTINE run
 
END PROGRAM sedov3d