poiseuille.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: poiseuille.f90                                                    #
!#                                                                           #
!# Copyright (C) 2006-2019                                                   #
!# Bjoern Sperling  <sperling@astrophysik.uni-kiel.de>                       #
!# 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 poiseuille
  USE fosite_mod
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER :: tsim  = 30.0                ! simulation time            !
  REAL, PARAMETER :: re    = 4.375              ! Reynolds Number            !
  REAL, PARAMETER :: pin   = 24.0               ! inflow pressure            !
  REAL, PARAMETER :: pout  = 23.0               ! outflow pressure           !
  REAL, PARAMETER :: rho0  = 1.0                ! initial density            !
  ! mesh settings
  INTEGER, PARAMETER :: mgeo = cylindrical
!  INTEGER, PARAMETER :: MGEO = CARTESIAN
  INTEGER, PARAMETER :: xres = 10               ! radial resolution          !
  INTEGER, PARAMETER :: yres = 10               ! azimuthal resolution       !
  INTEGER, PARAMETER :: zres = 10               ! resolution along the tube  !
  REAL, PARAMETER    :: ltube = 10.0            ! length of the tube         !
  REAL, PARAMETER    :: rtube = 1.0             ! radius of the tube         !
  !--------------------------------------------------------------------------!
  ! output file parameter
  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 = 'poiseuille'
  !--------------------------------------------------------------------------!
  CLASS(fosite),ALLOCATABLE   :: sim
  !--------------------------------------------------------------------------!
  ALLOCATE(sim)
  CALL sim%InitFosite()

  CALL makeconfig(sim, sim%config)

!  CALL PrintDict(config)

  CALL sim%Setup()

  ! set initial condition
  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)

  CALL sim%Run()
  CALL sim%Finalize()
  DEALLOCATE(sim)

CONTAINS

  SUBROUTINE makeconfig(Sim, config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Fosite)  :: Sim
    TYPE(Dict_TYP),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: bc(6)
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, sources, &
                               vis, timedisc, fluxes
    REAL              :: x1,x2,y1,y2,z1,z2
    REAL              :: dvis, bvis
    !------------------------------------------------------------------------!
    INTENT(INOUT)     :: sim
    !------------------------------------------------------------------------!
    ! geometry dependent setttings
    SELECT CASE(mgeo)
    CASE(cylindrical)
      x1 = 0.0
      x2 = rtube
      y1 = 0.0
      y2 = 2*pi
      z1 = 0.0
      z2 = ltube
      bc(west)  = axis
      bc(east)  = noslip
      bc(south) = periodic
      bc(north) = periodic
      bc(bottom)= fixed
      bc(top)   = fixed
    CASE(cartesian)
      x1 = 0.0
      x2 = rtube
      y1 = 0.0
      y2 = rtube
      z1 = 0.0
      z2 = ltube
      bc(west)  = noslip
      bc(east)  = noslip
      bc(south) = noslip
      bc(north) = noslip
      bc(bottom)= fixed
      bc(top)   = fixed
    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)

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

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

    ! viscosity source term
    ! dynamic viscosity
    dvis = sqrt(0.25 * (pin-pout) * rtube**3 * rho0 / ltube / re)
    ! bulk viscosity
    bvis = -2./3. * dvis
    vis => dict("stype"     / viscosity, &
          "vismodel"  / molecular, &
          "dynconst"  / dvis, &
          "bulkconst" / bvis)

    sources => dict("vis" / vis)

    ! time discretization settings
    timedisc => dict("method"   / modified_euler, &
               "order"    / 3, &
               "cfl"      / 0.4, &
               "stoptime" / tsim, &
               "dtlimit"  / 1.0e-8, &
               "maxiter"  / 1000000)

    ! initialize data input/output
    datafile => dict( &
               "fileformat" / xdmf, &
               "filename"   / (trim(odir) // trim(ofname)), &
               "count"      / onum)

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

  SUBROUTINE initdata(Mesh,Physics,Timedisc)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Physics_base) :: Physics
    CLASS(Mesh_base)    :: Mesh
    CLASS(Timedisc_base):: Timedisc
    !------------------------------------------------------------------------!
    INTENT(IN)        :: mesh,physics
    INTENT(INOUT)     :: timedisc
    !------------------------------------------------------------------------!
    ! initial condition
    timedisc%pvar%data4d(:,:,:,physics%DENSITY)   = rho0
    timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.
    timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.
    timedisc%pvar%data4d(:,:,:,physics%ZVELOCITY) = 0.
    timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = pin + (pout-pin)/ltube &
         * mesh%bccart(:,:,:,3)

    ! fixed boundary conditions
    ! inflow
    SELECT TYPE(bbottom => timedisc%Boundary%boundary(bottom)%p)
    CLASS IS (boundary_fixed)
       bbottom%data(:,:,:,physics%DENSITY)    = rho0
       bbottom%data(:,:,:,physics%XVELOCITY)  = 0.0
       bbottom%data(:,:,:,physics%YVELOCITY)  = 0.0
       bbottom%data(:,:,:,physics%ZVELOCITY)  = 0.0
       bbottom%data(:,:,:,physics%PRESSURE)   = pin
       ! imposed density, pressure and tangential velocities;
       ! extrapolated normal velocity
       bbottom%fixed(:,:,physics%DENSITY)   = .true.
       bbottom%fixed(:,:,physics%XVELOCITY) = .true.
       bbottom%fixed(:,:,physics%YVELOCITY) = .true.
       bbottom%fixed(:,:,physics%ZVELOCITY) = .false.
       bbottom%fixed(:,:,physics%PRESSURE)  = .true.
    END SELECT
    ! outflow
    SELECT TYPE(btop => timedisc%Boundary%boundary(top)%p)
    CLASS IS (boundary_fixed)
       btop%data(:,:,:,physics%DENSITY)    = rho0
       btop%data(:,:,:,physics%XVELOCITY)  = 0.0
       btop%data(:,:,:,physics%YVELOCITY)  = 0.0
       btop%data(:,:,:,physics%ZVELOCITY)  = 0.0
       btop%data(:,:,:,physics%PRESSURE)   = pout
       ! imposed pressure;
       ! extrapolated density, tangential and normal velocities
       btop%fixed(:,:,physics%DENSITY)   = .false.
       btop%fixed(:,:,physics%XVELOCITY) = .false.
       btop%fixed(:,:,physics%YVELOCITY) = .false.
       btop%fixed(:,:,physics%ZVELOCITY) = .false.
       btop%fixed(:,:,physics%PRESSURE)  = .true.
    END SELECT

    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
    CALL mesh%Info(" DATA-----> initial condition: " // &
         "tube with pressure gradient")

  END SUBROUTINE initdata

END PROGRAM poiseuille