foreseer_test_shock_tube Program

program foreseer_test_shock_tube
!< FORESEER test: shock tube tester, 1D Euler equation.

use flap, only : command_line_interface
use foodie, only : tvd_runge_kutta_integrator
use foreseer, only : conservative_compressible, primitive_compressible,         &
                     conservative_to_primitive_compressible, primitive_to_conservative_compressible, &
                     eos_compressible
use foreseer_euler_1d, only : euler_1d
use penf, only : cton, FR8P, I4P, R8P, str
use vecfor, only : ex, vector

implicit none
integer(I4P)                     :: weno_order                !< WENO reconstruction order.
character(len=:), allocatable    :: weno_variables            !< Variables set on which WENO reconstruction is done.
type(tvd_runge_kutta_integrator) :: rk_integrator             !< Runge-Kutta integrator.
integer(I4P)                     :: rk_stages_number          !< Runge-Kutta stages number.
type(euler_1d), allocatable      :: rk_stage(:)               !< Runge-Kutta stages.
real(R8P)                        :: dt                        !< Time step.
real(R8P)                        :: t                         !< Time.
integer(I4P)                     :: step                      !< Time steps counter.
type(euler_1d)                   :: domain                    !< Domain of Euler equations.
real(R8P)                        :: CFL                       !< CFL value.
character(3)                     :: BC_L                      !< Left boundary condition type.
character(3)                     :: BC_R                      !< Right boundary condition type.
integer(I4P)                     :: Ni                        !< Number of grid cells.
real(R8P)                        :: Dx                        !< Space step discretization.
real(R8P), allocatable           :: x(:)                      !< Cell center x-abscissa values.
integer(I4P)                     :: steps_max                 !< Maximum number of time steps.
real(R8P)                        :: t_max                     !< Maximum integration time.
character(99), allocatable       :: riemann_solver_schemes(:) !< Riemann Problem solver scheme(s).
character(99)                    :: s_scheme                  !< Space integration scheme.
character(99)                    :: t_scheme                  !< Time integration scheme.
logical                          :: results                   !< Flag for activating results saving.
logical                          :: time_serie                !< Flag for activating time serie-results saving.
logical                          :: verbose                   !< Flag for activating more verbose output.
integer(I4P)                     :: s                         !< Schemes counter.

call parse_command_line_interface
do s=1, size(riemann_solver_schemes, dim=1)
   if (verbose) print "(A)", 'Use Riemann Problem solver "'//trim(adjustl(riemann_solver_schemes(s)))//'"'
   call initialize(riemann_solver_scheme=riemann_solver_schemes(s))
   call save_time_serie(filename='euler_1D-'//&
                                 trim(adjustl(s_scheme))//'-'//&
                                 trim(adjustl(t_scheme))//'-'//&
                                 trim(adjustl(riemann_solver_schemes(s)))//'.dat', t=t)
   step = 0
   time_loop: do
      step = step + 1
      dt = domain%dt(steps_max=steps_max, t_max=t_max, t=t, CFL=CFL)
      call rk_integrator%integrate(U=domain, stage=rk_stage, dt=dt, t=t)
      t = t + dt
      call save_time_serie(t=t)
      if (verbose) print "(A)", 'step = '//str(n=step)//', time step = '//str(n=dt)//', time = '//str(n=t)
      if ((t == t_max).or.(step == steps_max)) exit time_loop
   enddo time_loop
enddo

contains
   subroutine initialize(riemann_solver_scheme)
   !< Initialize the test.
   character(*), intent(in)                  :: riemann_solver_scheme !< Riemann Problem solver scheme.
   type(primitive_compressible), allocatable :: initial_state(:)      !< Initial state of primitive variables.
   integer(I4P)                              :: i                     !< Space counter.

   if (allocated(rk_stage)) deallocate(rk_stage) ; allocate(rk_stage(1:rk_stages_number))
   call rk_integrator%init(stages=rk_stages_number)
   t = 0._R8P
   if (allocated(x)) deallocate(x) ; allocate(x(1:Ni))
   if (allocated(initial_state)) deallocate(initial_state) ; allocate(initial_state(1:Ni))
   Dx = 1._R8P / Ni
   ! Sod's problem
   BC_L = 'TRA'
   BC_R = 'TRA'
   do i=1, Ni / 2
      x(i) = Dx * i - 0.5_R8P * Dx
      initial_state(i)%density  = 1._R8P
      initial_state(i)%velocity = 0._R8P
      initial_state(i)%pressure = 1._R8P
   enddo
   do i=Ni / 2 + 1, Ni
      x(i) = Dx * i - 0.5_R8P * Dx
      initial_state(i)%density  = 0.125_R8P
      initial_state(i)%velocity = 0._R8P
      initial_state(i)%pressure = 0.1_R8P
   enddo
   call domain%initialize(Ni=Ni, Dx=Dx,                                         &
                          BC_L=BC_L, BC_R=BC_R,                                 &
                          initial_state=initial_state,                          &
                          eos=eos_compressible(cp=1040.004_R8P, cv=742.86_R8P), &
                          weno_order=weno_order,                                &
                          weno_variables=weno_variables,                        &
                          riemann_solver_scheme=riemann_solver_scheme)
   endsubroutine initialize

   subroutine parse_command_line_interface()
   !< Parse Command Line Interface (CLI).
   type(command_line_interface)  :: cli                   !< Command line interface handler.
   character(99)                 :: riemann_solver_scheme !< Riemann Problem solver scheme.
   integer(I4P)                  :: error                 !< Error handler.
   character(len=:), allocatable :: buffer                !< String buffer.

   call cli%init(description = 'FORESEER test: shock tube tester, 1D Euler equations', &
                 examples    = ["foreseer_test_shock_tube         ",                   &
                                "foreseer_test_shock_tube --tserie"])
   call cli%add(switch='--Ni', help='Number finite volumes used', required=.false., act='store', def='100')
   call cli%add(switch='--steps', help='Number time steps performed', required=.false., act='store', def='60')
   call cli%add(switch='--t-max', help='Maximum integration time', required=.false., act='store', def='0.')
   call cli%add(switch='--riemann', help='Riemann Problem solver', required=.false., act='store', def='all', &
                choices='all,exact,hllc,llf,pvl,roe')
   call cli%add(switch='--s-scheme', help='Space intergation scheme', required=.false., act='store', def='weno-char-1',           &
     choices='weno-char-1,weno-char-3,weno-char-5,weno-char-7,weno-char-9,weno-char-11,weno-char-13,weno-char-15,weno-char-17,'// &
             'weno-cons-1,weno-cons-3,weno-cons-5,weno-cons-7,weno-cons-9,weno-cons-11,weno-cons-13,weno-cons-15,weno-cons-17,'// &
             'weno-prim-1,weno-prim-3,weno-prim-5,weno-prim-7,weno-prim-9,weno-prim-11,weno-prim-13,weno-prim-15,weno-prim-17')
   call cli%add(switch='--t-scheme', help='Time intergation scheme', required=.false., act='store', def='tvd-rk-1', &
                choices='tvd-rk-1,tvd-rk-2,tvd-rk-3,tvd-rk-5')
   call cli%add(switch='--cfl', help='CFL value', required=.false., act='store', def='0.7')
   call cli%add(switch='--tserie', switch_ab='-t', help='Save time-serie-result', required=.false., act='store_true', def='.false.')
   call cli%add(switch='--verbose', help='Verbose output', required=.false., act='store_true', def='.false.')
   call cli%parse(error=error)
   call cli%get(switch='--Ni',       val=Ni,                    error=error) ; if (error/=0) stop
   call cli%get(switch='--steps',    val=steps_max,             error=error) ; if (error/=0) stop
   call cli%get(switch='--t-max',    val=t_max,                 error=error) ; if (error/=0) stop
   call cli%get(switch='--riemann',  val=riemann_solver_scheme, error=error) ; if (error/=0) stop
   call cli%get(switch='--s-scheme', val=s_scheme,              error=error) ; if (error/=0) stop
   call cli%get(switch='--t-scheme', val=t_scheme,              error=error) ; if (error/=0) stop
   call cli%get(switch='--cfl',      val=CFL,                   error=error) ; if (error/=0) stop
   call cli%get(switch='--tserie',   val=time_serie,            error=error) ; if (error/=0) stop
   call cli%get(switch='--verbose',  val=verbose,               error=error) ; if (error/=0) stop

   if (t_max > 0._R8P) steps_max = 0

   buffer = trim(adjustl(s_scheme))
   select case(buffer(6:9))
   case('char')
      weno_variables = 'characteristic'
   case('cons')
      weno_variables = 'conservative'
   case('prim')
      weno_variables = 'primitive'
   endselect
   weno_order = cton(buffer(11:), knd=1_I4P)

   select case(trim(adjustl(t_scheme)))
   case('tvd-rk-1')
      rk_stages_number = 1
   case('tvd-rk-2')
      rk_stages_number = 2
   case('tvd-rk-3')
      rk_stages_number = 3
   case('tvd-rk-5')
      rk_stages_number = 5
   endselect

   if (trim(adjustl(riemann_solver_scheme))=='all') then
      riemann_solver_schemes = ['exact', 'hllc ', 'llf  ', 'pvl  ', 'roe  ']
   else
      riemann_solver_schemes = [trim(adjustl(riemann_solver_scheme))]
   endif
   endsubroutine parse_command_line_interface

   subroutine save_time_serie(filename, finish, t)
   !< Save time-serie results.
   character(*), intent(in), optional :: filename  !< Output filename.
   logical,      intent(in), optional :: finish    !< Flag for triggering the file closing.
   real(R8P),    intent(in)           :: t         !< Current integration time.
   integer(I4P), save                 :: tsfile    !< File unit for saving time serie results.
   type(primitive_compressible)       :: primitive !< Primitive variables.
   integer(I4P)                       :: i         !< Counter.

   if (time_serie) then
      if (present(filename)) then
         open(newunit=tsfile, file=filename)
      endif
      write(tsfile, '(A)')'VARIABLES = "x" "rho" "u" "p"'
      write(tsfile, '(A)')'ZONE T="'//str(n=t)//'"'
      do i=1, Ni
         primitive = conservative_to_primitive_compressible(conservative=domain%U(i), eos=domain%eos)
         write(tsfile, '(4'//'('//FR8P//',1X))')x(i), primitive%density, primitive%velocity%x, primitive%pressure
      enddo
      if (present(finish)) then
         if (finish) close(tsfile)
      endif
   endif
   endsubroutine save_time_serie
endprogram foreseer_test_shock_tube