DCS.f90

Go to the documentation of this file.

program dcs
!-----------------------------------------------------------------------------------------------------------------------------------
use, intrinsic:: iso_fortran_env, only: stdout => output_unit, stderr => error_unit ! Standard output/error logical units.
USE ir_precision                                                                    ! Integers and reals precision definition.
USE data_type_cavity                                                                ! Definition of Type_Cavity.
USE data_type_conservative                                                          ! Definition of Type_Conservative.
!-----------------------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------------------
implicit none
type(type_cavity):: cavity
real(R8P)::         rtol = 0.0001_r8p
integer(I4P)::      nout = 1000_i4p
character(3)::      oform = "tec"
!-----------------------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------------------
call dcs_initialize ! initializing the simulation
call dcs_simulate   ! performing the simulation
call dcs_finalize   ! finalizing the simulation
stop
!-----------------------------------------------------------------------------------------------------------------------------------
contains
  integer function get_unit(Free_Unit)
  !---------------------------------------------------------------------------------------------------------------------------------
  implicit none
  integer, intent(OUT), optional:: free_unit
  integer::                        n1
  integer::                        ios
  logical::                        lopen
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  get_unit = -1
  n1=1
  do
    if ((n1/=stdout).AND.(n1/=stderr)) then
      inquire (unit=n1,opened=lopen,iostat=ios)
      if (ios==0) then
        if (.NOT.lopen) then
          get_unit = n1 ; if (present(free_unit)) free_unit = get_unit
          return
        endif
      endif
    endif
    n1=n1+1
  enddo
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endfunction get_unit

  subroutine print_usage()
  !---------------------------------------------------------------------------------------------------------------------------------
  ! Subroutine for printing the correct use of the program.
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  implicit none
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  write(stdout,'(A)')' DCS'
  write(stdout,'(A)')' Driven-Cavity Simulator code'
  write(stdout,'(A)')' Usage:'
  write(stdout,'(A)')' ./DCS -N number_of_cells [optionals parameters]'
  write(stdout,*)
  write(stdout,'(A)')' Optional parameters:'
  write(stdout,'(A)')' -Re Reynolds_number => default 500'
  write(stdout,'(A)')' -beta over_relaxation_parameter => default 0.6'
  write(stdout,'(A)')' -rtol residual_tolerance => default 10^-4'
  write(stdout,'(A)')' -nout standard_output_update_frequency => default 1000 iterations'
  write(stdout,'(A)')' -oform output_file_format => "tec" for Tecplot output or "gnu" for gnuplot one, default "tec"'
  write(stdout,*)
  write(stdout,'(A)')' Example: '
  write(stdout,'(A)')' ./DCS -Re 9.d2 -beta 0.5d0 -N 256 -nout 10000 -rtol 1.d-6'
  write(stdout,*)
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine print_usage

  subroutine dcs_initialize()
  !---------------------------------------------------------------------------------------------------------------------------------
  implicit none
  real(R8P)::      re
  real(R8P)::      beta
  integer(I4P)::   n
  integer(I4P)::   nca = 0
  integer(I4P)::   c
  character(6)::   ca_switch
  character(100):: sbuf
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  re = cavity%Re
  beta = cavity%beta
  n = cavity%mesh%N
  ! parsing command line arguments
  nca = command_argument_count()
  if (nca==0) then
    write(stderr,'(A)')' Error: no argument has been passed to command line'
    call print_usage
    stop
  else
    ! processing switch
    c = 0
    do while (c<nca)
      c = c + 1
      call get_command_argument(c, ca_switch)
      select case(adjustl(trim(ca_switch)))
      case('-Re')
        call get_command_argument(c+1,sbuf) ; c = c + 1
        read(sbuf,*)re
      case('-beta')
        call get_command_argument(c+1,sbuf) ; c = c + 1
        read(sbuf,*)beta
      case('-N')
        call get_command_argument(c+1,sbuf) ; c = c + 1
        read(sbuf,*)n
      case('-rtol')
        call get_command_argument(c+1,sbuf) ; c = c + 1
        read(sbuf,*)rtol
      case('-nout')
        call get_command_argument(c+1,sbuf) ; c = c + 1
        read(sbuf,*)nout
      case('-oform')
        call get_command_argument(c+1,oform) ; c = c + 1
      case default
        write(stderr,'(A)') ' Error: switch "'//adjustl(trim(ca_switch))//'" is unknown'
        call print_usage
        stop
      endselect
    enddo
  endif
  if (n==0_i4p) then
    write(stderr,'(A)')' Error: the number of cells N must be passed as command line argument'
    call print_usage
    stop
  endif
  write(stdout,'(A)')' Some information about the precision of runnig machine'
  call ir_print()
  write(stdout,*)
  write(stdout,'(A)')' Simulation parameters'
  write(stdout,'(A)')' Re='//trim(str(n=re))
  write(stdout,'(A)')' beta='//trim(str(n=beta))
  write(stdout,'(A)')' N='//trim(str(n=n))
  write(stdout,'(A)')' rtol='//trim(str(n=rtol))
  write(stdout,'(A)')' nout='//trim(str(n=nout))
  write(stdout,'(A)')' oform='//trim(oform)
  write(stdout,*)
  ! initializing cavity data
  call cavity%initialize(re=re,beta=beta,n=n)
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine dcs_initialize

  subroutine dcs_simulate()
  !---------------------------------------------------------------------------------------------------------------------------------
  implicit none
  type(type_conservative):: residual
  integer(I4P)::            n = 0_i4p
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! successive over-relaxation (SOR) loop
  call residual%set(s=maxr8p,v=maxr8p)
  sor: do while(((residual%s>rtol).AND.(residual%v>rtol)))
    ! updating time step counter
    n = n + 1_i4p
    ! computing new stream function
    call computestream(cons=cavity%cons,beta=cavity%beta,dh=cavity%mesh%dh)
    ! computing new vorticity
    call computevorticity(cons=cavity%cons,re=cavity%Re,beta=cavity%beta,dh=cavity%mesh%dh)
    ! computing residuals
    residual = computeresidual(cons=cavity%cons,re=cavity%Re,dh=cavity%mesh%dh)
    if (mod(n,nout)==0) then
      write(stdout,'(A)')' iteration n='//trim(str(n=n))//&
                         ' residuals: s='//trim(str(n=residual%s))//' v='//trim(str(n=residual%v))
    endif
  enddo sor
  write(stdout,'(A)')' Convergence achieved n='//trim(str(n=n))
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine dcs_simulate

  subroutine dcs_finalize()
  !---------------------------------------------------------------------------------------------------------------------------------
  implicit none
  integer(I4P):: u
  integer(I4P):: i,j
  !---------------------------------------------------------------------------------------------------------------------------------

  !---------------------------------------------------------------------------------------------------------------------------------
  ! writing solution
  select case(trim(oform))
  case('tec')
    ! tecplot format
    open(unit=get_unit(u),file='DCS_out.tec')
    write(u,'(A)')'VARIABLES="x" "y" "s" "v"'
    write(u,'(A)')'ZONE T="Driven Cavity" I='//trim(str(n=cavity%mesh%N+1))//' J='//trim(str(n=cavity%mesh%N+1))
    do j=0,cavity%mesh%N
      do i=0,cavity%mesh%N
        write(u,'(A)')trim(str(n=cavity%mesh%x(i)))//' '//trim(str(n=cavity%mesh%y(j)))//' '//&
                      trim(str(n=cavity%cons(i,j)%s))//' '//trim(str(n=cavity%cons(i,j)%v))
      enddo
    enddo
    close(u)
  case('gnu')
    ! gnuplot format
    open(unit=get_unit(u),file='DCS_out.gnu')
    do j=0,cavity%mesh%N
      do i=0,cavity%mesh%N
        write(u,'(A)')trim(str(n=cavity%mesh%x(i)))//' '//trim(str(n=cavity%mesh%y(j)))//' '//&
                      trim(str(n=cavity%cons(i,j)%s))//' '//trim(str(n=cavity%cons(i,j)%v))
      enddo
      write(u,*)
    enddo
    close(u)
  endselect
  ! finalizing cavity data
  call cavity%finalize
  return
  !---------------------------------------------------------------------------------------------------------------------------------
  endsubroutine dcs_finalize
endprogram dcs