mestel.f90

self-gravitating accretion disk

Author
Manuel Jung
Björn Sperling
Tobias Illenseer

2D simulation of a geometrically thin, self-gravitating accretion disk around a supermassive black hole in polar geometry with logarithmic radial spacing. The standard setup solves the non-isothermal inviscid Euler equations with thin-disk gray cooling (see sources_diskcooling_mod::lambda_gray ). Gravitational forces account for the central point mass as well as self-gravity of the disk.

The setup is based on those described in [britsch2006] .

Simulation parameters \( \quad \) --—
black hole mass \( 10^7\, \mathsf{M}_\odot \)
disk / black hole mass ratio \( 0.1 \)
mean molecular weight \( 6.02\cdot 10^{-4}\, \mathsf{kg/mol} \)
specific heat ratio \( 1.4 \)
inner radius \( 0.05\, \mathsf{pc} \)
outer radius \( 1\, \mathsf{pc} \)
Initial condition \( \quad \) --—
power law surface density, i. e. Mestel's disk \( \Sigma \propto 1/r + \mathsf{noise} \)
constant temperature \( 100\, \mathsf{K} \)
centrifugal balance \( v_\varphi^2 = -r \partial_r \Phi \)
column density
You can find a time-lapse movie showing the temporal evolution of the column density in the gallery.

References:

  • [britsch2006] M. Britsch, "Gravitational instability and fragmentation of self-gravitating accretion disks", PhD thesis, 2006
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: mestel.f90 #
5 !# #
6 !# Copyright (C) 2012-2019 #
7 !# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
8 !# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
10 !# #
11 !# This program is free software; you can redistribute it and/or modify #
12 !# it under the terms of the GNU General Public License as published by #
13 !# the Free Software Foundation; either version 2 of the License, or (at #
14 !# your option) any later version. #
15 !# #
16 !# This program is distributed in the hope that it will be useful, but #
17 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
18 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19 !# NON INFRINGEMENT. See the GNU General Public License for more #
20 !# details. #
21 !# #
22 !# You should have received a copy of the GNU General Public License #
23 !# along with this program; if not, write to the Free Software #
24 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25 !# #
26 !#############################################################################
27 
28 !----------------------------------------------------------------------------!
71 !----------------------------------------------------------------------------!
72 PROGRAM mestel
73  USE fosite_mod
74 #ifdef NECSXAURORA
75  USE asl_unified
76 #endif
77  IMPLICIT NONE
78  !--------------------------------------------------------------------------!
79  ! general constants
80  REAL, PARAMETER :: GN = 6.67384e-11 ! Newtons grav. constant [m^3/kg/s^2]
81  REAL, PARAMETER :: YEAR = 3.15576e+7 ! year [s]
82  REAL, PARAMETER :: PARSEC = 3.0857e+16 ! parsec [m]
83  REAL, PARAMETER :: AU = 1.49598e+11 ! astronomical unit [m]
84  REAL, PARAMETER :: MSUN = 1.989e+30 ! solar mass [kg]
85  REAL, PARAMETER :: RG = 8.31 ! gas constant
86  ! simulation parameters
87  REAL, PARAMETER :: SIMTIME = 1.0e+4*year ! simulation time [s]
88  REAL, PARAMETER :: MBH = 1.0e+7*msun ! initial black hole mass [kg]
89  REAL, PARAMETER :: MRATIO = 0.1 ! initial mdisk/mbh mass ratio
90  REAL, PARAMETER :: MDISK = mratio*mbh ! initial disk mass [kg]
91  REAL, PARAMETER :: TEMP = 100.0 ! initial temperature [K]
92  REAL, PARAMETER :: NOISE = 0.3 ! initial noise level
93  ! physics settings
94  INTEGER, PARAMETER :: PHYS = euler ! transport model
95  REAL, PARAMETER :: MU = 6.02e-4 ! mean molecular weight [kg/mol]
96  REAL, PARAMETER :: GAMMA = 1.4 ! ratio of specific heats
97  REAL, PARAMETER :: BETA_VIS = 1.0e-3 ! beta viscosity parameter
98  ! mesh settings
99  REAL, PARAMETER :: RMIN = 5.0e-2 * parsec ! inner radius [m]
100  REAL, PARAMETER :: RMAX = 1.e+0 * parsec ! outer radius [m]
101  REAL, PARAMETER :: RGEO = 1.0 * parsec ! geometry scaling constant
102  INTEGER, PARAMETER :: MGEO = logcylindrical ! mesh geometry
103  INTEGER, PARAMETER :: XRES = 64 ! mesh resolution (radial)
104  INTEGER, PARAMETER :: YRES = 128 ! mesh resolution (azimuthal)
105  INTEGER, PARAMETER :: ZRES = 1 ! mesh resolution (z)
106  ! output settings
107  INTEGER, PARAMETER :: ONUM = 100 ! number of output time steps
108  CHARACTER(LEN=256), PARAMETER :: &
109  OFNAME = 'markward', & ! data file name
110  odir = "./" ! output directory
111  !--------------------------------------------------------------------------!
112  CLASS(fosite), ALLOCATABLE :: Sim
113  !--------------------------------------------------------------------------!
114 
115 ALLOCATE(sim)
116 
117 #ifdef NECSXAURORA
118 CALL asl_library_initialize()
119 #endif
120 
121 CALL sim%InitFosite()
122 CALL makeconfig(sim%config)
123 CALL sim%Setup()
124 CALL initdata(sim%Timedisc,sim%Mesh,sim%Physics,sim%Fluxes,sim%Sources)
125 CALL sim%Run()
126 CALL sim%Finalize()
127 
128 #ifdef NECSXAURORA
129 CALL asl_library_finalize()
130 #endif
131 
132 DEALLOCATE(sim)
133 
134 CONTAINS
135  SUBROUTINE makeconfig(config)
136  IMPLICIT NONE
137  !------------------------------------------------------------------------!
138  TYPE(Dict_TYP),POINTER :: config
139  !------------------------------------------------------------------------!
140  ! Local variable declaration
141  TYPE(Dict_TYP),POINTER :: mesh,boundary,timedisc,datafile,&
142  sources,fluxes,grav,physics,&
143  vis,cooling
144  !------------------------------------------------------------------------!
145  physics => dict( &
146  "problem" / phys, &
147  "output/bccsound" / 1, &
148  "mu" / mu, &
149  "gamma" / gamma, &
150  "units" / si)
151 
152  fluxes => dict( &
153  "fluxtype" / kt, &
154  "order" / linear, &
155  "variables" / primitive, &
156  "limiter" / vanleer, &
157  "theta" / 1.2)
158 
159  mesh => dict( &
160  "meshtype" / midpoint, &
161  "geometry" / mgeo, &
162  "inum" / xres, &
163  "jnum" / yres, &
164  "knum" / zres, &
165  "xmin" / log(rmin/parsec), &
166  "xmax" / log(rmax/parsec), &
167  "ymin" / (-pi), &
168  "ymax" / ( pi), &
169  "zmin" / 0.0, &
170  "zmax" / 0.0, &
171  "gparam" / rgeo, &
172  "use_fargo" / 1, &
173  "fargo" / 1, &
174  "decomposition" / (/-1,1,1/), &
175  "output/volume" / 1 )
176 
177  boundary => dict( &
178  "western" / custom, &
179  "eastern" / custom, &
180  "southern" / periodic, &
181  "northern" / periodic, &
182  "bottomer" / reflecting, &
183  "topper" / reflecting)
184 
185  grav => dict( &
186  "stype" / gravity, &
187  "cvis" / 0.9, &
188  "energy" / 0, &
189  "output/height" / 1, &
190  "self/gtype" / spectral, &
191  "self/green" / 1, &
192  "pmass/gtype" / pointmass, &
193  "pmass/potential" / newton, &
194  "pmass/mass" / mbh, &
195  "pmass/outbound" / 1)
196 
197  ! cooling model
198  cooling => dict( &
199  "stype" / disk_cooling, &
200 ! "method" / GRAY, &
201  "method" / gammie, &
202  "b_cool" / 10.0, &
203  "Tmin" / 3.0, &
204  "rhomin" / 1.0e-30, &
205  "cvis" / 0.1)
206 
207  ! viscosity model
208  vis => dict( &
209  "stype" / viscosity, &
210  "vismodel" / beta, &
211  "dynconst" / beta_vis, &
212  "output/stress" / 1, &
213  "output/dynvis" / 1, &
214  "output/kinvis" / 1, &
215  "cvis" / 0.5)
216 
217  ! collect source terms
218  sources => dict( &
219  "diskcooling" / cooling, &
220 ! "viscosity" / vis, &
221  "grav" / grav)
222 
223  ! time discretization settings
224  timedisc => dict( &
225  "method" / ssprk, &
226  "tol_rel" / 1.0e-3, &
227  "cfl" / 0.3, &
228  "stoptime" / simtime, &
229  "dtlimit" / 1.0e-4, &
230  "tmin" / 3.0, &
231  "rhstype" / 1, &
232  "maxiter" / 2000000000)
233 
234  ! add absolute error bounds and output fields depending on physics
235  SELECT CASE(phys)
236  CASE(euler)
237  CALL setattr(timedisc, "tol_abs", (/1.0e-3, 1.0e-1, 1.0e-01, 1.0e+10/))
238  CALL setattr(timedisc, "output/xmomentum", 0)
239  CALL setattr(timedisc, "output/ymomentum", 0)
240  CALL setattr(timedisc, "output/energy", 0)
241  CALL setattr(timedisc, "output/rhs", 0)
242  CALL setattr(timedisc, "output/external_sources", 0)
243  CASE DEFAULT
244  CALL error(sim,"MakeConfig","Physics model not supported.")
245  END SELECT
246 
247  ! data file settings
248  datafile => dict( &
249  "fileformat" / vtk, &
250  "filename" / "mestel", &
251  "count" / onum)
252 
253  ! create configuration
254  config => dict( &
255  "physics" / physics, &
256  "fluxes" / fluxes, &
257  "mesh" / mesh, &
258  "boundary" / boundary, &
259  "sources" / sources, &
260  "timedisc" / timedisc, &
261  "datafile" / datafile)
262 
263  END SUBROUTINE makeconfig
264 
265 
266  SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes,Sources)
267  USE physics_euler_mod, ONLY : physics_euler, statevector_euler
268  IMPLICIT NONE
269  !------------------------------------------------------------------------!
270  CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
271  CLASS(mesh_base), INTENT(IN) :: Mesh
272  CLASS(physics_base), INTENT(INOUT) :: Physics
273  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
274  CLASS(sources_base), POINTER :: Sources
275  !------------------------------------------------------------------------!
276  ! Local variable declaration
277  CLASS(sources_base), POINTER :: sp
278  CLASS(sources_gravity), POINTER :: gp
279  CLASS(marray_base), ALLOCATABLE :: rands,Sigma
280 #ifdef PARALLEL
281  INTEGER :: ierror
282 #endif
283  REAL :: mass
284  CHARACTER(LEN=20) :: mdisk_str
285 #ifdef NECSXAURORA
286  INTEGER :: rng, n
287 #endif
288  !------------------------------------------------------------------------!
289  ALLOCATE(rands,sigma)
290  ! get gravitational acceleration
291  sp => sources
292  DO
293  IF (ASSOCIATED(sp).EQV..false.) RETURN
294  SELECT TYPE(sp)
295  CLASS IS(sources_gravity)
296  gp => sp
297  EXIT
298  END SELECT
299  sp => sp%next
300  END DO
301  IF (.NOT.ASSOCIATED(sp)) CALL physics%Error("mestel::InitData","no gravity term initialized")
302 
303  ! get random numbers for density noise
304  rands = marray_base()
305 #ifndef NECSXAURORA
306  CALL initrandseed(physics)
307  CALL random_number(rands%data1d)
308 #else
309  CALL asl_random_create(rng, asl_randommethod_mt19937_64)
310  CALL asl_random_distribute_uniform(rng)
311  n = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)
312  CALL asl_random_generate_d(rng, n, rands%data1d)
313  CALL asl_random_destroy(rng)
314 #endif
315  rands%data1d(:) = rands%data1d(:) * noise * 2.0 + (1.0 - noise)
316 
317  ! set surface density using radial power law (1/r) with a little noise
318  sigma = marray_base()
319  sigma%data1d(:) = rands%data1d(:) * rmin/mesh%radius%data2d(:,2)
320 
321  ! determine disk mass
322  mass = sum(mesh%volume%data1d(:)*sigma%data1d(:),mask=mesh%without_ghost_zones%mask1d(:))
323 #ifdef PARALLEL
324  CALL mpi_allreduce(mpi_in_place,mass,1,default_mpi_real,mpi_sum, &
325  mesh%comm_cart,ierror)
326 #endif
327 
328  ! setting for custom boundary conditions (western boundary)
329  SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
330  CLASS IS (boundary_custom)
331  CALL bwest%SetCustomBoundaries(mesh,physics, &
332  (/custom_nograd,custom_outflow,custom_kepler,custom_nograd/))
333  END SELECT
334 
335  ! setting for custom boundary conditions (eastern boundary)
336  SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
337  CLASS IS (boundary_custom)
338  CALL beast%SetCustomBoundaries(mesh,physics, &
339  (/custom_reflect,custom_reflect,custom_logexpol,custom_reflect/))
340  END SELECT
341 
342  SELECT TYPE (pvar => timedisc%pvar)
343  CLASS IS(statevector_euler)
344  ! 1. set surface density by rescaling Sigma
345  pvar%density%data1d(:) = sigma%data1d(:) * mdisk / mass
346  ! 2. initialize gravitational acceleration
347  CALL gp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
348  ! 3. set azimuthal velocity: balance initial radial gravitational
349  ! acceleration with centrifugal acceleration
350  pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
351  timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),gp%accel%data4d)
352  ! 4. transform velocities to rotating frame
353  ! ATTENTION: this works only if the second velocity is the azimuthal velocity
354  IF (mesh%OMEGA.GT.0.0) &
355  pvar%velocity%data2d(:,2) = pvar%velocity%data2d(:,2) - mesh%OMEGA*mesh%radius%data2d(:,2)
356  ! 5. set pressure using surface density and initial temperature
357  pvar%pressure%data1d(:) = physics%constants%RG/physics%MU * temp * pvar%density%data1d(:)
358  CLASS DEFAULT
359  CALL physics%Error("mestel::InitData","unsupported state vector")
360  END SELECT
361 
362  IF (mesh%FARGO.EQ.2) &
363  timedisc%w(:,:) = sqrt(physics%constants%GN*(mbh/mesh%radius%bcenter(:,mesh%JMIN,:)))
364 
365  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
366 
367  ! print some information
368  WRITE (mdisk_str, '(ES8.2)') mdisk/msun
369  CALL mesh%Info(" DATA-----> initial condition: Mestel's disk")
370  CALL mesh%Info(" disk mass: " // trim(mdisk_str) // " M_sun")
371 
372  CALL rands%Destroy()
373  CALL sigma%Destroy()
374  DEALLOCATE(rands,sigma)
375  END SUBROUTINE initdata
376 
377  SUBROUTINE initrandseed(Physics)
378  IMPLICIT NONE
379  !------------------------------------------------------------------------!
380  CLASS(physics_base), INTENT(IN) :: Physics
381  INTEGER :: i, n, clock
382  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
383  !------------------------------------------------------------------------!
384  ! Initialize random number generator with a seed based on the systems time
385  ! source: http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
386  CALL random_seed(size = n)
387  ALLOCATE(seed(n))
388  CALL system_clock(count=clock)
389  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
390 #ifdef PARALLEL
391  seed = seed + physics%GetRank()
392 #endif
393  CALL random_seed(put = seed)
394  DEALLOCATE(seed)
395  END SUBROUTINE initrandseed
396 END PROGRAM mestel
Generated by   doxygen 1.8.14