bondi3d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: bondi3d.f90                                                       #
!#                                                                           #
!# Copyright (C) 2006-2018                                                   #
!# 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 bondi3d
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  ! general constants
  REAL, PARAMETER    :: msun    = 1.989e+30   ! solar mass [kg]              !
  REAL, PARAMETER    :: gn      = 6.67408e-11 ! gravitional constant (SI)    !
  ! simulation parameters
  REAL, PARAMETER    :: tsim    = 20.0        ! simulation time [TAU]        !
  REAL, PARAMETER    :: accmass = 1.0*msun    ! mass of the accreting object !
  REAL, PARAMETER    :: gamma   = 1.4         ! ratio of specific heats      !
  ! boundary conditions
  REAL, PARAMETER    :: rhoinf  = 3.351e-17   ! density at infinity [kg/m^3] !
  REAL, PARAMETER    :: csinf   = 537.0       ! sound speed at infinity [m/s]!
  ! mesh settings
  INTEGER, PARAMETER :: mgeo    = spherical   ! geometry of the mesh         !
  INTEGER, PARAMETER :: xres    = 100          ! x-resolution                 !
  INTEGER, PARAMETER :: yres    = 4           ! y-resolution                 !
  INTEGER, PARAMETER :: zres    = 4           ! z-resolution                 !
  REAL, PARAMETER    :: rin     = 0.1         ! inner/outer radii in terms of!
  REAL, PARAMETER    :: rout    = 2.0         ! the Bondi radius RB, ROUT > 1!
  ! 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  = 'bondi3d'
  ! some derives quandities
  REAL               :: rb                    ! Bondi radius
  REAL               :: tau                   ! free fall time scale
  !--------------------------------------------------------------------------!
  CLASS(fosite), ALLOCATABLE :: sim
  REAL, DIMENSION(:), ALLOCATABLE :: sigma
  REAL :: sum_numer, sum_denom
  INTEGER :: n,den,vel,pre
  LOGICAL :: ok
  !--------------------------------------------------------------------------!

  ALLOCATE(sim)
  CALL sim%InitFosite()

#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_plan(4)
#ifdef PARALLEL
  END IF
#endif

  CALL makeconfig(sim, sim%config)
!  CALL PrintDict(Sim%config)
  CALL sim%Setup()
  ! set initial condition
  CALL initdata(sim%Mesh,sim%Physics,sim%Timedisc)
  ! compute numerical solution
  CALL sim%Run()
  ok = .NOT.sim%aborted
  ! compare with exact solution if requested
  IF (ASSOCIATED(sim%Timedisc%solution)) THEN
    ALLOCATE(sigma(sim%Physics%VNUM))
    DO n=1,sim%Physics%VNUM
      ! use L1 norm to estimate the deviation from the exact solution:
      !   Σ |pvar - pvar_exact| / Σ |pvar_exact|
      sum_numer = sum(abs(sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
                           sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n) &
                        -sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
                           sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
      sum_denom = sum(abs(sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
                           sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
#ifdef PARALLEL
      IF (sim%GetRank().GT.0) THEN
        CALL mpi_reduce(sum_numer,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
        CALL mpi_reduce(sum_denom,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
      ELSE
        CALL mpi_reduce(mpi_in_place,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
        CALL mpi_reduce(mpi_in_place,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
#endif
      sigma(n) = sum_numer / sum_denom
#ifdef PARALLEL
      END IF
#endif
    END DO
  ELSE
    sigma(:) = 0.0
  END IF

  den = sim%Physics%DENSITY
  vel = sim%Physics%XVELOCITY
  pre = sim%Physics%PRESSURE

#ifdef PARALLEL
  IF (sim%GetRank().EQ.0) THEN
#endif
tap_check(ok,"stoptime reached")
! These lines are very long if expanded. So we can't indent it or it will be cropped.
tap_check_small(sigma(den),4.0e-02,"density deviation < 4%")
tap_check_small(sigma(vel),4.0e-02,"radial velocity deviation < 4%")
! skip azimuthal velocity deviation, because exact value is 0
tap_check_small(sigma(pre),6.0e-02,"pressure deviation < 6%")
tap_done
#ifdef PARALLEL
  END IF
#endif

  IF (ALLOCATED(sigma)) DEALLOCATE(sigma)
  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)
    TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
                               grav, pmass, timedisc, fluxes
    REAL                    :: x1,x2,y1,y2,z1,z2
    !------------------------------------------------------------------------!
    INTENT(INOUT)           :: sim
    !------------------------------------------------------------------------!
    ! derived constants
    rb  = gn * accmass / csinf**2            ! bondi radius [m]              !
    tau = rb / csinf                         ! free fall time scale [s]      !

    ! geometry dependent setttings
    SELECT CASE(mgeo)
    CASE(spherical)
      x1 = rin * rb
      x2 = rout * rb
      y1 = 0.0
      y2 = pi
      z1 = 0.0
      z2 = 2*pi
      bc(west)   = absorbing
      bc(east)   = custom
      bc(south)  = axis
      bc(north)  = axis
      bc(bottom) = periodic
      bc(top)    = periodic
    CASE DEFAULT
       CALL sim%Physics%Error("bondi3d::MakeConfig","geometry currently not supported")
    END SELECT

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

    ! 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)

    ! flux calculation and reconstruction method
    fluxes => dict( &
          "order"       / linear, &
          "fluxtype"    / kt, &
          "variables"   / primitive, &
          "limiter"     / monocent, &
          "theta"       / 1.2)

    ! gravity term due to a point mass
    pmass => dict( &
          "gtype"       / pointmass, &
          "mass"        / accmass, &
          "outbound"    / 0)

   ! source term due to all gravity terms
    grav => dict( &
          "stype"       / gravity, &
          "pmass"       / pmass)

    ! time discretization settings
    timedisc => dict( &
          "method"      / modified_euler, &
          "order"       / 3, &
          "cfl"         / 0.4, &
          "stoptime"    / (tsim * tau), &
          "dtlimit"     / (1.0e-6 * tau), &
          "output/solution" / 1, &
          "maxiter"     / 20000)

    ! 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, &
          "sources/grav"  / grav, &
          "timedisc"    / timedisc, &
          "datafile"    / datafile)
  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
    REAL              :: r,rho,vr,cs2
    INTEGER           :: i,j,k,l
    CHARACTER(LEN=64) :: info_str
    !------------------------------------------------------------------------!
    ! initial condition
    SELECT TYPE(pvar => timedisc%pvar)
    TYPE IS(statevector_euler)
      ! constant density and pressure, vanishing velocity
      pvar%density%data1d(:)  = rhoinf
      pvar%velocity%data1d(:) = 0.
      pvar%pressure%data1d(:) = rhoinf * csinf**2 / gamma
    CLASS DEFAULT
      CALL physics%Error("bondi3d::InitData","only non-isothermal HD supported")
    END SELECT

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

    ! store exact stationary solution if requested,
    ! i.e. output/solution == 1 in timedisc config
    IF (ASSOCIATED(timedisc%solution)) THEN
      ! compute the stationary 2D planar bondi solution
      SELECT TYPE(pvar => timedisc%solution)
      TYPE IS(statevector_euler)
        DO k=mesh%KMIN,mesh%KMAX
          DO j=mesh%JMIN,mesh%JMAX
            DO i=mesh%IMIN,mesh%IMAX
              r = mesh%radius%bcenter(i,j,k)
              CALL bondi(r/rb,gamma,rhoinf,csinf,rho,vr)
              cs2 = csinf**2 * (rho/rhoinf)**(gamma-1.0)
              pvar%density%data3d(i,j,k)  = rho
              ! set the minimum of the exact velocity to some value > 0
              ! for comparison with the numerical results
!NEC$ UNROLL(3)
              DO l=1,physics%VDIM
                pvar%velocity%data4d(i,j,k,l) = sign(&
                  max(epsilon(cs2),abs(vr*mesh%posvec%bcenter(i,j,k,l)/r)), &
                    vr*mesh%posvec%bcenter(i,j,k,l))
              END DO
              pvar%pressure%data3d(i,j,k) = rho * cs2 / gamma
            END DO
          END DO
        END DO
      END SELECT
    END IF

    ! boundary condition: subsonic inflow according to Bondi's solution
    ! calculate Bondi solution for y=ymin..ymax at xmax
    SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
    CLASS IS (boundary_custom)
      ! this tells the boundary routine the actual boundary condition for each variable
      CALL beast%SetCustomBoundaries(mesh,physics, &
        (/custom_fixed,custom_nograd,custom_fixed,custom_fixed,custom_fixed/))
      DO k=mesh%KMIN,mesh%KMAX
        DO j=mesh%JMIN,mesh%JMAX
          DO i=1,mesh%GNUM
            ! get distance to the origin for each boundary cell
            r = mesh%radius%bcenter(mesh%IMAX+i,j,k)
            CALL bondi(r/rb,gamma,rhoinf,csinf,rho,vr)
            cs2 = csinf**2 * (rho/rhoinf)**(gamma-1.0)
            beast%data(i,j,k,physics%DENSITY)   = rho
            DO l=1,physics%VDIM
              beast%data(i,j,k,physics%XVELOCITY+l-1) &
                = vr*mesh%posvec%bcenter(mesh%IMAX+i,j,k,l)/r
            END DO
            beast%data(i,j,k,physics%PRESSURE)  = rho * cs2 / gamma
          END DO
        END DO
      END DO
    END SELECT

    WRITE(info_str,"(ES9.3)") rb
    CALL mesh%Info("                               " // "Bondi radius:       " &
         // trim(info_str) // " m")
    WRITE(info_str,"(ES9.3)") tau
    CALL mesh%Info("                               " // "Free fall time:     " &
         // trim(info_str) // " s")
  END SUBROUTINE initdata


  SUBROUTINE bondi(r,gamma,rhoinf,csinf,rho,vr)
    USE roots
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    ! computes the Bondi solution for spherically symmetric accretion        !
    ! uses funcd, GetRoot                                                    !
    ! INPUT paramter:                                                        !
    !   r     : radius in units of the Bondi radius r_b = G*M/csinf**2       !
    !   gamma : ratio of specific heats (1 < gamma < 5/3)                    !
    !   rhoinf: density at infinity                                          !
    !   csinf : speed of sound at infinity                                   !
    ! OUTPUT paramter:                                                       !
    !   rho   : density @ r                                                  !
    !   vr    : radial velocity @ r                                          !
    !------------------------------------------------------------------------!
    REAL, INTENT(IN)  :: r,gamma,rhoinf,csinf
    REAL, INTENT(OUT) :: rho,vr
    !------------------------------------------------------------------------!
    INTERFACE
       PURE SUBROUTINE funcd(x,fx,dfx,plist)
         IMPLICIT NONE
         REAL, INTENT(IN)  :: x
         REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
         REAL, INTENT(OUT) :: fx,dfx
       END SUBROUTINE funcd
    END INTERFACE
    !------------------------------------------------------------------------!
    REAL, PARAMETER :: xacc = 1.0e-6     ! accuracy for root finding
    REAL :: gp1,gm1,g35,rc,chi,lambda,psi,gr
    INTEGER :: err
    COMMON /funcd_parameter/ gm1, gr
    !------------------------------------------------------------------------!
    ! for convenience
    gm1 = gamma - 1.0
    gp1 = gamma + 1.0
    g35 = 0.5 * (5.0-3.0*gamma)

    ! critical radius
    rc = 0.5 * g35
    ! critical dimensionless accretion rate
    lambda = 0.25 * g35**(-g35/gm1)

    ! Newton-Raphson to solve Bondis equations for psi
    chi = r**2 / lambda
    gr  = chi**(2.*gm1/gp1) * (1./gm1 + 1./r)
    IF (r.LT.rc) THEN
       CALL getroot_newton(funcd,1.0,gr,psi,err)
    ELSE
       CALL getroot_newton(funcd,1.0e-6,1.0,psi,err)
    END IF

    ! return values
    rho = rhoinf * chi**(-2./gp1) / psi        ! density
    vr  = -csinf * chi**(-gm1/gp1) * psi       ! radial velocity
  END SUBROUTINE bondi


END PROGRAM bondi3d


! for exact Bondi solution at the outer boundary
PURE SUBROUTINE funcd(y,fy,dfy,plist)
  IMPLICIT NONE
  !------------------------------------------------------------------------!
  REAL, INTENT(IN)  :: y
  REAL, INTENT(OUT) :: fy,dfy
  REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
  !------------------------------------------------------------------------!
  REAL :: gm1,gr
  COMMON /funcd_parameter/ gm1,gr
  !------------------------------------------------------------------------!
  fy  = 0.5*y*y + y**(-gm1) / gm1 - gr
  dfy = y - y**(-gm1-1.)
END SUBROUTINE funcd