ode.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sedov2d.f90                                                       #
!#                                                                           #
!# Copyright (C) 2006-2012, 2013                                             #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Björn Sperling   <sperling@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.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!

!! References:
!! \cite sedov1945 Sedov, L. I.: Unsteady motions of compressible fluids,
!!     J. Appl. Math. Mech. 9 (1945)
!! \cite sedov1959  Sedov, L. I.: Similarity and Dimensional Methods in Mechanics
!!     Academic Press Ltd., New York (1959)
!----------------------------------------------------------------------------!
PROGRAM ode_test
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 0.05     ! 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   = 3.0e-2
  ! mesh settings
  INTEGER, PARAMETER :: mgeo = cartesian   ! geometry
!  INTEGER, PARAMETER :: MGEO = CYLINDRICAL   ! geometry
!!$  INTEGER, PARAMETER :: MGEO = POLAR
!!$  INTEGER, PARAMETER :: MGEO = LOGPOLAR
!!$  INTEGER, PARAMETER :: MGEO = TANPOLAR
!!$  INTEGER, PARAMETER :: MGEO = SINHPOLAR
  INTEGER, PARAMETER :: xres = 100         ! x-resolution
  INTEGER, PARAMETER :: yres = 100         ! y-resolution
  INTEGER, PARAMETER :: zres = 1           ! z-resolution
  REAL, PARAMETER    :: rmin = 1.0e-3      ! inner radius for polar grids
  REAL, PARAMETER    :: rmax = 0.3         ! outer radius
  REAL, PARAMETER    :: gpar = 0.2         ! geometry scaling parameter     !
  ! output parameters
  INTEGER, PARAMETER :: onum = 5          ! number of output data sets
  CHARACTER(LEN=256), PARAMETER &          ! output data dir
                     :: odir = './'
  CHARACTER(LEN=256), PARAMETER &          ! output data file name
                     :: ofname = 'ode'
  !--------------------------------------------------------------------------!
  CLASS(fosite),ALLOCATABLE      :: sim
  !--------------------------------------------------------------------------!
  INTEGER, PARAMETER    :: num = 5 ! last ODE-solver
  INTEGER               :: ode, i
  CHARACTER(LEN=256),DIMENSION(NUM) :: odename
  REAL,DIMENSION(NUM)   :: time
  CHARACTER(LEN=256)    :: prefix, buffer
  !--------------------------------------------------------------------------!

  tap_plan(num)
DO ode=1,num
  ALLOCATE(sim)
  CALL sim%InitFosite()

  write(prefix,'((A),(I2.2),(A),(I2.2))') "_ode-",ode,"_geo-",mgeo
  CALL makeconfig(sim, sim%config,ode)

!  CALL PrintDict(Sim%config)
  CALL sim%Setup()
  ! set initial condition
  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)

  CALL sim%Run()

  CALL sim%ComputeRunTime()
  time(ode) = sim%run_time
  odename(ode) = sim%Timedisc%GetName()

  tap_check(.true.,"Simulation finished")
  IF(ode.NE.num) THEN
    CALL sim%Finalize()
    Deallocate(sim)
  END IF
END DO
  CALL sim%Info(repeat("*",64))
  DO i=1,num
    WRITE(buffer,'((a),(I2),(A),(F5.2))') "#", i, ". RUN, ODE solver:  " //trim(odename(i)) //", runtime: ", time(i)
    CALL sim%Info(buffer)
  END DO
  CALL sim%Info(repeat("*",64))

  CALL sim%Finalize()
  Deallocate(sim)

  tap_done

CONTAINS

 SUBROUTINE makeconfig(Sim, config, i)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Fosite)  :: Sim
    TYPE(Dict_TYP),POINTER :: config
    INTEGER           :: i
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: bc(6)
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
                               timedisc, fluxes
    REAL              :: x1,x2,y1,y2,z1,z2
    !------------------------------------------------------------------------!
    INTENT(IN)        :: i
    INTENT(INOUT)     :: sim
    !------------------------------------------------------------------------!
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(cartesian)
       x1 =-rmax
       x2 = rmax
       y1 =-rmax 
       y2 = rmax
       z1 = 0.0
       z2 = 0.0
!     CASE(POLAR)
!       x1 = RMIN
!       x2 = RMAX
!       y1 = 0.0
!       y2 = 2*PI   
     CASE(cylindrical)
       x1 = rmin
       x2 = rmax
       y1 = 0.0
       y2 = 2*pi
       z1 = 0.0
       z2 = 0.0
!    CASE(LOGPOLAR)
!       x1 = LOG(RMIN/GPAR)
!       x2 = LOG(RMAX/GPAR)
!       y1 = 0.0 
!       y2 = 2*PI       
!    CASE(TANPOLAR)
!       x1 = ATAN(RMIN/GPAR)
!       x2 = ATAN(RMAX/GPAR)
!       y1 = 0.0 
!       y2 = 2*PI       
!    CASE(SINHPOLAR)
!       x1 = RMIN/GPAR
!       x1 = LOG(x1+SQRT(1.0+x1*x1))  ! = ASINH(RMIN/GPAR))
!       x2 = RMAX/GPAR
!       x2 = LOG(x2+SQRT(1.0+x2*x2))  ! = ASINH(RMAX/GPAR))
!       y1 = 0.0 
!       y2 = 2*PI       
    CASE DEFAULT
       CALL sim%Error("InitProgram","mesh geometry not supported for 2D 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)

    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(cartesian)
       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(cylindrical)
       bc(west)  = no_gradients
       bc(east)  = no_gradients
       bc(south) = periodic
       bc(north) = periodic
       bc(bottom)= no_gradients
       bc(top)   = no_gradients
!     CASE(POLAR)
!       bc(WEST)  = NO_GRADIENTS
!       bc(EAST)  = NO_GRADIENTS
!       bc(SOUTH) = PERIODIC
!       bc(NORTH) = PERIODIC
!    CASE(LOGPOLAR)
!       bc(WEST)  = NO_GRADIENTS
!       bc(EAST)  = NO_GRADIENTS
!       bc(SOUTH) = PERIODIC
!       bc(NORTH) = PERIODIC
!    CASE(TANPOLAR)
!       bc(WEST)  = NO_GRADIENTS
!       bc(EAST)  = NO_GRADIENTS
!       bc(SOUTH) = PERIODIC
!       bc(NORTH) = PERIODIC
!    CASE(SINHPOLAR)
!       bc(WEST)  = NO_GRADIENTS
!       bc(EAST)  = NO_GRADIENTS
!       bc(SOUTH) = PERIODIC
!       bc(NORTH) = PERIODIC
    CASE DEFAULT
       CALL sim%Error("InitProgram","mesh geometry not supported for 2D Sedov explosion")
    END SELECT

    ! 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)                 ! ratio of specific heats        !

    ! flux calculation and reconstruction method
    fluxes => dict("order"     / linear, &
             "fluxtype"  / kt, &
             "variables" / conservative, &        ! vars. to use for reconstruction!
             "limiter"   / monocent, &    ! one of: minmod, monocent,...   !
             "theta"     / 1.2)          ! optional parameter for limiter !

    ! time discretization settings
    timedisc => dict( &
           "method"   / i, &
           "cfl"      / 0.3, &
           "ShowButcherTableau" / 1, &
           "stoptime" / tsim, &
           "dtlimit"  / 1.0e-15, &
           "tol_rel"  / 0.01, &
           "tol_abs"  / (/0.0,1e-3,1e-3,0.0/), &
           "output/error" / 1, &
           "maxiter"  / 1000000)

    ! initialize data input/output
    datafile => dict("fileformat" / vtk, &
!    datafile => Dict("fileformat" / GNUPLOT, "filecycles" / 0, &
               "filename"   / (trim(odir) // trim(ofname) //trim(prefix)), &
               "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) :: Physics
    CLASS(Mesh_base)    :: Mesh
    CLASS(Timedisc_base):: Timedisc
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: n
    REAL              :: P1
    !------------------------------------------------------------------------!
    INTENT(IN)        :: mesh,physics
    INTENT(INOUT)     :: timedisc
    !------------------------------------------------------------------------!
    ! isothermal modules are excluded
    SELECT TYPE (phys => physics)
    CLASS IS(physics_euler)
      ! peak pressure
      n  = 2 ! 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.
    ! pressure
    WHERE ((mesh%bccart(:,:,:,1)**2 + mesh%bccart(:,:,:,2)**2 + mesh%bccart(:,:,:,3)**2).LE.r0**2)
       ! 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: 2D Sedov explosion")

  END SUBROUTINE initdata


END PROGRAM ode_test