u/erica/AccretionModelingBlog: problem.f90_steadystate

File problem.f90_steadystate, 8.6 KB (added by Erica Kaminski, 7 years ago)
Line 
1!#########################################################################
2!
3! Copyright (C) 2003-2012 Department of Physics and Astronomy,
4! University of Rochester,
5! Rochester, NY
6!
7! problem.f90 of module Bondi is part of AstroBEAR.
8!
9! AstroBEAR is free software: you can redistribute it and/or modify
10! it under the terms of the GNU General Public License as published by
11! the Free Software Foundation, either version 3 of the License, or
12! (at your option) any later version.
13!
14! AstroBEAR is distributed in the hope that it will be useful,
15! but WITHOUT ANY WARRANTY; without even the implied warranty of
16! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17! GNU General Public License for more details.
18!
19! You should have received a copy of the GNU General Public License
20! along with AstroBEAR. If not, see <http://www.gnu.org/licenses/>.
21!
22!#########################################################################
23!Bondi Module
24
25MODULE Problem
26
27 USE DataDeclarations
28 USE GlobalDeclarations
29 USE PhysicsDeclarations
30 USE ParticleDeclarations
31 USE CommonFunctions
32 USE Bondi
33 IMPLICIT NONE
34 PRIVATE
35
36 PUBLIC ProblemModuleInit, ProblemGridInit, &
37 ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
38 TYPE(PointGravityDef), POINTER :: PointGravityObj
39 REAL(KIND=qprec) :: namb, tamb, ibs, obs, mcent, r_bh, ramb, mdot_bondi, camb, r_sonic!, mdot_ruffert
40CONTAINS
41
42 SUBROUTINE ProblemModuleInit
43 INTEGER :: iErr
44 TYPE(ParticleDef), POINTER :: Particle
45 NAMELIST /ProblemData/ namb, tamb, ibs, obs, mcent
46
47
48 OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data')
49 READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
50 CLOSE(PROBLEM_DATA_HANDLE, IOSTAT=iErr)
51
52 IF (.not. lRestart) THEN
53 NULLIFY(Particle)
54 CALL CreateParticle(Particle)
55 Particle%q(1)=mcent*MSolar/rScale/lScale**3
56 Particle%xloc=0
57 Particle%iAccrete=2 !0=no accretion, 1=fed, 2=!KRUMHOLZ_ACCRETION
58 Particle%lFixed=.true.
59 particle%buffer(0:MaxLevel)=ceiling(ibs/levels(0:MaxLevel)%dx)
60 CALL AddSinkParticle(Particle)
61 CALL CreatePointGravityObject(Particle%PointGravityObj)
62 Particle%PointGravityObj%soft_length=1d0*levels(MaxLevel)%dx
63 Particle%PointGravityObj%soft_function=SPLINESOFT
64 Particle%PointGravityObj%Mass=mcent*MSolar/rScale/lScale**3 !Central object mass in computational units
65
66 END IF
67
68 !Calculate Bondi radius of initial particle mass
69 r_BH=ScaleGrav*Particle%q(1)/(gamma*tamb/TempScale)
70 IF (MPI_ID == 0) THEN
71 write(*,*) 'Bondi radius (CU)= ', r_BH
72 END IF
73
74 r_BH=(G*mcent*MSolar*xmu*amu)/(gamma*Boltzmann*tamb)
75 r_BH=r_BH/lscale !scaled into CU
76! IF (MPI_ID == 0) THEN
77! write(*,*) 'Bondi radius (CU)= ', r_BH
78! END IF
79
80 !Calculate theoretical accretion rate
81 ramb = (namb*amu*xmu)/rscale
82 camb=sqrt(gamma*tamb/tempscale)
83 mdot_bondi=(4*pi*bondi_lambda*ramb*(ScaleGrav*Particle%q(1))**2) / (camb**3)
84 !mdot_ruffert = (4*pi*ramb*r_bh**2)*(bondi_lambda**2*)
85 r_sonic = ((5.0-3.0*gamma)/4.0)*((ScaleGrav*particle%q(1))/camb**2)
86 IF (MPI_ID == 0) THEN
87 PRINT *, 'cs_inf=', camb
88 PRINT *, 'sonic radius=', r_sonic
89 PRINT *, 'T_inf=', Tamb/TempScale
90 PRINT *, "rho_inf=", ramb
91 PRINT *, "mdot_bondi (CU) = ", mdot_bondi
92 PRINT *, 'M_star=', particle%q(1)
93 END IF
94
95 END SUBROUTINE ProblemModuleInit
96
97 !> Initial Conditions
98 !! @param Info Info object
99 SUBROUTINE ProblemGridInit(Info)
100 !! @brief Initializes the grid data according to the requirements of the problem.
101 !! @param Info A grid structure.
102 TYPE (InfoDef) :: Info
103 INTEGER :: i,j,k
104 INTEGER :: rmbc,zrmbc,level
105 INTEGER :: mx, my, mz
106 INTEGER :: iErr
107 REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:) :: q
108 REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,r
109 REAL(KIND=qPREC) :: x_, y_, z_, rho, v, temp
110 level=Info%level
111 q=>Info%q
112 ! Calculating the number of ghost cells on each side of the grid.
113 rmbc=levels(level)%gmbc(levels(level)%step)
114 dx=levels(level)%dX
115 dy=dx
116 SELECT CASE(nDim)
117 CASE(2)
118 zrmbc=0
119 dz=0d0
120 CASE(3)
121 zrmbc=rmbc
122 dz=dx
123 END SELECT
124 mx = Info%mX(1)
125 my = Info%mX(2)
126 mz = Info%mX(3)
127 xl=Info%xBounds(1,1)
128 yl=Info%xBounds(2,1)
129 zl=Info%xBounds(3,1)
130
131 !All of the values below are in computational units
132 !print *, 'levels(0)%dx=', levels(0)%dx
133
134
135 !the was the following is calculated, need the original to be at 0,0,0:
136 DO i=1-rmbc, mx+rmbc
137 x = (xl+(REAL(i,xPrec)-half)*dx)
138 DO j=1-rmbc, my+rmbc
139 y = (yl+(REAL(j,xPrec)-half)*dy)
140 DO k=1-zrmbc, mz+zrmbc
141 z = (zl+(REAL(k,xPrec)-half)*dz)
142 r = sqrt(x**2 + y**2+z**2)
143 !IF (mpi_id==0) THEN
144 ! print *, 'level=', level
145 !END IF
146 ! Calculate non-dimensional radius(x_), density(z_), and velocity(y_)
147 x_=max(r, ibs-5d0*levels(0)%dx)/r_BH
148 !print *, 'r=', r
149 !print *, 'r/r_BH=', r/r_BH
150 !print *, 'x_=', x_
151 z_=BH_alpha(x_)
152 !print *, 'z_=', z_
153 y_=Bondi_lambda/(z_*x_**(myDim-1))
154 !print *, 'y_=', y_
155
156 ! Then calculate physical values
157 rho=z_*namb/nScale
158 v=-y_*sqrt(gamma*tamb/TempScale)
159 temp=tamb*(z_**(gamma1))/TempScale
160
161 q(i,j,k,1)=rho
162 q(i,j,k,ivx)=rho*v*x/r
163 q(i,j,k,ivy)=rho*v*y/r
164 !erica added the following for the z component:
165 q(i,j,k,ivz)=rho*v*z/r
166 q(i,j,k,iE)=half*rho*v**2+gamma7*rho*temp
167
168 END DO
169 END DO
170 END DO
171 END SUBROUTINE ProblemGridInit
172
173 !! @param Info Info object
174 SUBROUTINE ProblemBeforeStep(Info)
175 !! @brief Performs any tasks required before the advance step.
176 !! @param Info A grid structure.
177 TYPE (InfoDef) :: Info
178 INTEGER :: i,j,k
179 INTEGER :: rmbc,zrmbc,level
180 INTEGER :: mx, my, mz
181 REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:) :: q
182 REAL(KIND=xprec) :: x,y,z,xl,yl,zl,dx,dy,dz,r
183 REAL(KIND=qPREC) :: x_, y_, z_, rho, v, temp
184
185 level=Info%level
186 q=>Info%q
187 ! Calculating the number of ghost cells on each side of the grid.
188 rmbc=levels(level)%gmbc(levels(level)%step)
189 dx=levels(level)%dX
190 dy=dx
191 SELECT CASE(nDim)
192 CASE(2)
193 zrmbc=0
194 dz=0d0
195 CASE(3)
196 zrmbc=rmbc
197 dz=dx
198 END SELECT
199 mx = Info%mX(1)
200 my = Info%mX(2)
201 mz = Info%mX(3)
202 xl=Info%xBounds(1,1)
203 yl=Info%xBounds(2,1)
204 zl=Info%xBounds(3,1)
205
206 !All of the values below are in computational units
207
208 DO i = 1-rmbc, mx+rmbc
209 x = (xl+(REAL(i,xPrec)-half)*dx)
210 DO j = 1-rmbc, my+rmbc
211 y = (yl+(REAL(j,xPrec)-half)*dy)
212 DO k = 1-zrmbc,mz+zrmbc
213 z = (zl+(REAL(k,xPrec)-half)*dz)
214
215 r = sqrt(x**2 + y**2 +z**2)
216! r = sqrt(x**2 + y**2) -- this was original line, changed to the previous
217
218 IF (r < ibs .OR. r > obs) THEN !Set cells to analytic solution
219
220 x_=max(ibs-5d0*levels(0)%dx, r)/r_BH !Adjust values deep inside inner region to avoid extremely high velocities etc...
221 z_=BH_alpha(x_)
222 y_=Bondi_lambda/(z_*x_**(myDim-1))
223 rho=z_*namb/nScale
224 v=-y_*sqrt(gamma*tamb/TempScale)
225 temp=tamb*(z_**(gamma1))/TempScale
226
227 q(i,j,k,1)=rho
228 q(i,j,k,ivx)=rho*v*x/r
229 q(i,j,k,ivy)=rho*v*y/r
230 !erica added z component:
231 q(i,j,k,ivz)=rho*v*z/r
232 q(i,j,k,iE)=half*rho*v**2+gamma7*rho*temp
233
234 END IF
235
236
237
238 END DO
239 END DO
240 END DO
241
242 END SUBROUTINE ProblemBeforeStep
243
244 !> Does nothing
245 !! @param Info Info object
246 SUBROUTINE ProblemAfterStep(Info)
247 !! @brief Performs any post-step corrections that are required.
248 !! @param Info A grid structure.
249 TYPE (InfoDef) :: Info
250 CALL ProblemBeforeStep(Info)
251 END SUBROUTINE ProblemAfterStep
252
253 !> Does nothing
254 !! @param Info Info object
255 SUBROUTINE ProblemSetErrFlag(Info)
256 !! @brief Sets error flags according to problem-specific conditions..
257 !! @param Info A grid structure.
258 TYPE (InfoDef) :: Info
259 END SUBROUTINE ProblemSetErrFlag
260
261 SUBROUTINE ProblemBeforeGlobalStep(n)
262 INTEGER :: n
263 END SUBROUTINE ProblemBeforeGlobalStep
264
265END MODULE Problem
266
267