bondi2d.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: bondi2d.f90 #
5!# #
6!# Copyright (C) 2006-2018 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
9!# #
10!# This program is free software; you can redistribute it and/or modify #
11!# it under the terms of the GNU General Public License as published by #
12!# the Free Software Foundation; either version 2 of the License, or (at #
13!# your option) any later version. #
14!# #
15!# This program is distributed in the hope that it will be useful, but #
16!# WITHOUT ANY WARRANTY; without even the implied warranty of #
17!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18!# NON INFRINGEMENT. See the GNU General Public License for more #
19!# details. #
20!# #
21!# You should have received a copy of the GNU General Public License #
22!# along with this program; if not, write to the Free Software #
23!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24!# #
25!#############################################################################
26
27!----------------------------------------------------------------------------!
43!----------------------------------------------------------------------------!
44PROGRAM bondi2d
45 USE fosite_mod
46#include "tap.h"
47 IMPLICIT NONE
48 INTERFACE
49 PURE SUBROUTINE funcd(x,fx,dfx,plist)
50 IMPLICIT NONE
51 REAL, INTENT(IN) :: x
52 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
53 REAL, INTENT(OUT) :: fx,dfx
54 END SUBROUTINE funcd
55 END INTERFACE
56 !--------------------------------------------------------------------------!
57 ! general constants
58 REAL, PARAMETER :: msun = 1.989e+30 ! solar mass [kg]
59 REAL, PARAMETER :: gn = 6.67408e-11 ! gravitional constant (SI)
60 ! simulation parameters
61 REAL, PARAMETER :: tsim = 20.0 ! simulation time [TAU] (free fall)
62 REAL, PARAMETER :: accmass = 1.0*msun ! mass of the accreting object
63 REAL, PARAMETER :: gamma = 1.4 ! ratio of specific heats
64 ! boundary conditions
65 REAL, PARAMETER :: rhoinf = 1.0e-20 ! density at infinity [kg/m^3]
66 REAL, PARAMETER :: csinf = 1.0e+04 ! sound speed at infinity [m/s]
67 ! mesh settings
68 INTEGER, PARAMETER :: mgeo = cylindrical ! geometry
69 INTEGER, PARAMETER :: xres = 100 ! x-resolution
70 INTEGER, PARAMETER :: yres = 6 ! y-resolution
71 INTEGER, PARAMETER :: zres = 1 ! z-resolution
72 REAL, PARAMETER :: rin = 0.1 ! inner/outer radii in terms of
73 REAL, PARAMETER :: rout = 2.0 ! the Bondi radius RB, ROUT > 1
74 REAL, PARAMETER :: gpar = 1.0 ! geometry scaling parameter in [RB]
75 ! output parameters
76 INTEGER, PARAMETER :: onum = 10 ! number of output data sets
77 CHARACTER(LEN=256), PARAMETER & ! output data dir
78 :: odir = './'
79 CHARACTER(LEN=256), PARAMETER & ! output data file name
80 :: ofname = 'bondi2d'
81 ! some derives quandities
82 REAL :: rb ! Bondi radius
83 REAL :: tau ! free fall time scale
84 !--------------------------------------------------------------------------!
85 CLASS(fosite), ALLOCATABLE :: sim
86 REAL, DIMENSION(:), ALLOCATABLE :: sigma
87 REAL :: sum_numer, sum_denom
88 INTEGER :: n,den,vel,pre
89 LOGICAL :: ok
90 !--------------------------------------------------------------------------!
91
92 ALLOCATE(sim)
93 CALL sim%InitFosite()
94
95#ifdef PARALLEL
96 IF (sim%GetRank().EQ.0) THEN
97#endif
98tap_plan(4)
99#ifdef PARALLEL
100 END IF
101#endif
102
103 CALL makeconfig(sim, sim%config)
104! CALL PrintDict(config)
105 CALL sim%Setup()
106 ! set initial condition
107 CALL initdata(sim%Mesh, sim%Physics, sim%Fluxes, sim%Timedisc)
108 ! run the simulation
109 CALL sim%Run()
110 ok = .NOT.sim%aborted
111 ! compare with exact solution if requested
112 IF (ASSOCIATED(sim%Timedisc%solution)) THEN
113 ALLOCATE(sigma(sim%Physics%VNUM))
114 DO n=1,sim%Physics%VNUM
115 ! use L1 norm to estimate the deviation from the exact solution:
116 ! Σ |pvar - pvar_exact| / Σ |pvar_exact|
117 sum_numer = sum(abs(sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
118 sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n) &
119 -sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
120 sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
121 sum_denom = sum(abs(sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
122 sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
123#ifdef PARALLEL
124 IF (sim%GetRank().GT.0) THEN
125 CALL mpi_reduce(sum_numer,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
126 CALL mpi_reduce(sum_denom,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
127 ELSE
128 CALL mpi_reduce(mpi_in_place,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
129 CALL mpi_reduce(mpi_in_place,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
130#endif
131 sigma(n) = sum_numer / sum_denom
132#ifdef PARALLEL
133 END IF
134#endif
135 END DO
136 ELSE
137 sigma(:) = 0.0
138 END IF
139
140 den = sim%Physics%DENSITY
141 vel = sim%Physics%XVELOCITY
142 pre = sim%Physics%PRESSURE
143
144#ifdef PARALLEL
145 IF (sim%GetRank().EQ.0) THEN
146#endif
147tap_check(ok,"stoptime reached")
148! These lines are very long if expanded. So we can't indent it or it will be cropped.
149tap_check_small(sigma(den),1.0e-02,"density deviation < 1%")
150tap_check_small(sigma(vel),1.0e-02,"radial velocity deviation < 1%")
151! skip azimuthal velocity deviation, because exact value is 0
152tap_check_small(sigma(pre),1.0e-02,"pressure deviation < 1%")
153tap_done
154#ifdef PARALLEL
155 END IF
156#endif
157
158 IF (ALLOCATED(sigma)) DEALLOCATE(sigma)
159 CALL sim%Finalize()
160 DEALLOCATE(sim)
161
162CONTAINS
163
164 SUBROUTINE makeconfig(Sim, config)
165 IMPLICIT NONE
166 !------------------------------------------------------------------------!
167 CLASS(fosite), INTENT(INOUT) :: Sim
168 TYPE(dict_typ), POINTER :: config
169 !------------------------------------------------------------------------!
170 ! Local variable declaration
171 INTEGER :: bc(6)
172 TYPE(dict_typ), POINTER :: mesh, physics, boundary, datafile, &
173 grav, pmass, timedisc, fluxes
174 REAL :: x1,x2,y1,y2,z1,z2
175 !------------------------------------------------------------------------!
176 ! derived constants
177 rb = gn * accmass / csinf**2 ! bondi radius [m] !
178 tau = rb / csinf ! free fall time scale [s] !
179
180 ! mesh settings
181 SELECT CASE(mgeo)
182 CASE(cylindrical)
183 x1 = rin * rb
184 x2 = rout * rb
185 y1 = 0.0
186 y2 = 2*pi
187 z1 = 0.0
188 z2 = 0.0
189 bc(west) = no_gradients
190 bc(east) = custom
191 bc(south) = periodic
192 bc(north) = periodic
193 bc(bottom) = no_gradients
194 bc(top) = no_gradients
195 CASE DEFAULT
196 CALL sim%Physics%Error("bondi2d::MakeConfig","geometry currently not supported")
197 END SELECT
198
199 mesh => dict( &
200 "meshtype" / midpoint, &
201 "geometry" / mgeo, &
202 "inum" / xres, &
203 "jnum" / yres, &
204 "knum" / zres, &
205 "xmin" / x1, &
206 "xmax" / x2, &
207 "ymin" / y1, &
208 "ymax" / y2, &
209 "zmin" / z1, &
210 "zmax" / z2, &
211 "gparam" / (gpar*rb))
212
213 ! boundary conditions
214 boundary => dict( &
215 "western" / bc(west), &
216 "eastern" / bc(east), &
217 "southern" / bc(south), &
218 "northern" / bc(north), &
219 "bottomer" / bc(bottom), &
220 "topper" / bc(top))
221
222 ! physics settings
223 physics => dict( &
224 "problem" / euler, &
225 "gamma" / gamma)
226
227 ! flux calculation and reconstruction method
228 fluxes => dict( &
229 "order" / linear, &
230 "fluxtype" / kt, &
231 "variables" / conservative, &
232 "limiter" / monocent, &
233 "theta" / 1.2)
234
235 ! gravity term due to a point mass
236 pmass => dict( &
237 "gtype" / pointmass, &
238 "mass" / accmass, &
239 "outbound" / 0)
240
241 ! source term due to all gravity terms
242 grav => dict( &
243 "stype" / gravity, &
244 "pmass" / pmass)
245
246 ! time discretization settings
247 timedisc => dict( &
248 "method" / modified_euler, &
249 "order" / 3, &
250 "cfl" / 0.4, &
251 "stoptime" / (tsim * tau), &
252 "dtlimit" / (1.0e-6 * tau), &
253 "output/solution" / 1, &
254 "maxiter" / 1000000)
255
256 ! initialize data input/output
257 datafile => dict( &
258 "fileformat" / vtk, &
259 "filename" / (trim(odir) // trim(ofname)), &
260 "count" / onum)
261
262 config => dict( &
263 "mesh" / mesh, &
264 "physics" / physics, &
265 "boundary" / boundary, &
266 "fluxes" / fluxes, &
267 "sources/grav" / grav, &
268 "timedisc" / timedisc, &
269 "datafile" / datafile)
270 END SUBROUTINE makeconfig
271
272
273 SUBROUTINE initdata(Mesh,Physics,Fluxes,Timedisc)
274 IMPLICIT NONE
275 !------------------------------------------------------------------------!
276 CLASS(physics_base), INTENT(IN) :: Physics
277 CLASS(mesh_base), INTENT(IN) :: Mesh
278 CLASS(fluxes_base), INTENT(IN) :: Fluxes
279 CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
280 !------------------------------------------------------------------------!
281 ! Local variable declaration
282 REAL :: r,rho,vr,cs2
283 INTEGER :: i,j,k,l
284 CHARACTER(LEN=64) :: info_str
285 !------------------------------------------------------------------------!
286 ! initial condition
287 SELECT TYPE(pvar => timedisc%pvar)
288 TYPE IS(statevector_euler)
289 ! constant density and pressure, vanishing velocity
290 pvar%density%data1d(:) = rhoinf
291 pvar%velocity%data1d(:) = 0.
292 pvar%pressure%data1d(:) = rhoinf * csinf**2 / gamma
293 CLASS DEFAULT
294 CALL physics%Error("bondi2d::InitData","only non-isothermal HD supported")
295 END SELECT
296
297 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
298
299 ! store exact stationary solution if requested,
300 ! i.e. output/solution == 1 in timedisc config
301 IF (ASSOCIATED(timedisc%solution)) THEN
302 ! compute the stationary 2D planar bondi solution
303 SELECT TYPE(pvar => timedisc%solution)
304 TYPE IS(statevector_euler)
305 DO k=mesh%KMIN,mesh%KMAX
306 DO j=mesh%JMIN,mesh%JMAX
307 DO i=mesh%IMIN,mesh%IMAX
308 r = mesh%radius%bcenter(i,j,k)
309 CALL bondi(r/rb,gamma,rhoinf,csinf,rho,vr)
310 cs2 = csinf**2 * (rho/rhoinf)**(gamma-1.0)
311 pvar%density%data3d(i,j,k) = rho
312 ! set the minimum of the exact velocity to some value > 0
313 ! for comparison with the numerical results
314!NEC$ UNROLL(3)
315 DO l=1,physics%VDIM
316 pvar%velocity%data4d(i,j,k,l) = sign(&
317 max(epsilon(cs2),abs(vr*mesh%posvec%bcenter(i,j,k,l)/r)), &
318 vr*mesh%posvec%bcenter(i,j,k,l))
319 END DO
320 pvar%pressure%data3d(i,j,k) = rho * cs2 / gamma
321 END DO
322 END DO
323 END DO
324 END SELECT
325 END IF
326
327 ! boundary condition: subsonic inflow according to Bondi's solution
328 ! calculate Bondi solution for y=ymin..ymax at xmax
329 SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
330 CLASS IS (boundary_custom)
331 ! this tells the boundary routine the actual boundary condition for each variable
332 CALL beast%SetCustomBoundaries(mesh,physics, &
333 (/custom_fixed,custom_nograd,custom_fixed,custom_fixed/))
334 DO k=mesh%KMIN,mesh%KMAX
335 DO j=mesh%JMIN,mesh%JMAX
336 DO i=1,mesh%GINUM
337 ! get distance to the origin for each boundary cell
338 r = mesh%radius%bcenter(mesh%IMAX+i,j,k)
339 CALL bondi(r/rb,gamma,rhoinf,csinf,rho,vr)
340 cs2 = csinf**2 * (rho/rhoinf)**(gamma-1.0)
341 beast%data(i,j,k,physics%DENSITY) = rho
342 DO l=1,physics%VDIM
343 beast%data(i,j,k,physics%XVELOCITY+l-1) &
344 = vr*mesh%posvec%bcenter(mesh%IMAX+i,j,k,l)/r
345 END DO
346 beast%data(i,j,k,physics%PRESSURE) = rho * cs2 / gamma
347 END DO
348 END DO
349 END DO
350 END SELECT
351
352 ! print some information
353 CALL mesh%Info(" DATA-----> initial condition: 2D Bondi accretion")
354 WRITE(info_str,"(ES11.3)") rb
355 CALL mesh%Info(" " // "Bondi radius: " &
356 // trim(info_str) // " m")
357 WRITE(info_str,"(ES11.3)") tau
358 CALL mesh%Info( " " // "Free fall time: " &
359 // trim(info_str) // " s")
360 END SUBROUTINE initdata
361
362 SUBROUTINE bondi(r,gamma,rhoinf,csinf,rho,vr)
363 USE roots
364 IMPLICIT NONE
365 !------------------------------------------------------------------------!
366 ! computes the Bondi solution for rotationally symmetric accretion in !
367 ! planar geometry !
368 ! uses funcd, GetRoot_newton !
369 ! INPUT paramter: !
370 ! r : radius in units of the Bondi radius r_b = G*M/csinf**2 !
371 ! gamma : ratio of specific heats (1 < gamma < 5/3) !
372 ! rhoinf: density at infinity !
373 ! csinf : speed of sound at infinity !
374 ! OUTPUT paramter: !
375 ! rho : density @ r !
376 ! vr : radial velocity @ r !
377 !------------------------------------------------------------------------!
378 REAL, INTENT(IN) :: r,gamma,rhoinf,csinf
379 REAL, INTENT(OUT) :: rho,vr
380 !------------------------------------------------------------------------!
381 REAL, PARAMETER :: xacc = 1.0e-6 ! accuracy for root finding
382 REAL :: gp1,gm1,rc,chi,lambda,psi,gr
383 REAL :: params(2)
384 INTEGER :: err
385 !------------------------------------------------------------------------!
386 ! for convenience
387 gm1 = gamma - 1.0
388 gp1 = gamma + 1.0
389
390 ! critical radius
391 rc = (3.0-gamma) / 2.0
392
393 ! critical dimensionless accretion rate
394 lambda = rc**(-rc/gm1)
395
396 ! Newton-Raphson to solve Bondis equations for psi
397 chi = r / lambda
398 gr = chi**(2.*gm1/gp1) * (1./gm1 + 1./r)
399 params(1) = gm1
400 params(2) = gr
401 IF (r.LT.rc) THEN
402 CALL getroot_newton(funcd,1.0,gr,psi,err,plist=params)
403 ELSE
404 CALL getroot_newton(funcd,1.0e-5,1.0,psi,err,plist=params)
405 END IF
406
407 ! return values
408 rho = rhoinf * chi**(-2./gp1) / psi ! density
409 vr = -csinf * chi**(-gm1/gp1) * psi ! radial velocity
410 END SUBROUTINE bondi
411
412END PROGRAM bondi2d
413
414! find the root of fy to compute the exact Bondi solution
415PURE SUBROUTINE funcd(y,fy,dfy,plist)
416 IMPLICIT NONE
417 !------------------------------------------------------------------------!
418 REAL, INTENT(IN) :: y
419 REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
420 REAL, INTENT(OUT) :: fy,dfy
421 !------------------------------------------------------------------------!
422 ! plist(1) = gamma - 1
423 ! plist(2) = gr
424 !------------------------------------------------------------------------!
425 fy = 0.5*y*y + y**(-plist(1)) / plist(1) - plist(2)
426 dfy = y - y**(-plist(1)-1.)
427END SUBROUTINE funcd
428
429
430
431
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
pure subroutine funcd(y, fy, dfy, plist)
Definition: bondi2d.f90:416
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
program bondi2d
Program and data initialization for 2D Bondi accretion.
Definition: bondi2d.f90:44
subroutine bondi(r, gamma, rhoinf, csinf, rho, vr)
Definition: bondi2d.f90:363
root finding subroutines
Definition: roots.f90:45
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)
Definition: roots.f90:247
main fosite class
Definition: fosite.f90:71