collapse.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: collapse.f90                                                      #
!#                                                                           #
!# Copyright (C) 2008-2024                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Björn Sperling   <sperling@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 collapse
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  REAL, PARAMETER    :: gn        = 6.6742e-11     ! [m^3/kg/s^2]
  ! simulation parameters
  REAL, PARAMETER    :: tsim      = 1.0e-0 ! simulation time in terms of the
                                           !   free-fall time [TAU]
  REAL, PARAMETER    :: gamma     = 1.4    ! ratio of specific heats
  REAL, PARAMETER    :: mass      = 1.0e+0 ! mass of spheroid
  REAL, PARAMETER    :: centmass  = 1.0e+2 ! central pointmass
  REAL, PARAMETER    :: rsph      = 30.0   ! semi-minor axis of the spheroid
  REAL, PARAMETER    :: ecc       = 0.0    ! eccentricity (0.0 is sphere)
  REAL, PARAMETER    :: vol0      = 4*pi/3 * rsph*rsph*rsph / (1.-ecc*ecc)
                                           ! volume of the spheroid
  REAL, PARAMETER    :: omega     = 1.0e-7 ! angular velocity (0.0 to disable)
  REAL, PARAMETER    :: omega_frame = 0.0*omega ! angular velocity (0.0 to disable)
  REAL, PARAMETER    :: eta_p     = 100.0  ! ratio of p_(hydro_static) to p 
                                           ! (in case of self-gravity) approx
                                           ! 100 (free fall); approx 1 stable
                                           ! (without pointmass)
  REAL, PARAMETER    :: eta_rho   = 1.0e-6 ! density ratio rho / rho_inf 
  ! mesh settings
!  INTEGER, PARAMETER :: MGEO = SPHERICAL   ! geometry of the mesh
  INTEGER, PARAMETER :: mgeo = cylindrical
!  INTEGER, PARAMETER :: MGEO = TANCYLINDRICAL
  INTEGER, PARAMETER :: xres = 150          ! x-resolution
  INTEGER, PARAMETER :: yres = 1 !XRES          ! y-resolution
  INTEGER, PARAMETER :: zres = xres           ! z-resolution
  REAL, PARAMETER    :: rmax = 1.5         ! size of comput. domain in [RSPH]
  REAL, PARAMETER    :: gpar = 0.5*rsph    ! geometry scaling parameter
  ! 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 = 'collapse_cyl_rhs1_150'
  !--------------------------------------------------------------------------!
  REAL               :: tau                ! free-fall time scale
  REAL               :: rho0               ! initial density of the sphere
  REAL               :: p0                 ! initial hydrostatic pressure
  !--------------------------------------------------------------------------!
  CLASS(fosite),ALLOCATABLE   :: sim
  LOGICAL  :: ok
  !--------------------------------------------------------------------------!
  ALLOCATE(sim)
  CALL sim%InitFosite()
 
#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_plan(1)
#ifdef PARALLEL
  END IF
#endif
 
  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()
  ok = .NOT.sim%aborted
 
#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
  tap_check(ok,"Finished simulation")
  tap_done
#ifdef PARALLEL
  END IF
#endif
 
  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),sgbc
    TYPE(dict_typ), POINTER :: mesh,physics, boundary, datafile, rotframe, &
                               sources, timedisc, fluxes, grav, pmass, selfgravity
    REAL              :: x1,x2,y1,y2,z1,z2
    !------------------------------------------------------------------------!
    INTENT(INOUT)     :: sim
    !------------------------------------------------------------------------!
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(spherical)
       x2 = rmax*rsph
       x1 = 2.*x2 / (xres+2)        ! x_min = 2*dx
       y1 = 0.005
       y2 = pi-0.005
       z1 = 0.0
       z2 = 2.*pi
    CASE(cylindrical)
       x2 = rmax*rsph
       x1 = 2.*x2 / (xres+2)        ! x_min = 2*dx
       y1 = 0.0
       y2 = 2.*pi
       z1 = -rmax*rsph
       z2 = rmax*rsph
!    CASE(OBLATE_SPHEROIDAL)
!       x2 = RMAX*RSPH/GPAR
!       x2 = LOG(x2+SQRT(x2**2-1.0)) ! = ACOSH(RMAX*RSPH/GPAR)
!       x1 = 2.*x2/(XRES+2)          ! x_min = 2*dx
!       y1 = -0.5*PI
!       y2 = 0.5*PI
    CASE(tancylindrical)
       x1 = atan(-rmax*rsph/gpar)
       x2 = atan(rmax*rsph/gpar)
       y2 = rmax*rsph
       y1 = 2.*y2 / (yres+2)        ! x_min = 2*dx
       z1 = 0.0
       z2 = 2.*pi
!    CASE(SINHSPHERICAL)
!       x2 = RMAX*RSPH/GPAR
!       x2 = LOG(x2+SQRT(x2**2+1.0)) ! = ASINH(RMAX*RSPH/GPAR)
!       x1 = 2.*x2/(XRES+2)          ! x_min = 2*dx
!       y1 = 0.0
!       y2 = PI
    CASE DEFAULT
       CALL sim%Error("InitProgram","geometry not supported for collapse simulation")
    END SELECT
 
    ! mesh settings
    mesh => dict("meshtype" / midpoint, &
           "geometry" / mgeo, &
           "inum"     / xres, &
           "jnum"     / yres, &
           "knum"     / zres, &
           "omega"    / omega_frame, &
           "use_fargo"/ 0,        &
           "fargo"    / 0,        &
           "xmin"     / x1, &
           "xmax"     / x2, &
           "ymin"     / y1, &
           "ymax"     / y2, &
           "zmin"     / z1, &
           "zmax"     / z2, &
           "output/radius" / 1, &
           "output/volume" / 1, &
           "gparam"   / gpar)
    ! mesh settings and boundary conditions
    SELECT CASE(mgeo)
    CASE(spherical)
       bc(west)  = reflecting !ABSORBING
       bc(east)  = reflecting !ABSORBING
       bc(south) = reflecting !AXIS
       bc(north) = reflecting !AXIS
       bc(bottom)= periodic
       bc(top)   = periodic
!       sgbc = SPHERMULTEXPAN        ! use spherical multipole expansion for BC
                                    !   in the multigrid poisson solver
    CASE(cylindrical)
       bc(west)  = axis       !ABSORBING
       bc(east)  = reflecting !ABSORBING
       bc(south) = periodic
       bc(north) = periodic
       bc(bottom)= reflecting
       bc(top)   = reflecting
!       sgbc = CYLINMULTEXPAN        ! cylindrical multipole expansion
!    CASE(OBLATE_SPHEROIDAL)
!       bc(WEST)  = ABSORBING
!       bc(EAST)  = ABSORBING
!       bc(SOUTH) = AXIS
!       bc(NORTH) = AXIS
!       sgbc = CYLINMULTEXPAN        ! cylindrical multipole expansion
    CASE(tancylindrical)
       bc(west)  = reflecting !ABSORBING
       bc(east)  = reflecting !ABSORBING
       bc(south) = reflecting !AXIS
       bc(north) = reflecting !NO_GRADIENTS
       bc(bottom)= periodic
       bc(top)   = periodic
!       sgbc = CYLINMULTEXPAN        ! cylindrical multipole expansion
!    CASE(SINHSPHERICAL)
!       bc(WEST)  = ABSORBING
!       bc(EAST)  = ABSORBING
!       bc(SOUTH) = AXIS
!       bc(NORTH) = AXIS
!       sgbc = SPHERMULTEXPAN        ! spherical multipole expansion
    CASE DEFAULT
       CALL sim%Error("InitProgram","geometry not supported for collapse simulation")
    END SELECT
 
    ! 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"   / gamma)                 ! ratio of specific heats        !
 
     ! compute some derived simulation parameters
    rho0 = mass / vol0                 ! initial density within the spheroid !
    ! "hydrostatic" pressure * ETA_P 
    !     => with ETA_P approx 100 => free-fall (in case of self-gravity)
    p0 = 4.0/3.0*pi*gn*rho0**2*rsph**2 / eta_p
    ! free-fall time (at radius RSPH) with contributions from both
    ! the selfgravitating spheroid and the central point mass
    tau = sqrt((rsph**3)/gn/(4./3.*pi*rsph**3*rho0 + centmass))
 
    ! flux calculation and reconstruction method
    fluxes => dict("order"     / linear, &
             "fluxtype"  / kt, &
             "variables" / primitive, &   ! vars. to use for reconstruction!
             "limiter"   / vanleer, &    ! one of: minmod, monocent,...   !
             "theta"     / 1.3)          ! optional parameter for limiter !
 
    ! source term due to a point mass
    pmass => dict("gtype" / pointmass, &        ! grav. accel. of a point mass   !
            "mass"  / centmass)          ! mass [kg]                      !
 
    ! source term due to self-gravity
!     selfgravity => Dict( &
!           "gtype"        / SPECTRAL)!, &   ! poisson solver for self-gravity!
!          "gtype"        / MULTIGRID, &   ! poisson solver for self-gravity!
!          "maxmult"      / 5, &           ! number of (spher.) multipol moments
!          "maxresidnorm" / 1.0E-7, &     ! accuracy of multigrid solver (max error)
!          "relaxtype"    / BLOCK_GAUSS_SEIDEL, &      ! relaxation method
!          "relaxtype"    / RED_BLACK_GAUSS_SEIDEL, &
!          "relaxtype"    / GAUSS_SEIDEL, &
!          "npre"         / 1, &           ! number of pre smoothings
!          "npost"        / 1, &           ! and post smoothings
!          "minres"       / 3, &           ! resolution of coarsest grid
!          "nmaxcycle"    / 250, &         ! limit for iterations
!          "bndrytype"    / sgbc)         ! multipole expansion (see above)
 
    ! enable gravity
    grav => dict("stype" / gravity, &
                "point" / pmass)
!     CALL SetAttr(grav, "self", selfgravity)  ! currently not supported because the multigird solver is not working
 
    ! initialize source term with gravitz module
    sources => dict("grav" / grav)
 
    ! check for rotating frame
    IF (omega_frame.GT.tiny(omega_frame)) THEN
      rotframe => dict("stype" / rotating_frame, &
               "x"     / 0.0, &
               "y"     / 0.0)
      CALL setattr(sources, "rotframe", rotframe)
    END IF
 
 
    ! time discretization settings
    timedisc => dict( &
         "method"    / modified_euler, &
         "order"     / 3, &
         "cfl"       / 0.4, &
         "stoptime"  / (tsim * tau), &
         "dtlimit"   / 1.0e-9, &
         "rhstype"   / 1, &
         "maxiter"   / 10000000)
 
    ! initialize data input/output
    datafile => dict( &
           "fileformat" / vtk, &
!            "fileformat" / XDMF, &
           "filename"   / (trim(odir) // trim(ofname)), &
           "count"      / onum)
 
    config => dict("mesh" / mesh, &
             "physics"  / physics, &
             "boundary" / boundary, &
             "fluxes"   / fluxes, &
             "timedisc" / timedisc, &
             "sources" / sources, &
             "datafile" / datafile)
  END SUBROUTINE makeconfig
 
 
  SUBROUTINE initdata(Mesh,Physics,Timedisc)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(mesh_base)    :: Mesh
    CLASS(physics_base) :: Physics
    CLASS(timedisc_base):: Timedisc
    !------------------------------------------------------------------------!
    CHARACTER(LEN=64) :: value
    INTEGER           :: i,j,k
    TYPE(marray_base) :: ez,posvec
    !------------------------------------------------------------------------!
    INTENT(IN)        :: mesh,physics
    INTENT(INOUT)     :: timedisc
    !------------------------------------------------------------------------!
    ! compute curvilinear components of vertical unit vector e_z
    ez = marray_base(3)
    ! set cartesian components
    ez%data2d(:,1:2) = 0.0
    ez%data2d(:,3)   = 1.0
    ! convert to curvilinear components
    CALL mesh%geometry%Convert2Curvilinear(mesh%bcenter,ez%data4d,ez%data4d)
 
    ! bary center position vector
    posvec = marray_base(3)
    posvec%data4d = mesh%posvec%bcenter
 
    SELECT TYPE(pvar => timedisc%pvar)
    TYPE IS(statevector_euler) ! non-isothermal HD
      ! density
      WHERE (mesh%radius%data2d(:,2).LE.rsph) ! bary center values collapsed to data1d
        pvar%density%data1d(:) = rho0
      ELSEWHERE
        pvar%density%data1d(:) = rho0 * eta_rho
      END WHERE
      ! velocities
      pvar%velocity = (omega-omega_frame) * (ez.x.posvec)
      ! pressure
      pvar%pressure%data1d(:) = p0
    CLASS DEFAULT
      CALL sim%Error("collapse::InitData","physics not supported in this setup")
    END SELECT
 
    CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
 
    CALL mesh%Info(" DATA-----> initial condition: " // &
         "Homogenious density with uniform angular frequency")
    IF (centmass.GT.tiny(centmass)) THEN
       WRITE (value,"(E12.4)") sqrt(rsph**3/(centmass*physics%constants%GN))
       CALL mesh%Info("                               " // &
          "timescale of pointmass:    " //trim(value))
    END IF
    IF (mass .GT. tiny(mass)) THEN
       WRITE(value,"(E12.4)") sqrt(3.0*pi/(4.0*rho0*physics%constants%GN))
       CALL mesh%Info("                               " // &
          "timescale of self-gravity: " //trim(value))
    END IF
    IF (omega .GT. tiny(omega)) THEN
       WRITE (value,"(E12.4)") omega
       CALL mesh%Info("                               " // &
          "angular velocity:          " // trim(value))
    END IF
 
    WRITE (value,"(E12.4)") tau
    CALL mesh%Info("                               " // &
          "free-fall time:            " //trim(value))
 
    SELECT TYPE(bwest => timedisc%Boundary%Boundary(west)%p)
    CLASS IS(boundary_custom)
      IF(mgeo.EQ.spherical) THEN
        CALL bwest%SetCustomBoundaries(mesh,physics, &
          (/custom_nograd,custom_reflect,custom_reflect,custom_fixed,custom_reflect/))
        bwest%data(:,:,:,physics%ZVELOCITY) = 0.0
      ELSE IF(mgeo.EQ.cylindrical) THEN
        CALL bwest%SetCustomBoundaries(mesh,physics, &
!          (/CUSTOM_NOGRAD,CUSTOM_REFLECT,CUSTOM_FIXED,CUSTOM_REFLECT/))
          (/custom_nograd,custom_reflect,custom_fixed,custom_nograd,custom_reflect/))
        bwest%data(:,:,:,physics%YVELOCITY) = 0.0
      END IF
    END SELECT
 
    SELECT TYPE(bbottom => timedisc%Boundary%Boundary(bottom)%p)
    CLASS IS(boundary_fixed)
      IF(mgeo.EQ.cylindrical) THEN
        bbottom%fixed(:,:,:) = .true.
        bbottom%data(:,:,1,:) = timedisc%pvar%data4d(:,:,mesh%KMIN,:)
        bbottom%data(:,:,2,:) = timedisc%pvar%data4d(:,:,mesh%KMIN,:)
      END IF
    CLASS IS(boundary_farfield)
        bbottom%data(:,:,1,:) = timedisc%pvar%data4d(:,:,mesh%KMIN,:)
        bbottom%data(:,:,2,:) = timedisc%pvar%data4d(:,:,mesh%KMIN,:)
    END SELECT
 
    SELECT TYPE(btop => timedisc%Boundary%Boundary(top)%p)
    CLASS IS(boundary_fixed)
      IF(mgeo.EQ.cylindrical) THEN
        btop%fixed(:,:,:) = .true.
        btop%data(:,:,1,:) = timedisc%pvar%data4d(:,:,mesh%KMAX,:)
        btop%data(:,:,2,:) = timedisc%pvar%data4d(:,:,mesh%KMAX,:)
      END IF
    CLASS IS(boundary_farfield)
        btop%data(:,:,1,:) = timedisc%pvar%data4d(:,:,mesh%KMAX,:)
        btop%data(:,:,2,:) = timedisc%pvar%data4d(:,:,mesh%KMAX,:)
    END SELECT
 
  END SUBROUTINE initdata
 
END PROGRAM collapse