keplerianvortex.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: keplerianvortex.f90                                               #
!#                                                                           #
!# Copyright (C) 2014-2018                                                   #
!# Manuel Jung <mjung@astrophysik.uni-kiel.de>                               #
!# Jannes Klee <jklee@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 keplerianvortex
  USE fosite_mod
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER :: gn          = 1.0
  REAL, PARAMETER :: mcentral    = 1.0
  REAL, PARAMETER :: sigma0      = 1.0
  REAL, PARAMETER :: m           = 10.
  REAL, PARAMETER :: omega       = 0.

  REAL, PARAMETER :: p0          = 2.0*pi/1.0
  REAL, PARAMETER :: p           = 3
  REAL, PARAMETER :: t           = p*p0

  REAL, PARAMETER :: gamma       = 5./3.

  REAL, PARAMETER :: rmin        = 0.5
  REAL, PARAMETER :: rmax        = 1.5

  INTEGER, PARAMETER :: gravenergysource = 1
  INTEGER, PARAMETER :: xres     = 64
  INTEGER, PARAMETER :: yres     = 4*xres
  INTEGER, PARAMETER :: onum     = p * 10
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE          :: sim
  !--------------------------------------------------------------------------!

ALLOCATE(sim)
CALL sim%InitFosite()
CALL makeconfig(sim%config)
CALL sim%Setup()
CALL initdata(sim%Timedisc, sim%Mesh, sim%Physics, sim%Fluxes)
CALL sim%Run()
CALL sim%Finalize()
DEALLOCATE(sim)

CONTAINS

  SUBROUTINE makeconfig(config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(Dict_TYP),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(Dict_TYP),POINTER :: mesh,boundary,timedisc,datafile,&
                              sources,fluxes,pmass,physics,damping,vis,pot
    !------------------------------------------------------------------------!
    physics => dict( &
        "problem" / euler, &
        "gamma" / gamma, &
        "units" / geometrical)

    fluxes => dict( &
        "fluxtype" / kt, &
        "order" / linear, &
        "variables" / primitive, &
        "limiter" / vanleer)

    mesh => dict( &
        "meshtype" / midpoint, &
        "geometry" / cylindrical, &
        "omega" / omega, &
        "inum" / xres, &
        "jnum" / yres, &
        "knum" / 1, &
        "xmin" / rmin, &
        "xmax" / rmax, &
        "ymin" / (-pi), &
        "ymax" / (pi), &
        "zmin" / 0.0, &
        "zmax" / 0.0, &
        "decomposition" / (/ 1, -1/), &
        "output/volume" / 1,&
        "output/bh" / 1,&
        "output/rotation" / 1,&
        "output/dl" / 1)

    boundary => dict( &
        "western" / custom, &
        "eastern" / custom, &
        "southern" / periodic, &
        "northern" / periodic)

    pmass => dict( &
        "gtype" / pointmass, &
        "potential" / newton, &
        "mass" / mcentral, &
        "cvis" / 0.9, &
        "outbound" / 0)!, &
!        "output/accel" / 1 )

!    damping => Dict( &
!        "stype" / WAVE_DAMPING, &
!        "zone0/rmin" / RMIN, &
!        "zone0/rmax" / (0.5), &
!        "zone0/tau"  / (0.1*2.0*PI / SQRT(GN*MCENTRAL/(0.4)**3)), &
!        "zone1/rmin" / (1.9), &
!        "zone1/rmax" / RMAX, &
!        "zone1/tau"  / (0.1*2.0*PI / SQRT(GN*MCENTRAL/(2.0)**3)))

    sources => dict( &
!        "damping" / damping, &
        "grav/stype" / gravity, &
!        "grav/energy" / GravEnergySource, &
        "grav/output/accel" / 1, &
        "grav/pmass"   / pmass)

    timedisc => dict( &
        "method" / ssprk, &
!         "fargo" / 1, &
        "tol_rel" / 1.0e-3, &
        "tol_abs" / (/ 0., 1.e-3, 1.e-3, 0. /), &
        "cfl" / 0.4, &
        "stoptime" / t, &
        "dtlimit" / 1.0e-10, &
        "maxiter" / 2000000000, &
        "output/xmomentum" / 1, &
        "output/ymomentum" / 1, &
        "output/energy" / 1, &
!        "rhstype" / 1, &
        "output/geometrical_sources" / 1, &
        "output/external_sources" / 1, &
        !"output/bflux" / 1, &
        "output/fluxes" / 1, &
        "output/rhs" / 1)


    datafile => dict( &
        "fileformat" / xdmf, &
        "filename" / "keplerianvortex", &
        "count" / onum)

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

  END SUBROUTINE makeconfig


  SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(INOUT) :: Physics
    CLASS(fluxes_base),   INTENT(INOUT) :: Fluxes
    !------------------------------------------------------------------------!
    ! Local variable declaration
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
                      :: dvr, dvphi, x, y, phi, r
    REAL              :: h, kappa
    INTEGER           :: dir,j,ig,i,k
    CLASS(sources_base), POINTER :: sp
    CLASS(sources_gravity), POINTER :: gp
    !------------------------------------------------------------------------!
    x = mesh%bccart(:,:,:,1)
    y = mesh%bccart(:,:,:,2)

    IF(gravenergysource.EQ.0) THEN
      CALL setpotential_euler2dia(physics,mesh,-gn * mcentral / mesh%radius%bcenter)
      CALL setpotential_euler2dia(physics,mesh,-gn * mcentral / mesh%radius%faces)
    END IF

    ! get gravitational acceleration
    sp => sim%Sources
    DO
      IF (ASSOCIATED(sp).EQV..false.) RETURN
      SELECT TYPE(sp)
      CLASS IS(sources_gravity)
        gp => sp
        EXIT
      END SELECT
      sp => sp%next
    END DO

    SELECT TYPE (pvar => timedisc%pvar)
    CLASS IS(statevector_euler)
      SELECT TYPE (cvar => timedisc%cvar)
      CLASS IS(statevector_euler)
        cvar%density%data3d(:,:,:) = sigma0
        cvar%momentum%data4d(:,:,:,1:physics%VDIM) = 0.

        CALL physics%Convert2Primitive(cvar,pvar)

        pvar%pressure%data3d(:,:,:) = 1./(gamma*m**2)

        pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
          pvar%velocity%data4d(:,:,:,1:physics%VDIM) &
          + timedisc%GetCentrifugalVelocity(mesh,physics,sim%Fluxes,sim%Sources,(/0.,0.,1./))

      END SELECT
    END SELECT

    ! boundary conditions
    ! custom boundary conditions at western boundary if requested
    SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
    CLASS IS (boundary_custom)
      CALL bwest%SetCustomBoundaries(mesh,physics, &
        (/custom_reflect,custom_reflneg,custom_kepler,custom_reflect/))
    END SELECT
    SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
    CLASS IS (boundary_custom)
      CALL beast%SetCustomBoundaries(mesh,physics, &
        (/custom_reflect,custom_reflneg,custom_kepler,custom_reflect/))
    END SELECT

    h = 1./(2.*m)
    kappa = -1.
    phi = mesh%center(:,:,:,2)
    x = x - 1./sqrt(2.)
    y = y - 1./sqrt(2.)
    r = sqrt(x**2 + y**2)

    dvr   = kappa * exp(-r**2/h**2) * (-y*cos(phi) + x*sin(phi))
    dvphi = kappa * exp(-r**2/h**2) * ( y*sin(phi) + x*cos(phi))

    SELECT TYPE (pvar => timedisc%pvar)
    CLASS IS(statevector_euler)
      SELECT TYPE (cvar => timedisc%cvar)
      CLASS IS(statevector_euler)
        pvar%velocity%data4d(:,:,:,1) = pvar%velocity%data4d(:,:,:,1) + dvr
        pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) + dvphi
      END SELECT
    END SELECT

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

  END SUBROUTINE initdata
END PROGRAM keplerianvortex