Changes between Version 3 and Version 4 of BondiModule

Timestamp:

04/23/18 13:28:31 (7 years ago)

Author:

Erica Kaminski

Comment:

—

BondiModule

@@ -109 +109 @@
 
 
-
+{{{
+!#########################################################################
+!               
+!    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
+
+
+}}}