bondi2d.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: bondi2d.f90                                                       #
!#                                                                           #
!# Copyright (C) 2006-2018                                                   #
!# Tobias Illenseer <tillense@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 bondi2d
  USE fosite_mod
#include "tap.h"
  IMPLICIT NONE
  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
  !--------------------------------------------------------------------------!
  ! 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] (free fall)
  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  = 1.0e-20     ! density at infinity [kg/m^3]
  REAL, PARAMETER    :: csinf   = 1.0e+04     ! sound speed at infinity [m/s]
  ! mesh settings
  INTEGER, PARAMETER :: mgeo    = cylindrical ! geometry
  INTEGER, PARAMETER :: xres    = 100         ! x-resolution
  INTEGER, PARAMETER :: yres    = 6           ! y-resolution
  INTEGER, PARAMETER :: zres    = 1           ! z-resolution
  REAL, PARAMETER    :: rin     = 0.1         ! inner/outer radii in terms of
  REAL, PARAMETER    :: rout    = 2.0         !   the Bondi radius RB, ROUT > 1
  REAL, PARAMETER    :: gpar    = 1.0         ! geometry scaling parameter in [RB]
  ! 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  = 'bondi2d'
  ! 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(config)
  CALL sim%Setup()
  ! set initial condition
  CALL initdata(sim%Mesh, sim%Physics, sim%Fluxes, sim%Timedisc)
  ! run the simulation
  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),1.0e-02,"density deviation < 1%")
tap_check_small(sigma(vel),1.0e-02,"radial velocity deviation < 1%")
! skip azimuthal velocity deviation, because exact value is 0
tap_check_small(sigma(pre),1.0e-02,"pressure deviation < 1%")
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), INTENT(INOUT) :: 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
    !------------------------------------------------------------------------!
    ! derived constants
    rb  = gn * accmass / csinf**2           ! bondi radius [m]               !
    tau = rb / csinf                        ! free fall time scale [s]       !

    ! mesh settings
    SELECT CASE(mgeo)
    CASE(cylindrical)
       x1 = rin * rb
       x2 = rout * rb
       y1 = 0.0
       y2 = 2*pi
       z1 = 0.0
       z2 = 0.0
       bc(west)   = no_gradients
       bc(east)   = custom
       bc(south)  = periodic
       bc(north)  = periodic
       bc(bottom) = no_gradients
       bc(top)    = no_gradients
    CASE DEFAULT
       CALL sim%Physics%Error("bondi2d::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, &
             "gparam"    / (gpar*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" / conservative, &
             "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"   / 1000000)

    ! 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,Fluxes,Timedisc)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(physics_base),  INTENT(IN)    :: Physics
    CLASS(mesh_base),     INTENT(IN)    :: Mesh
    CLASS(fluxes_base),   INTENT(IN)    :: Fluxes
    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("bondi2d::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/))
      DO k=mesh%KMIN,mesh%KMAX
        DO j=mesh%JMIN,mesh%JMAX
          DO i=1,mesh%GINUM
            ! 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

    ! print some information
    CALL mesh%Info(" DATA-----> initial condition: 2D Bondi accretion")
    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 rotationally symmetric accretion in    !
    !   planar geometry                                                      !
    ! uses funcd, GetRoot_newton                                             !
    ! 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
    !------------------------------------------------------------------------!
    REAL, PARAMETER :: xacc = 1.0e-6     ! accuracy for root finding
    REAL :: gp1,gm1,rc,chi,lambda,psi,gr
    REAL :: params(2)
    INTEGER :: err
    !------------------------------------------------------------------------!
    ! for convenience
    gm1 = gamma - 1.0
    gp1 = gamma + 1.0

    ! critical radius
    rc = (3.0-gamma) / 2.0

    ! critical dimensionless accretion rate
    lambda = rc**(-rc/gm1)

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

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

END PROGRAM bondi2d

! find the root of fy to compute the exact Bondi solution
PURE SUBROUTINE funcd(y,fy,dfy,plist)
  IMPLICIT NONE
  !------------------------------------------------------------------------!
  REAL, INTENT(IN)  :: y
  REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
  REAL, INTENT(OUT) :: fy,dfy
  !------------------------------------------------------------------------!
  ! plist(1) = gamma - 1
  ! plist(2) = gr
  !------------------------------------------------------------------------!
  fy  = 0.5*y*y + y**(-plist(1)) / plist(1) - plist(2)
  dfy = y - y**(-plist(1)-1.)
END SUBROUTINE funcd