orbitingcylinders.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: orbitingcylinders.f90                                             #
!#                                                                           #
!# Copyright (C) 2013-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.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
!#define HAVE_FGSL
PROGRAM orbitingcylinders
  USE fosite_mod
#ifdef HAVE_FGSL
  USE fgsl
#endif
  USE functions, ONLY : ei,bessel_i0,bessel_i1,bessel_k0,bessel_k1,bessel_k0e
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! problem setup
  REAL, PARAMETER :: sigma       = 0.1
  REAL, DIMENSION(2), PARAMETER :: x1 = (/ 1.0, 0. /)
  REAL, DIMENSION(2), PARAMETER :: x2 = (/ 1.0, pi /)
  REAL, DIMENSION(2), PARAMETER :: x3 = (/ 0.9, 3./4.*pi /)
  REAL, DIMENSION(2), PARAMETER :: x4 = (/ 1.0, -pi/2. /)

  ! general constants
  REAL, PARAMETER :: gn          = 1.
  REAL, PARAMETER :: p0          = 1.
  REAL, PARAMETER :: p           = 10.
  REAL, PARAMETER :: t           = p*p0

  REAL, PARAMETER :: flaring     = 0.05
  REAL, PARAMETER :: rg          = 8.31
  REAL, PARAMETER :: mu          = 6.02e-04
  REAL, PARAMETER :: gamma       = 5./3.

  REAL, PARAMETER :: rmin        = 0.4
  REAL, PARAMETER :: rmax        = 2.0
  INTEGER, PARAMETER :: yres     = 128
  INTEGER, PARAMETER :: xres     = yres/4
  INTEGER, PARAMETER :: onum     = p * 1
  REAL, PARAMETER    :: omega    = 1.0
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  INTEGER :: green
  REAL, DIMENSION(:,:,:), POINTER :: numpot, anapot
  !--------------------------------------------------------------------------!

  ALLOCATE(sim)
  DO green=3,3


    CALL sim%InitFosite()
    CALL makeconfig(sim, sim%config)

    SELECT CASE(green)
    CASE(3)
      CALL setattr(sim%config, "/mesh/rmin", 0.2)
      CALL setattr(sim%config, "/mesh/rmax", 1.8)
      CALL setattr(sim%config, "/mesh/inum", 64)
      CALL setattr(sim%config, "/mesh/jnum", 192)
    END SELECT

    CALL sim%Setup()
    CALL initdata(sim%Timedisc,sim%Mesh,sim%Physics,sim%Sources)

    SELECT CASE(green)
    CASE(1,2)
      CALL sim%FirstStep()
      print *,"max local relative error: ",maxval(abs((numpot-anapot)/anapot))
      print *,"global relative error: ",sum(abs(numpot-anapot))/sum(abs(anapot))
    CASE(3)
      CALL sim%Run()
      CALL sim%Finalize()
    END SELECT
  END DO
  DEALLOCATE(sim)


CONTAINS

  SUBROUTINE makeconfig(Sim, config)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(fosite), INTENT(INOUT) :: Sim
    TYPE(Dict_TYP),POINTER :: config
    !------------------------------------------------------------------------!
    ! Local variable declaration
    TYPE(Dict_TYP),POINTER :: mesh,boundary,timedisc,datafile,&
                              fluxes,physics,grav,rotframe,sources
    !------------------------------------------------------------------------!
    CHARACTER(LEN=1)  :: str
    !------------------------------------------------------------------------!
    physics => dict( &
        "problem" / euler,&
        "mu" / mu, &
        "gamma" / gamma, &
        "units" / geometrical)

    fluxes => dict( &
        "fluxtype"  / kt, &!HLLC, &
        "order" / linear, &
        "variables" / primitive, &
        "limiter" / monocent, &
        "theta" / 1.2)

    mesh => dict( &
        "meshtype" / midpoint, &
        "geometry" / cylindrical, &
        "xmin" / rmin, &
        "xmax" / rmax, &
!        "geometry" / LOGPOLAR, &
!        "xmin" / LOG(RMIN/1.), &
!        "xmax" / LOG(RMAX/1.), &
        "gparam" / 1., &
        "omega" / omega, &
        "inum" / xres, &
        "jnum" / yres, &
        "knum" / 1, &
        "ymin" / (0.), &
        "ymax" / (2.*pi), &
        "zmin" / (0.), &
        "zmax" / (0.), &
        "decomposition" / (/ -1, 1/))!, &
!        "output/dl" / 1, &
!        "output/radius" / 1, &
!        "output/volume" / 1 )

    boundary => dict( &
        "western" / noslip, &
        "eastern" / noslip, &
        "southern" / periodic, &
        "northern" / periodic, &
        "bottomer" / reflecting, &
        "topper" / reflecting)

    grav => dict( &
        "stype" / gravity, &
        "output/accel" / 1, &
        "self/gtype" / spectral, &
        "self/green" / green, &
        "self/sigma" / sigma)
!        "self/output/potential" / 1)

    rotframe => dict( &
               "stype" / rotating_frame, &
               "x"     / 0.0, &
               "y"     / 0.0, &
               "z"     / 0.0)

    sources => dict( &
!              "rotframe"        / rotframe, &
              "grav"            / grav)

    timedisc => dict( &
        "method" / ssprk, &
        "rhstype" / 1, &
!        "dtmax" / 8., &
!        "fargo" / 1, &
!        "method" / RK_FEHLBERG, &
!        "order" / 5, &
        "tol_rel" / 1.0e-3, &
        "cfl" / 0.3, &
        "stoptime" / t, &
        "dtlimit" / 1.0e-10, &
        "maxiter" / 2000000000, &
        "output/solution" / 1)

    WRITE(str,'(I1)')green

    datafile => dict( &
        "fileformat" / vtk, &
        "filename" / ("orbitingcylinders_" // str), &
        "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,Sources)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(physics_base),  INTENT(INOUT) :: Physics
    CLASS(sources_base),  POINTER       :: Sources
    !------------------------------------------------------------------------!
    ! Local variable declaration
    INTEGER           :: i,j,k,dir,ig
    REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: r,r1,r2,r3
    CLASS(sources_base), POINTER :: sp
    CLASS(sources_gravity), POINTER :: gp
    !------------------------------------------------------------------------!
    r = mesh%radius%bcenter
    sp => 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


    anapot => timedisc%solution%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1)
    numpot => gp%glist%pot(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,1)

    SELECT CASE(green)
    CASE(1)
      r1 = rsq(mesh,mesh%radius%bcenter,x1)
      r2 = rsq(mesh,mesh%radius%bcenter,x3)
      r3 = rsq(mesh,mesh%radius%bcenter,x4)
      timedisc%pvar%data4d(:,:,:,physics%DENSITY) &
        = 2.0 / (2.*pi*sigma**2) * exp(-sqrt(r1)/sigma) &
        + 0.5 / (2.*pi*sigma**2) * exp(-sqrt(r2)/sigma) &
        + 1.0 / (2.*pi*sigma**2) * exp(-sqrt(r3)/sigma)

      r1 = sqrt(rsq(mesh,mesh%radius%faces(:,:,:,1),x1)) / (2.*sigma)
      r2 = sqrt(rsq(mesh,mesh%radius%faces(:,:,:,1),x3)) / (2.*sigma)
      r3 = sqrt(rsq(mesh,mesh%radius%faces(:,:,:,1),x4)) / (2.*sigma)
      DO k=mesh%KGMIN,mesh%KGMAX
        DO j=mesh%JGMIN,mesh%JGMAX
          DO i=mesh%IGMIN,mesh%IGMAX
            timedisc%solution%data4d(i,j,k,1) &
#ifdef HAVE_FGSL
           = -2.0 * r1(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r1(i,j,k))*fgsl_sf_bessel_kc1(r1(i,j,k))&
                    -fgsl_sf_bessel_ic1(r1(i,j,k))*fgsl_sf_bessel_kc0(r1(i,j,k))) &
             -0.5 * r2(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r2(i,j,k))*fgsl_sf_bessel_kc1(r2(i,j,k))&
                    -fgsl_sf_bessel_ic1(r2(i,j,k))*fgsl_sf_bessel_kc0(r2(i,j,k))) &
             -1.0 * r3(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r3(i,j,k))*fgsl_sf_bessel_kc1(r3(i,j,k))&
                    -fgsl_sf_bessel_ic1(r3(i,j,k))*fgsl_sf_bessel_kc0(r3(i,j,k)))
#else
           = -2.0 * r1(i,j,k)/sigma*(bessel_i0(r1(i,j,k))*bessel_k1(r1(i,j,k))&
                    -bessel_i1(r1(i,j,k))*bessel_k0(r1(i,j,k))) &
             -0.5 * r2(i,j,k)/sigma*(bessel_i0(r2(i,j,k))*bessel_k1(r2(i,j,k))&
                    -bessel_i1(r2(i,j,k))*bessel_k0(r2(i,j,k))) &
             -1.0 * r3(i,j,k)/sigma*(bessel_i0(r3(i,j,k))*bessel_k1(r3(i,j,k))&
                    -bessel_i1(r3(i,j,k))*bessel_k0(r3(i,j,k)))
#endif
          END DO
        END DO
      END DO

      timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = 1.0
      timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.0
      timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.0
    CASE(2)
      r1 = rsq(mesh,mesh%radius%bcenter,x1)
      r2 = rsq(mesh,mesh%radius%bcenter,x3)
      r3 = rsq(mesh,mesh%radius%bcenter,x4)
      timedisc%pvar%data4d(:,:,:,physics%DENSITY) &
        = 2.0 / (2.*pi*sigma**2) * exp(-r1/(2.*sigma**2)) &
        + 0.5 / (2.*pi*sigma**2) * exp(-r2/(2.*sigma**2)) &
        + 1.0 / (2.*pi*sigma**2) * exp(-r3/(2.*sigma**2))

      r1 = rsq(mesh,mesh%radius%faces(:,:,:,1),x1)
      r2 = rsq(mesh,mesh%radius%faces(:,:,:,1),x3)
      r3 = rsq(mesh,mesh%radius%faces(:,:,:,1),x4)

      timedisc%solution%data4d(:,:,:,1) &
        = - 2.0 / sqrt(r1) * erf(sqrt(r1 / 2.) / sigma) &
          - 0.5 / sqrt(r2) * erf(sqrt(r2 / 2.) / sigma) &
          - 1.0 / sqrt(r3) * erf(sqrt(r3 / 2.) / sigma)
      timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = 1.0
      timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.0
      timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = 0.0
    CASE(3)
      r1 = rsq(mesh,mesh%curv%bcenter,x1)
      r2 = rsq(mesh,mesh%curv%bcenter,x2)
      timedisc%pvar%data4d(:,:,:,physics%DENSITY)  = 0.02/(3.2*pi) &
        + 0.99 * (  exp(-0.5*(r1/sigma**2))/(2.*pi*sigma**2) &
                  + exp(-0.5*(r2/sigma**2))/(2.*pi*sigma**2))

      timedisc%pvar%data4d(:,:,:,physics%PRESSURE) = 0.02/(3.2*pi) &
        + gn / (2.*pi*sigma**2) &
           * (  ei(-(r1/sigma**2)) - ei(-0.5*(r1/sigma**2)) &
              + ei(-(r2/sigma**2)) - ei(-0.5*(r2/sigma**2)))

      timedisc%pvar%data4d(:,:,:,physics%XVELOCITY) = 0.0
      timedisc%pvar%data4d(:,:,:,physics%YVELOCITY) = r - r*omega
    END SELECT

    DO dir=west,east
      SELECT TYPE(bound => timedisc%Boundary%boundary(dir)%p)
      CLASS IS (boundary_noslip)
        DO k=mesh%KMIN,mesh%KMAX
          DO j=mesh%JMIN,mesh%JMAX
            DO ig=1,mesh%GNUM
              SELECT CASE (dir)
              CASE(west)
                  i = mesh%IMIN-ig
              CASE(east)
                  i = mesh%IMAX+ig
              END SELECT
              bound%data(ig,j,k,physics%YVELOCITY) &
                = timedisc%pvar%data4d(i,j,k,physics%YVELOCITY)
            END DO
          END DO
        END DO
      END SELECT
    END DO

    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
  END SUBROUTINE initdata

  FUNCTION rsq(Mesh,r,x) RESULT(res)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base) :: mesh
    REAL, DIMENSION(2) :: x
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: r
    REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: res
    !------------------------------------------------------------------------!
    INTENT(IN) :: x
    !------------------------------------------------------------------------!
    res = r**2 + x(1)**2 - 2.*r*x(1) * cos(mesh%curv%bcenter(:,:,:,2)-x(2))
  END FUNCTION rsq

END PROGRAM orbitingcylinders