riemann1d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: riemann1d.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 riemann1d
  USE fosite_mod
  USE solutions
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  INTEGER, PARAMETER :: testnum = 8 ! number of available tests (do not modify)
  CHARACTER(LEN=32), PARAMETER :: teststr(testnum) = (/ &
                              "Sod shock tube                  ", &
                              "Toro test no. 2                 ", &
                              "Toro test no. 3                 ", &
                              "Toro test no. 4                 ", &
                              "Toro test no. 5                 ", &
                              "Noh problem                     ", &
                              "Isothermal shock tube           ", &
                              "Isothermal Noh problem          " /)
  !--------------------------------------------------------------------------!
  ! mesh settings
  CHARACTER(LEN=3), PARAMETER :: testdir  = "all" ! direction: x,y,z or all
  INTEGER, PARAMETER :: res  = 400       ! resolution
  ! output parameters
  INTEGER, PARAMETER :: onum = 1         ! number of output data sets
  CHARACTER(LEN=8), PARAMETER :: ofname(testnum) = (/ & ! output file name
    "sod     ","toro2   ","toro3   ","toro4   ","toro5   ","noh     ", &
    "sodiso  ","noiso   " /)
  CHARACTER(LEN=256), PARAMETER &        ! output data dir
                     :: odir = './'
  ! some global constants
  REAL               :: gamma            ! ratio of specific heats
  REAL               :: tsim             ! simulation time
  REAL               :: csiso = 0.0      ! isothermal sound speed (test no. 6)
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE     :: sim
  CLASS(marray_compound), POINTER :: pvar_exact
  CHARACTER(LEN=64)  :: verbose_tap_output
  INTEGER            :: ic,sd,dir_min,dir_max
  REAL, DIMENSION(:,:), ALLOCATABLE :: sigma
  REAL, PARAMETER    :: sigma_tol(testnum) &
                        = (/ 1.8e-2, 1.4e-2, 7.9e+0, 3.5e+1, 4.3e+0, 7.9e-2, 1.0e+0, 1.0e+0 /)
  !--------------------------------------------------------------------------!

  ! check whether we perform tests in all directions
  SELECT CASE (trim(testdir))
  CASE("all")
    dir_min = 1
    dir_max = 3
  CASE("y")
    dir_min = 2
    dir_max = 2
  CASE("z")
    dir_min = 3
    dir_max = 3
  CASE DEFAULT
    ! x-direction is also the default
    dir_min = 1
    dir_max = 1
  END SELECT
  
  tap_plan(testnum*(dir_max-dir_min+1))

  ALLOCATE(sigma(testnum,dir_min:dir_max))
  ! loop over all tests
  DO ic=1,testnum
    ! loop over selected directions
    DO sd=dir_min,dir_max
    
      ! initialize Fosite
      ALLOCATE(sim)
      CALL sim%InitFosite()
      
      ! get configuration for the particular test and direction
      CALL makeconfig(sim, sim%config,ic,sd)
!       CALL PrintDict(config)

      ! setup the simulation
      CALL sim%Setup()

      ! create state vector for exact solution
      CALL sim%Physics%new_statevector(pvar_exact,primitive)

      ! set initial conditions and run the simulation
      CALL run(sim,pvar_exact,ic,sd)

      IF (sim%aborted) THEN
        sigma(ic,sd) = huge(1.0)
      ELSE
        IF (ASSOCIATED(sim%Timedisc%solution)) THEN
          sigma(ic,sd) = sqrt(sum((sim%Timedisc%pvar%data1d(:) &
                          -sim%Timedisc%solution%data1d(:))**2)/SIZE(sim%Timedisc%pvar%data1d))
        ELSE
          sigma(ic,sd) = 0.0
        END IF
      END IF

      CALL pvar_exact%Destroy()
      CALL sim%Finalize()
      DEALLOCATE(sim)

    END DO
  END DO

  ! check results
  ! loop over all tests
  DO ic=1,testnum
    ! loop over selected directions
    DO sd=dir_min,dir_max
      WRITE (verbose_tap_output,'(A,ES8.2)') trim(teststr(ic)) // " / " // achar(119+sd) &
        // " sigma = ", sigma(ic,sd)
! This line is long if expanded. So we can't indent it or it will be cropped.
tap_check_small(sigma(ic,sd),sigma_tol(ic),trim(verbose_tap_output))
    END DO
  END DO

tap_done

CONTAINS

  SUBROUTINE makeconfig(Sim,config,ic,dir)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Fosite)  :: Sim
    TYPE(Dict_TYP),POINTER :: config
    INTEGER           :: ic,dir
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
                               timedisc, fluxes
    INTEGER           :: xres,yres,zres
    REAL              :: xmax,ymax,zmax
    !------------------------------------------------------------------------!
    INTENT(INOUT)     :: sim
    INTENT(IN)        :: ic,dir
    !------------------------------------------------------------------------!
    ! mesh settings depends on the direction selected above
    SELECT CASE(dir)
    CASE(1)
      xmax = 1.0
      ymax = 0.0
      zmax = 0.0
      xres = res
      yres = 1
      zres = 1
    CASE(2)
      xmax = 0.0
      ymax = 1.0
      zmax = 0.0
      xres = 1
      yres = res
      zres = 1
    CASE(3)
      xmax = 0.0
      ymax = 0.0
      zmax = 1.0
      xres = 1
      yres = 1
      zres = res
    CASE DEFAULT
      CALL sim%Error("riemann1d::InitData","direction must be one of 1,2,3")
    END SELECT

    ! mesh settings
    mesh => dict("meshtype" / midpoint, &
           "geometry" / cartesian, &
               "inum" / xres, &            ! resolution in x and            !
               "jnum" / yres, &            !   y direction                  !
               "knum" / zres, &            !   z direction                  !
               "xmin" / (0.0), &
               "xmax" / xmax, &
               "ymin" / (-0.0), &
               "ymax" / ymax, &
               "zmin" / (0.0), &
               "zmax" / zmax  &
               )
    
    SELECT CASE(ic)
    CASE(1)
       gamma  = 1.4
       tsim   = 0.25
    CASE(2)
       gamma  = 1.4
       tsim   = 0.15
    CASE(3)
       gamma  = 1.4
       tsim   = 0.012
    CASE(4)
       gamma  = 1.4
       tsim   = 0.035
    CASE(5)
       gamma  = 1.4
       tsim   = 0.012
    CASE(6)
       gamma  = 5./3.
       tsim   = 1.0
    CASE(7)
       csiso  = 1.0
       tsim   = 0.25
    CASE(8)
       csiso  = 1.0
       tsim   = 0.5
    CASE DEFAULT
       CALL sim%Mesh%Error("riemann1d::InitData", "Test problem number should be 1,2,3,4,5,6,7 or 8")
    END SELECT


    ! physics settings
    IF (csiso.NE.0.0) THEN
       physics => dict("problem" / euler_isotherm, &
                 "cs"      / csiso)
    ELSE   
       physics => dict("problem" / euler, &
                 "gamma"   / gamma)         ! ratio of specific heats        !
    END IF

    ! flux calculation and reconstruction method
    fluxes => dict(&
             "fluxtype"  / kt, &
             "order"     / linear, &
             "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"   / modified_euler, &
           "order"    / 3, &
           "cfl"      / 0.4, &
           "stoptime" / tsim, &
           "dtlimit"  / 1.0e-10, &
!           "output/xmomentum" / 1, &
!           "output/ymomentum" / 1, &
!           "output/zmomentum" / 1, &
!           "output/energy" / 1, &
!           "output/rhs" / 1 &
           "output/solution" / 0, &
           "maxiter"  / 100000 &
           )

    ! compute and output exact solution in non-isothermal simulations
    IF (csiso.EQ.0.0) CALL setattr(timedisc,"output/solution",1)

    ! boundary conditions
    boundary => dict( &
           "western"  / no_gradients, &
           "eastern"  / no_gradients, &
           "southern" / no_gradients, &
           "northern" / no_gradients, &
           "bottomer" / no_gradients, &
           "topper"   / no_gradients)

    ! initialize data input/output
    datafile => dict( &
           "fileformat" / vtk, &
           "filename"   / (trim(odir) // trim(ofname(ic)) // "-" // achar(119+dir)), &
!           "filecycles" / 0, &
           "count"      / onum)

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

  SUBROUTINE run(this,pvar,ic,dir)
    USE solutions
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite), INTENT(INOUT) :: this
    CLASS(marray_compound), POINTER :: pvar
    INTEGER, INTENT(IN)          :: ic,dir
    !------------------------------------------------------------------------!
    ! Local variable declaration
    REAL    :: x0, rho_l, rho_r, u_l, u_r, p_l, p_r
    REAL, DIMENSION(:), ALLOCATABLE :: x
    !------------------------------------------------------------------------!

    SELECT CASE(ic)
    CASE(1) ! Sod shock tube (Toro's test no. 1)
       x0    = 0.5
       rho_l = 1.0
       rho_r = 0.125
       u_l   = 0.0
       u_r   = 0.0
       p_l   = 1.0
       p_r   = 0.1
    CASE(2) ! Toro's test no. 2
       x0    = 0.5
       rho_l = 1.0
       rho_r = 1.0
       u_l   = -2.0
       u_r   = 2.0
       p_l   = 0.4
       p_r   = 0.4
    CASE(3) ! Toro's test no. 3
       x0    = 0.5
       rho_l = 1.0
       rho_r = 1.0
       u_l   = 0.0
       u_r   = 0.0
       p_l   = 1000.
       p_r   = 0.01
    CASE(4) ! Toro's test no. 4
       x0    = 0.4
       rho_l = 5.99924
       rho_r = 5.99924
       u_l   = 19.5975
       u_r   = -6.19633
       p_l   = 460.894
       p_r   = 46.0950
    CASE(5) ! Toro's test no. 5
       x0    = 0.8
       rho_l = 1.0
       rho_r = 1.0
       u_l   = -19.5975
       u_r   = -19.5975
       p_l   = 1000.0
       p_r   = 0.01
    CASE(6) ! Noh problem
       x0    = 0.5
       rho_l = 1.0
       rho_r = 1.0
       u_l   = 1.0
       u_r   = -1.0
       p_l   = 1.0e-05
       p_r   = 1.0e-05
    CASE(7) ! isothermal shock tube
       x0    = 0.5
       rho_l = 1.0
       rho_r = 0.125
       u_l   = 0.0
       u_r   = 0.0
    CASE(8) ! isothermal Noh problem
       x0    = 0.5
       rho_l = 1.0
       rho_r = 1.0
       u_l   = 1.0
       u_r   = -1.0
    CASE DEFAULT
       CALL this%Error("InitData", "Test problem should be 1,2,3,4,5,6,7 or 8")
    END SELECT

    ! initial condition
    ! always set density and first velocity
    SELECT TYPE(pvar => this%Timedisc%pvar)
    CLASS IS(statevector_eulerisotherm)
      ! set all velocities to 0
      pvar%velocity%data1d(:) = 0.0
      WHERE (this%Mesh%bcenter(:,:,:,dir).LT.x0)
        pvar%density%data3d(:,:,:)    = rho_l
        pvar%velocity%data4d(:,:,:,1) = u_l
      ELSEWHERE
        pvar%density%data3d(:,:,:)    = rho_r
        pvar%velocity%data4d(:,:,:,1) = u_r
      END WHERE
    END SELECT
    
    ! set pressure for non-isothermal HD
    SELECT TYPE(pvar => this%Timedisc%pvar)
    TYPE IS(statevector_euler)
      WHERE (this%Mesh%bcenter(:,:,:,dir).LT.x0)
        pvar%pressure%data3d(:,:,:)   = p_l
      ELSEWHERE
        pvar%pressure%data3d(:,:,:)   = p_r
      END WHERE
    END SELECT
   
    CALL this%Physics%Convert2Conservative(this%Timedisc%pvar,this%Timedisc%cvar)
    CALL this%Info(" DATA-----> initial condition: " &
      // "1D Riemann problem: " // trim(teststr(ic)))

    IF (ASSOCIATED(this%Timedisc%solution)) THEN
      ! store initial condition in exact solution array
      this%Timedisc%solution = this%Timedisc%pvar

      ! set 1D coordinate array for exact solution
      ALLOCATE(x(SIZE(this%Mesh%cart%data3d(:,2,dir),1)), &
        source=this%Mesh%cart%data3d(:,2,dir))

      DO WHILE((this%Timedisc%maxiter.LE.0).OR.(this%iter.LE.this%Timedisc%maxiter))
        SELECT TYPE(phys => this%Physics)
        TYPE IS (physics_euler)
          ! compute exact solution of the 1D Riemann problem (only non-isothermal HD)
          CALL riemann(x0,phys%gamma,rho_l,u_l,p_l,rho_r,u_r,p_r, &
                      this%Timedisc%time,x,this%Timedisc%solution%data2d(:,:))
        END SELECT
          IF(this%Step()) EXIT
       END DO
    ELSE
       CALL this%Run()
    END IF
  END SUBROUTINE run

END PROGRAM riemann1d