vortex3d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: vortex2d.f90                                                      #
!#                                                                           #
!# Copyright (C) 2006-2014                                                   #
!# 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 vortex3d
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 30.0     ! simulation stop time
  REAL, PARAMETER    :: gamma   = 1.4      ! ratio of specific heats
  REAL, PARAMETER    :: csiso   = &
!                                   0.0      ! non-isothermal simulation
                                  1.127    ! isothermal simulation
                                           !   with CSISO as sound speed
  ! initial condition (dimensionless units)
  REAL, PARAMETER    :: rhoinf  = 1.       ! ambient density
  REAL, PARAMETER    :: pinf    = 1.       ! ambient pressure
  REAL, PARAMETER    :: vstr    = 5.0      ! nondimensional vortex strength
  REAL, PARAMETER    :: uinf    = 0.0      ! cartesian components of constant
  REAL, PARAMETER    :: vinf    = 0.0      !   global velocity field
  REAL, PARAMETER    :: winf    = 0.0
  REAL, PARAMETER    :: x0      = 0.0      ! vortex position (cart. coords.)
  REAL, PARAMETER    :: y0      = 0.0
  REAL, PARAMETER    :: z0      = 0.0
  REAL, PARAMETER    :: r0      = 1.0      ! size of vortex
  REAL, PARAMETER    :: omega   = 0.0      ! angular speed of rotational frame
                                           ! around [X0,Y0]
  ! mesh settings
  INTEGER, PARAMETER :: mgeo = cylindrical
  INTEGER, PARAMETER :: xres = 40          ! x-resolution
  INTEGER, PARAMETER :: yres = 40          ! y-resolution
  INTEGER, PARAMETER :: zres = 1           ! z-resolution
  REAL, PARAMETER    :: rmin = 1.0e-2      ! inner radius for polar grids
  REAL, PARAMETER    :: rmax = 5.0         ! outer radius
  REAL, PARAMETER    :: zmin = -1.0        ! extent in z-direction
  REAL, PARAMETER    :: zmax = 1.0
  REAL, PARAMETER    :: gpar = 1.0         ! geometry scaling parameter     !
  ! physics settings
!!$  LOGICAL, PARAMETER :: WITH_IAR = .TRUE.  ! use EULER2D_IAMROT
  LOGICAL, PARAMETER :: with_iar = .false.

  ! output parameters
  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 = 'vortex3d'
  !--------------------------------------------------------------------------!
  CLASS(fosite),ALLOCATABLE :: sim
  CLASS(marray_compound), POINTER :: pvar_init
  REAL :: sigma
  !--------------------------------------------------------------------------!

  tap_plan(1)

  ALLOCATE(sim)
  CALL sim%InitFosite()
  CALL makeconfig(sim, sim%config)
  CALL sim%Setup()
  CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc)
  CALL sim%Physics%new_statevector(pvar_init,primitive)
  pvar_init = sim%Timedisc%pvar
  CALL sim%Run()
  ! compare with initial condition
  sigma = sqrt(sum((sim%Timedisc%pvar%data1d(:)-pvar_init%data1d(:))**2)/SIZE(pvar_init%data1d))
  CALL sim%Finalize()
  CALL pvar_init%Destroy()
  DEALLOCATE(sim,pvar_init)
  tap_check_small(sigma,3.8e-2,"vortex3d")
  tap_done

CONTAINS

  SUBROUTINE makeconfig(Sim, config)
    USE functions, ONLY : asinh
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(fosite)            :: Sim
    TYPE(Dict_TYP),POINTER  :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER                 :: bc(6)
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
                               timedisc, fluxes, sources, rotframe
    REAL                    :: x1,x2,y1,y2,z1,z2
    !------------------------------------------------------------------------!
    INTENT(INOUT)           :: sim
    !------------------------------------------------------------------------!
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(cartesian)
       x1 = -rmax
       x2 =  rmax
       y1 = -rmax
       y2 =  rmax
       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)
       x1 = rmin
       x2 = rmax
       y1 =  0.0
       y2 =  2.0*pi
       bc(west)   = no_gradients
       bc(east)   = no_gradients
       bc(south)  = periodic
       bc(north)  = periodic
       bc(bottom) = no_gradients
       bc(top)    = no_gradients
    CASE DEFAULT
       CALL sim%Error("vortex3d::MakeConfig","geometry currently not supported")
    END SELECT
    IF (zres.GT.1) THEN
      z1 =  zmin
      z2 =  zmax
    ELSE
      z1 = 0
      z2 = 0
    END IF

    !mesh settings
    mesh => dict( &
              "meshtype" / midpoint, &
              "geometry" / mgeo,     &
              "omega"    / omega,    &
              "inum"     / xres,     &
              "jnum"     / yres,     &
              "knum"     / zres,     &
              "xmin"     / x1,       &
              "xmax"     / x2,       &
              "ymin"     / y1,       &
              "ymax"     / y2,       &
              "zmin"     / z1,       &
              "zmax"     / z2,       &
              "gparam"   / gpar      )

    ! boundary conditions
    boundary => dict( &
                "western"   / bc(west),   &
                "eastern"   / bc(east),   &
                "southern"  / bc(south),  &
                "northern"  / bc(north),  &
                "bottomer"  / bc(bottom), &
                "topper"    / bc(top)     )

    ! physics settings
    IF (csiso.GT.tiny(csiso)) THEN
       physics => dict("problem" / euler_isotherm, &
                 "cs"      / csiso)                      ! isothermal sound speed  !
    ELSE
!       IF (WITH_IAR) THEN
!          ! REMARK: the optimal softening parameter depends on mesh geometry, limiter and
!          ! possibly other settings; modify this starting with the default of 1.0, if the
!          ! results show odd behaviour near the center of rotation; larger values increase
!          ! softening; 0.5 give reasonable results for PP limiter on cartesian mesh
!          physics => Dict("problem" / EULER2D_IAMT,        &
!                     "centrot_x"    / X0,"centrot_y" / Y0, & ! center of rotation      !
!                     "softening"    / 0.5,                 & ! softening parameter     !
!                     "gamma"        / GAMMA                ) ! ratio of specific heats !
!       ELSE
          physics => dict("problem"   / euler, &
                    "gamma"     / gamma)
!       END IF
    END IF

    ! flux calculation and reconstruction method
    fluxes => dict( &
              "fluxtype"  / kt,        &
!              "order"     / CONSTANT, &
              "order"     / linear,    &
              "variables" / primitive, &
              "limiter"   / vanleer,   &                     ! PP limiter gives better results than
              "theta"     / 1.0e-20,   &                     ! should be < EPS for PP limiter
              "output/slopes" / 0)

    ! activate inertial forces due to rotating frame if OMEGA > 0
    NULLIFY(sources)
    IF (omega.GT.tiny(omega)) THEN
       rotframe => dict("stype" / rotating_frame, &
               "x"     / x0, &
               "y"     / y0)
       sources => dict("rotframe" / rotframe)
    END IF

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

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

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

    ! add sources terms
    IF (ASSOCIATED(sources)) &
        CALL setattr(config, "sources", sources)

  END SUBROUTINE makeconfig


  SUBROUTINE initdata(Mesh,Physics,Timedisc)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(physics_base),INTENT(IN)    :: Physics
    CLASS(mesh_base),INTENT(IN)       :: Mesh
    CLASS(timedisc_base),INTENT(INOUT):: Timedisc
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: n
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
                      :: radius,dist_axis,domega
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) &
                      :: posvec,ephi
    REAL              :: csinf
    !------------------------------------------------------------------------!
    IF (abs(x0).LE.tiny(x0).AND.abs(y0).LE.tiny(y0).AND.abs(z0).LE.tiny(z0)) THEN
       ! no shift of rotational axis
       radius(:,:,:) = mesh%radius%bcenter(:,:,:)
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
    ELSE
       ! axis is assumed to be parallel to the cartesian z-axis
       ! but shifted in the plane perpendicular to it
       ! compute curvilinear components of shift vector
       posvec(:,:,:,1) = x0
       posvec(:,:,:,2) = y0
       posvec(:,:,:,3) = 0.0
       CALL mesh%geometry%Convert2Curvilinear(mesh%bcenter,posvec,posvec)
       ! subtract the result from the position vector:
       ! this gives you the curvilinear components of all vectors pointing
       ! from the point mass to the bary center of any cell on the mesh
       posvec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - posvec(:,:,:,:)
       ! compute its absolute value
       radius(:,:,:) = sqrt(posvec(:,:,:,1)**2+posvec(:,:,:,2)**2+posvec(:,:,:,3)**2)
    END IF
    ! distance to axis of rotation
    dist_axis(:,:,:) = sqrt(posvec(:,:,:,1)**2+posvec(:,:,:,2)**2)

    ! curvilinear components of azimuthal unit vector
    ! (maybe with respect to shifted origin)
    ! from ephi = ez x er = ez x posvec/radius = ez x (rxi*exi + reta*eeta)/r
    !             = rxi/r*(ez x exi) + reta/r*(ez x eeta) = rxi/r*eeta - reta/r*exi
    ! because (ez,exi,eeta) is right handed orthonormal set of basis vectors
    ephi(:,:,:,1) = -posvec(:,:,:,2)/radius(:,:,:)
    ephi(:,:,:,2) =  posvec(:,:,:,1)/radius(:,:,:)
    ephi(:,:,:,3) = 0.0

    ! initial condition
    ! angular velocity
    domega(:,:,:) = 0.5*vstr/pi*exp(0.5*(1.-(dist_axis(:,:,:)/r0)**2))
    SELECT TYPE(pvar => timedisc%pvar)
    TYPE IS(statevector_eulerisotherm) ! isothermal HD
      ! density
      pvar%density%data3d(:,:,:)  = rhoinf * exp(-0.5*(r0*domega(:,:,:)/csiso)**2)
      ! velocities
      DO n=1,physics%VDIM
        pvar%velocity%data4d(:,:,:,n) = domega(:,:,:)*dist_axis(:,:,:)*ephi(:,:,:,n)
      END DO
    TYPE IS(statevector_euler) ! non-isothermal HD
      csinf = sqrt(gamma*pinf/rhoinf) ! sound speed at infinity (isentropic vortex)
      ! density
      ! ATTENTION: there's a factor of 1/PI missing in the density
      ! formula  eq. (3.3) in [1]
      pvar%density%data3d(:,:,:)  = rhoinf * (1.0 - &
              0.5*(gamma-1.0)*(r0*domega/csinf)**2  )**(1./(gamma-1.))
      ! pressure
      pvar%pressure%data1d(:) = pinf * (pvar%density%data1d(:)/rhoinf)**gamma
      ! velocities
      DO n=1,physics%VDIM
        pvar%velocity%data4d(:,:,:,n) = domega(:,:,:)*dist_axis(:,:,:)*ephi(:,:,:,n)
      END DO
    END SELECT

    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)

!     ! compute curvilinear components of constant background velocity field
!     ! and add to the vortex velocity field
!     IF (ABS(UINF).GT.TINY(UINF).OR.ABS(VINF).GT.TINY(VINF).AND.ABS(WINF).LE.TINY(WINF)) THEN
!        v0(:,:,:,1) = UINF
!        v0(:,:,:,2) = VINF
!        v0(:,:,:,3) = WINF
!        CALL Mesh%geometry%Convert2Curvilinear(Mesh%bcenter,v0,v0)
!        Timedisc%pvar%data4d(:,:,:,Physics%XVELOCITY:Physics%ZVELOCITY) = &
!            Timedisc%pvar%data4d(:,:,:,Physics%XVELOCITY:Physics%ZVELOCITY) + v0(:,:,:,1:3)
!     END IF
! 
! 
    CALL mesh%Info(" DATA-----> initial condition: 3D vortex")
  END SUBROUTINE initdata

END PROGRAM vortex3d