vortex3d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: vortex2d.f90                                                      #
!#                                                                           #
!# Copyright (C) 2006-2024                                                   #
!# 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
  INTEGER, PARAMETER :: fargo= 0           ! 0: disables fargo transport
                                           ! 1: dynamic background velocity field,
                                           ! 2: fixed background velocity field
  ! 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.1      ! 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 => null()
  REAL, DIMENSION(:), ALLOCATABLE :: sum_numer, sum_denom, sigma
  LOGICAL :: ok
  INTEGER :: n,den,vel
  !--------------------------------------------------------------------------!
 
  ALLOCATE(sim)
  CALL sim%InitFosite()
 
#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_plan(3)
#ifdef PARALLEL
  END IF
#endif
 
  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%data1d(:) = sim%Timedisc%pvar%data1d(:)
  CALL sim%Run()
  ok = .NOT.sim%aborted
 
  ALLOCATE(sum_numer(sim%Physics%VNUM),sum_denom(sim%Physics%VNUM),sigma(sim%Physics%VNUM))
  ! use L1 norm to estimate the deviation from the exact stationary solution given
  ! in the initial condition:
  !   Σ |pvar - pvar_exact| / Σ |pvar_exact|
  DO n=1,sim%Physics%VNUM
    sum_numer(n) = sum(abs(sim%Timedisc%pvar%data2d(:,n)-pvar_init%data2d(:,n)),mask=sim%Mesh%without_ghost_zones%mask1d)
    sum_denom(n) = sum(abs(pvar_init%data2d(:,n)),mask=sim%Mesh%without_ghost_zones%mask1d)
  END DO
! PRINT *,"Rank: ",Sim%GetRank(),ACHAR(10),"numer=",sum_numer,ACHAR(10),"denom=",sum_denom,ACHAR(10),"sigma=",sigma
 
#ifdef PARALLEL
  IF (sim%GetRank().GT.0) THEN
    CALL mpi_reduce(sum_numer,sum_numer,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
    CALL mpi_reduce(sum_denom,sum_denom,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
  ELSE
    CALL mpi_reduce(mpi_in_place,sum_numer,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
    CALL mpi_reduce(mpi_in_place,sum_denom,sim%Physics%VNUM,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
  END IF
#endif
  WHERE (abs(sum_denom(:)).LE.2*tiny(sum_denom))
    sigma(:) = 0.0
  ELSEWHERE
    sigma(:) = sum_numer(:) / sum_denom(:)
  END WHERE
 
  den = sim%Physics%DENSITY
  vel = sim%Physics%YVELOCITY
 
#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
! PRINT *,"numer=",sum_numer,ACHAR(10)," denom=",sum_denom,ACHAR(10)," sigma=",sigma
tap_check_small(sigma(den),0.01,"density deviation < 1%")
tap_check_small(sigma(vel),0.1,"azimuthal velocity deviation < 10%")
tap_check(ok,"stoptime reached")
tap_done
! skip radial velocity deviation, because the exact value is 0
#ifdef PARALLEL
  END IF
#endif
 
  CALL sim%Finalize()
  DEALLOCATE(sim,pvar_init,sum_numer,sum_denom,sigma)
 
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(cylindrical)
      x1 = rmin
      x2 = rmax
      y1 =  0.0
      y2 =  2.0*pi
      IF (zres.GT.1) THEN
        z1 =  zmin
        z2 =  zmax
      ELSE
        z1 = 0
        z2 = 0
      END IF
      bc(west)   = axis
      bc(east)   = reflecting
      bc(south)  = periodic
      bc(north)  = periodic
      bc(bottom) = periodic
      bc(top)    = periodic
    CASE DEFAULT
       CALL sim%Error("vortex3d::MakeConfig","geometry currently not supported")
    END SELECT
 
    !mesh settings
    mesh => dict( &
              "meshtype" / midpoint, &
              "geometry" / mgeo,     &
              "omega"    / omega,    &
              "fargo/method" / fargo, &
              "decomposition"   / (/ -1, 1, -1/), & ! do not decompose along 2nd dimension with FARGO!
              "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,          &
         "rhstype"  / 0)
 
    ! 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(:,:,:)-mesh%OMEGA)*dist_axis(:,:,:)*ephi(:,:,:,n)
      END DO
      SELECT CASE(mesh%fargo%GetType())
      CASE(1,2)
        SELECT CASE(mesh%fargo%GetDirection())
        CASE(1)
          timedisc%w(:,:) = pvar%velocity%data4d(mesh%IMIN,:,:,1)
        CASE(2)
          timedisc%w(:,:) = pvar%velocity%data4d(:,mesh%JMIN,:,2)
        CASE(3)
          timedisc%w(:,:) = pvar%velocity%data4d(:,mesh%KMIN,:,3)
        END SELECT
      END SELECT
    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(:,:,:)-mesh%OMEGA)*dist_axis(:,:,:)*ephi(:,:,:,n)
      END DO
      SELECT CASE(mesh%fargo%GetType())
      CASE(1,2)
        SELECT CASE(mesh%fargo%GetDirection())
        CASE(1)
          timedisc%w(:,:) = pvar%velocity%data4d(mesh%IMIN,:,:,1)
        CASE(2)
          timedisc%w(:,:) = pvar%velocity%data4d(:,mesh%JMIN,:,2)
        CASE(3)
          timedisc%w(:,:) = pvar%velocity%data4d(:,:,mesh%KMIN,3)
        END SELECT
      END SELECT
    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