wiki:BondiModule

Context Navigation

Version 5 (modified by Erica Kaminski, 8 years ago) ( diff )
—

Bondi Module

Description

The Bondi module (located in astrobear/src/modules/) was originally designed to be used by the accretion module to enable sink particles to accrete. It has since been expanded to include the ability to initialize the domain with a Bondi solution, as well as impose the Bondi solution at an inner spherical boundary surrounding a sink particle, or over a spherical boundary at the edge of the domain. For the details of the Bondi solution see Bondi's 1952 paper (1952MNRAS.112..195B).

Initializing the Bondi Module

The problem.data file contains the following:

ProblemData

mcent: The mass of the sink particle in fraction of solar masses (e.g. mcent=0.4 would be 40% solar mass)
namb: The number density of the ambient medium in units of cm^-3.
tamb: The temperature of the ambient medium in units of K.
IBS: The radius of an inner spherical boundary in computational units, within which the Bondi solution will be pasted (this can be shut off in the problem module, as described below).
OBS: The radius of an outer boundary sphere in computational units, outside of which the Bondi solution will be pasted.

Initializing the Computational Domain

The Bondi module takes as input an inner and outer spherical boundary radius (IBS and OBS), the temperature and density of the gas at infinity (Tamb and namb), and the percentage of solar mass the sink particle should have (mcent). It then calculates the nondimensional radius for each zone of the grid, which is defined as:

\chi = r/r_{BH}

where r_{BH} is the bondi radius calculated using the sound speed of the gas at infinity and the mass of the central object:

r_{BH}=\frac{GM*amu*\mu}{\gamma K_B T_{\infty}}

Note, strictly speaking the Bondi radius formula in the code also has a v_\infty term, the velocity of the gas at infinity. This is zero in the present case so it is ignored in the above expression.

Next, the non dimensional density profile (z=n/n_{amb}) is calculated as a function of \chi

Using this, the density in each cell is set using:

info\%q(i,j,k,1)=z*n_{amb}

Once z is known, the nondimensional velocity follows from Eqn. 8 in Bondi 1952:

y=\frac{\lambda}{\chi^2 z}

According to Bondi's scaling relations, y= v/C_{\infty}, where v is the local inward velocity of the gas, and C_\infty is the sound speed at infinity. The inward radial momentum is then set in each zone using:

Info\%(:,:,:,px)=-y C_{\infty} \rho \frac{x}{r}

Since the gas is a polytrope (P\propto \rho ^\gamma), the energy of each cell is then easily calculated (the total energy is a sum of radial kinetic energy + internal energy).

The above solution is pasted everywhere in the grid (including ghost zones) except close to the origin where the speeds become very large. Within the inner spherical boundary the solution is just set to be constant — using the Bondi solution at nondimensional radius \chi=IBS/r_{BH}. Note, in the problem module copied below this functionality is actually turned off, and the Bondi solution is copied to the mesh everywhere in the domain. This is because under most circumstances the code should be able to handle the high speeds generated by the Bondi solution close to the origin.

At each subsequent cycle after initialization, the bondi flow is only set in the ghost zones (again with the option of stepping on the zones within IBS if desired).

Initializing the Central Sink Particle

The Bondi module initializes a sink particle at the center of the Bondi flow corresponding to the origin. This particle is set to be fixed at the origin so that it cannot move over the course of the simulation (this is a safety check — given the spherical symmetry of the problem the particle should not move on its own). Its mass is set using the input variable mcent defined above. Its accretion routine is set Krumholz accretion, which is AstroBEAR's closest prescription to Bondi accretion.

Multiphysics

The only physics turned on in this module on top of the hydrodynamics is point gravity, associated with the sink particle. The gas in the mesh is not self-gravitating, but rather, only feels the force of the particle which causes its collapse.

Module

Without further adieu, the Bondi module in its most recent form is copied below. I have tried to comment the module heavily, but if the reader should have any questions, feel free to email the astrobear developer team with questions (astrobear_dev@pas.rochester.edu).

!#########################################################################
!		
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    problem.f90 of module Bondi is part of AstroBEAR.
!
!    AstroBEAR 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 3 of the License, or    
!    (at your option) any later version.
!
!    AstroBEAR 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.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with AstroBEAR.  If not, see <http://www.gnu.org/licenses/>.
!
!#########################################################################
!Bondi Module

MODULE Problem

  USE DataDeclarations
  USE GlobalDeclarations
  USE PhysicsDeclarations
  USE ParticleDeclarations
  USE CommonFunctions
  USE Bondi
  USE Winds
  IMPLICIT NONE
  PRIVATE 

  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  TYPE(PointGravityDef), POINTER :: PointGravityObj
  REAL(KIND=qprec) ::   namb, tamb, ibs, obs, mcent, r_bh, ramb, mdot_bondi, camb, r_sonic!, mdot_ruffert


CONTAINS

  SUBROUTINE ProblemModuleInit
    INTEGER :: iErr
    TYPE(ParticleDef), POINTER :: Particle
    NAMELIST /ProblemData/ namb, tamb, ibs, obs, mcent

    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data')
    READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
    CLOSE(PROBLEM_DATA_HANDLE, IOSTAT=iErr)

    ramb = (namb*amu*xmu)/rscale !convert from namb to mass density in computational units
    camb=sqrt(gamma*tamb/tempscale) !calculate scaled sound speed of gas at infinity


    IF (.not. lRestart) THEN
       NULLIFY(Particle)
       CALL CreateParticle(Particle)
       Particle%q(1)=mcent*MSolar/rScale/lScale**3.0
       Particle%xloc=0
       Particle%iAccrete=2 !0=no accretion, 1=fed, 2=!KRUMHOLZ_ACCRETION
       Particle%lFixed=.true.
       particle%buffer(0:MaxLevel)=ceiling(ibs/levels(0:MaxLevel)%dx)
       CALL AddSinkParticle(Particle)
       CALL CreatePointGravityObject(Particle%PointGravityObj)
       Particle%PointGravityObj%soft_length=1d0*levels(MaxLevel)%dx
       Particle%PointGravityObj%soft_function=SPLINESOFT
       Particle%PointGravityObj%Mass=mcent*MSolar/rScale/lScale**3.0 !Central object mass in computational units
    ELSE !the next line is needed for restarts using only one particle, such as this module:
       particle=>SinkParticles%self  !must redefine the particle pointer upon restarts -- for multiple particles, will need to iterate over them (see the binary problem for ex)
    END IF

    !Calculate theoretical accretion rate (computational units):
    mdot_bondi=(4.0*pi*bondi_lambda*ramb*(ScaleGrav*Particle%q(1))**2.0) / (camb**3.0)
    !Calculate the sonic radius (computational units):
    r_sonic = ((5.0-3.0*gamma)/4.0)*((ScaleGrav*particle%q(1))/camb**2.0)

    IF (MPI_ID == 0) THEN
       PRINT *, 'cs_inf=', camb
       PRINT *, 'sonic radius=', r_sonic
       PRINT *, 'T_inf=', Tamb/TempScale
       PRINT *, "rho_inf=", ramb
       PRINT *, "mdot_bondi (CU) = ", mdot_bondi
       PRINT *, 'M_star=', particle%q(1)
    END IF


    !Calculate Bondi radius of initial particle mass
    r_BH=ScaleGrav*Particle%q(1)/(gamma*tamb/TempScale)
    IF (MPI_ID == 0) THEN
       write(*,*) 'Bondi radius (CU)= ', r_BH
    END IF

  END SUBROUTINE ProblemModuleInit

  !> Initial Conditions
  !! @param Info Info object
  SUBROUTINE ProblemGridInit(Info)
    !! @brief Initializes the grid data according to the requirements of the problem.
    !! @param Info A grid structure.	
    TYPE (InfoDef) :: Info
    INTEGER :: i,j,k
    INTEGER :: rmbc,zrmbc,level
    INTEGER :: mx, my, mz
    INTEGER :: iErr
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:) :: q
    REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,r
    REAL(KIND=qPREC) :: x_, y_, z_, rho, v, temp
    level=Info%level
    q=>Info%q
    ! Calculating the number of ghost cells on each side of the grid.
    rmbc=levels(level)%gmbc(levels(level)%step)
    dx=levels(level)%dX
    dy=dx
    SELECT CASE(nDim)
    CASE(2)
       zrmbc=0
       dz=0d0
    CASE(3)
       zrmbc=rmbc
       dz=dx
    END SELECT
    mx = Info%mX(1)
    my = Info%mX(2)
    mz = Info%mX(3)
    xl=Info%xBounds(1,1)
    yl=Info%xBounds(2,1)
    zl=Info%xBounds(3,1)

    !All of the values below are in computational units
    ! need the origin to be at 0,0,0:
    DO i=1-rmbc, mx+rmbc  
       x = (xl+(REAL(i,xPrec)-half)*dx) 
       DO j=1-rmbc, my+rmbc  
          y = (yl+(REAL(j,xPrec)-half)*dy)
          DO k=1-zrmbc, mz+zrmbc 
             z = (zl+(REAL(k,xPrec)-half)*dz)
             r = sqrt(x**2 + y**2+z**2)
            
             ! Calculate non-dimensional radius(x_), density(z_), and velocity(y_
             !x_=max(r, ibs-5d0*levels(0)%dx)/r_BH --- commented to avoid Bondi solution being copied within IBS
             x_=r/r_BH 
             z_=BH_alpha(x_) !nondimensional density (normalized to rho_inf)
             y_=Bondi_lambda/(z_*x_**(myDim-1)) !nondimensional velocity (normalized to Cs_inf)
          

             ! Then calculate physical values
             rho=z_*namb/nScale
             v=-y_*sqrt(gamma*tamb/TempScale)
             temp=tamb*(z_**(gamma1))/TempScale

             q(i,j,k,1)=rho
             q(i,j,k,ivx)=rho*v*x/r
             q(i,j,k,ivy)=rho*v*y/r
             q(i,j,k,ivz)=rho*v*z/r
             q(i,j,k,iE)=half*rho*v**2+gamma7*rho*temp
          END DO
       END DO
    END DO
  END SUBROUTINE ProblemGridInit

  !! @param Info Info object
  SUBROUTINE ProblemBeforeStep(Info)
    !! @brief Performs any tasks required before the advance step.
    !! @param Info A grid structure.	
    TYPE (InfoDef) :: Info
    TYPE(WindDef), POINTER :: Wind
    INTEGER :: i,j,k, iErr
    INTEGER :: rmbc,zrmbc,level
    INTEGER :: mx, my, mz, edge
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:) :: q
    REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,r
    REAL(KIND=qPREC) :: x_, y_, z_, rho, v, temp


  !Set up wind object which prevents any material from leaving the grid (i.e. Ghost zones setup):
    DO i=1,nDim
       DO edge=1,2
          IF (Gmthbc(i,edge) == 1) THEN 
             CALL CreateWind(Wind)
             Wind%density=ramb
             wind%temperature=tamb/tempscale
             wind%type=OUTFLOW_ONLY
             Wind%dir=i
             Wind%edge=edge
          END IF
       END DO
    END DO

    level=Info%level
    q=>Info%q
    ! Calculating the number of ghost cells on each side of the grid.
    rmbc=levels(level)%gmbc(levels(level)%step)
    dx=levels(level)%dX
    dy=dx
    SELECT CASE(nDim)
    CASE(2)
       zrmbc=0
       dz=0d0
    CASE(3)
       zrmbc=rmbc
       dz=dx
    END SELECT
    mx = Info%mX(1)
    my = Info%mX(2)
    mz = Info%mX(3)
    xl=Info%xBounds(1,1)
    yl=Info%xBounds(2,1)
    zl=Info%xBounds(3,1)

    !All of the values below are in computational units

    DO i = 1-rmbc, mx+rmbc
       x = (xl+(REAL(i,xPrec)-half)*dx) 
       DO j = 1-rmbc, my+rmbc  
          y = (yl+(REAL(j,xPrec)-half)*dy)
          DO k = 1-zrmbc,mz+zrmbc 
             z = (zl+(REAL(k,xPrec)-half)*dz)

             r = sqrt(x**2 + y**2 +z**2)

             !Commented following line so that only solution copied to outer boundary sphere:
             !IF (r < ibs .OR. r > obs) THEN !Set cells to analytic solution
             IF (r > obs) THEN !Set cells to analytic solution
                !same with following line:
                !x_=max(ibs-5d0*levels(0)%dx, r)/r_BH !Adjust values deep inside inner region to avoid extremely high velocities etc...
                x_=r/r_BH 
                z_=BH_alpha(x_)
                y_=Bondi_lambda/(z_*x_**(myDim-1))
                rho=z_*namb/nScale
                v=-y_*sqrt(gamma*tamb/TempScale)
                temp=tamb*(z_**(gamma1))/TempScale

                q(i,j,k,1)=rho
                q(i,j,k,ivx)=rho*v*x/r
                q(i,j,k,ivy)=rho*v*y/r
                !erica added z component:
                q(i,j,k,ivz)=rho*v*z/r
                q(i,j,k,iE)=half*rho*v**2+gamma7*rho*temp

             END IF



          END DO
       END DO
    END DO

  END SUBROUTINE ProblemBeforeStep

  !> Does nothing
  !! @param Info Info object
  SUBROUTINE ProblemAfterStep(Info)
    !! @brief Performs any post-step corrections that are required.
    !! @param Info A grid structure.	
    TYPE (InfoDef) :: Info
    CALL ProblemBeforeStep(Info)
  END SUBROUTINE ProblemAfterStep

  !> Does nothing
  !! @param Info Info object
  SUBROUTINE ProblemSetErrFlag(Info)
    !! @brief Sets error flags according to problem-specific conditions..
    !! @param Info A grid structure.	
    TYPE (InfoDef) :: Info
  END SUBROUTINE ProblemSetErrFlag

  SUBROUTINE ProblemBeforeGlobalStep(n)
    INTEGER :: n
  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem

Note: See TracWiki for help on using the wiki.

Download in other formats:

Plain Text