| Version 6 (modified by , 7 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 ( and ), the temperature and density of the gas at infinity ( and ), and the percentage of solar mass the sink particle should have (). It then calculates the nondimensional radius for each zone of the grid, which is defined as:
where is the bondi radius calculated using the sound speed of the gas at infinity and the mass of the central object:
Note, strictly speaking the Bondi radius formula in the code also has a 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 () is calculated as a function of
Using this, the density in each cell is set using:
Once z is known, the nondimensional velocity follows from Eqn. 8 in Bondi 1952:
According to Bondi's scaling relations, , where is the local inward velocity of the gas, and is the sound speed at infinity. The inward radial momentum is then set in each zone using:
Since the gas is a polytrope (), 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 . 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 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 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
