pringle3d.f90

Jim Pringles famous viscous spreading ring solution in 3D

Author
Björn Sperling
Tobias Illenseer

This variant simulates the ring with vertical extent with and without rotational symmetry. The equations are solved in non-dimensional units, using the Keplerian velocity at the location of the initial ring at R=1 as the velocity scale.

Remarks
(1) This test is very sensitive to mesh settings. If there are not enough mesh points at small radii instabilities grow starting at the inner boundary.
(2) There is an unsolved stability issue with the viscous source term. The stability condition doesn't seem to be sufficient. Thus in case of high Reynolds numbers one has to adjust the viscous CFL number cvis to enforce smaller time steps and preserve stability.

References:

  • [lynden1974] Lynden-Bell, D. & Pringle, J. E. Evolution of Viscous Disks and Origin of Nebular Variables, M.N.R.A.S vol. 168(3) (1974) pp. 603-637 http://adsabs.harvard.edu/abs/1974MNRAS.168..603L
  • [pringle1981] Pringle, J. E. Accretion discs in astrophysics Annual review of astronomy and astrophysics. vol. 19 (1981) pp. 137-162 DOI: 10.1146/annurev.aa.19.090181.001033
  • [speith2003] Speith, R. and Kley, W. Stability of the viscously spreading ring Astron. Astrophys. vol. 399 (2003), pp. 395-407 DOI: 10.1051/0004-6361:20021783
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: pringle3d.f90 #
5 !# #
6 !# Copyright (C) 2008-2019 #
7 !# Bjoern Sperling <sperling@astrophysik.uni-kiel.de> #
8 !# Tobias Illenseer <tillense@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 !----------------------------------------------------------------------------!
60 !----------------------------------------------------------------------------!
61 PROGRAM pringle3d
62  USE fosite_mod
63  USE common_dict
64  IMPLICIT NONE
65  !--------------------------------------------------------------------------!
66  ! simulation parameters
67  REAL, PARAMETER :: TSIM = 5.0e-3 ! simulation time [TVIS] see below
68  ! should be < 1 ~ viscous time scale
69  ! these are the basic parameters; for stable solutions one requires
70  ! MA >> 1 and RE > MA
71  REAL, PARAMETER :: RE = 1.0e+4 ! Reynolds number (at R=1)
72  REAL, PARAMETER :: MA = 1.0e+2 ! Mach number (at R=1)
73  REAL, PARAMETER :: MDISK = 1.0e-2 ! disk mass << 1 = central mass
74  ! lower limit for nitial density
75  REAL, PARAMETER :: RHOMIN = 1.0e-8 ! minimal initial density
76  ! viscosity prescription
77  INTEGER, PARAMETER :: VISTYPE = beta
78 ! INTEGER, PARAMETER :: VISTYPE = POWERLAW
79  REAL, PARAMETER :: PL_EXP = 0.0 ! exponent for power law viscosity
80  REAL, PARAMETER :: TVIS = 4./3.*re ! viscous time scale
81  REAL, PARAMETER :: TINIT = 1.0e-3 ! time for initial condition [TVIS]
82  ! mesh settings
83 ! INTEGER, PARAMETER :: MGEO = CYLINDRICAL
84 ! INTEGER, PARAMETER :: MGEO = LOGCYLINDRICAL
85  INTEGER, PARAMETER :: MGEO = spherical !!! ATTENTION: not applicable in 1D
86  INTEGER, PARAMETER :: XRES = 200 ! x-resolution
87  INTEGER, PARAMETER :: YRES = 16 ! y-resolution
88  INTEGER, PARAMETER :: ZRES = 6 ! z-resolution
89  REAL, PARAMETER :: RMIN = 0.1 ! min radius of comp. domain
90  REAL, PARAMETER :: RMAX = 2.0 ! max radius of comp. domain
91  REAL, PARAMETER :: GPAR = 1.0 ! geometry scaling parameter
92  ! output parameters
93  INTEGER, PARAMETER :: ONUM = 100 ! number of output data sets
94  CHARACTER(LEN=256), PARAMETER & ! output data dir
95  :: ODIR = './'
96  CHARACTER(LEN=256), PARAMETER & ! output data file name
97  :: OFNAME = 'pringle3d'
98  ! derived parameters
99  REAL, PARAMETER :: CSISO = 1./ma ! isothermal speed of sound
100  !--------------------------------------------------------------------------!
101  CLASS(fosite), ALLOCATABLE :: Sim
102  !--------------------------------------------------------------------------!
103 
104  ALLOCATE(sim)
105  CALL sim%InitFosite()
106  CALL makeconfig(sim, sim%config)
107  CALL sim%Setup()
108  CALL initdata(sim,sim%Mesh, sim%Physics, sim%Timedisc)
109  CALL sim%Run()
110  CALL sim%Finalize()
111  DEALLOCATE(sim)
112 
113 CONTAINS
114 
115  SUBROUTINE makeconfig(Sim, config)
116  USE functions, ONLY : asinh
117  IMPLICIT NONE
118  !------------------------------------------------------------------------!
119  CLASS(fosite) :: Sim
120  TYPE(Dict_TYP),POINTER :: config
121  !------------------------------------------------------------------------!
122  ! Local variable declaration
123  INTEGER :: bc(6)
124  TYPE(Dict_TYP), POINTER :: mesh, physics, boundary, datafile, &
125  grav, vis, timedisc, fluxes, sources
126  REAL :: x1,x2,y1,y2,z1,z2,cvis
127  !------------------------------------------------------------------------!
128  INTENT(INOUT) :: sim
129  !------------------------------------------------------------------------!
130  ! mesh settings and boundary conditions
131  SELECT CASE(mgeo)
132  CASE(cylindrical)
133  x1 = rmin
134  x2 = rmax
135  y1 = -pi
136  y2 = pi
137 ! z1 = 0.0
138 ! z2 = 0.0
139  z1 = -0.1
140  z2 = 0.1
141 ! bc(WEST) = AXIS
142  bc(west) = no_gradients
143  bc(east) = no_gradients
144 ! bc(WEST) = ABSORBING
145 ! bc(EAST) = ABSORBING
146 ! bc(WEST) = CUSTOM
147 ! bc(EAST) = CUSTOM
148  bc(south) = periodic
149  bc(north) = periodic
150 ! bc(BOTTOM) = ABSORBING
151 ! bc(TOP) = ABSORBING
152  bc(bottom) = reflecting
153  bc(top) = reflecting
154  !!! ATTENTION: empirical formula found for standard parameters
155  !!! seems to work for RMIN=0.05 and XRES <= 1000
156  cvis = (xres/1207.0)**1.9
157  CASE(logcylindrical)
158  x1 = log(rmin/gpar)
159  x2 = log(rmax/gpar)
160  y1 = -pi
161  y2 = pi
162  z1 = 0.0
163  z2 = 0.0
164  bc(west) = no_gradients
165  bc(east) = no_gradients
166  bc(south) = periodic
167  bc(north) = periodic
168  bc(bottom) = no_gradients
169  bc(top) = no_gradients
170  cvis = 0.1
171  CASE(spherical)
172  x1 = rmin
173  x2 = rmax
174  y1 = pi/2 - 0.02*pi
175  y2 = pi/2 + 0.02*pi
176  z1 = -pi
177  z2 = pi
178 ! bc(WEST) = ABSORBING
179  bc(east) = absorbing
180 ! bc(WEST) = NO_GRADIENTS
181 ! bc(EAST) = NO_GRADIENTS
182  bc(west) = custom
183 ! bc(EAST) = CUSTOM
184 ! bc(SOUTH) = NO_GRADIENTS
185 ! bc(NORTH) = NO_GRADIENTS
186  bc(south) = absorbing
187  bc(north) = absorbing
188  bc(bottom) = periodic
189  bc(top) = periodic
190  !!! ATTENTION: empirical formula found for standard parameters
191  !!! seems to work for RMIN=0.05 and XRES <= 1000
192  cvis = (xres/1207.0)**1.9
193  CASE DEFAULT
194  CALL sim%Error("pringle::MakeConfig","mesh geometry not supported for Pringle disk")
195  END SELECT
196 
197  ! mesh settings
198  mesh => dict( &
199  "meshtype" / midpoint, &
200  "geometry" / mgeo, &
201  "inum" / xres, &
202  "jnum" / yres, &
203  "knum" / zres, &
204  "xmin" / x1, &
205  "xmax" / x2, &
206  "ymin" / y1, &
207  "ymax" / y2, &
208  "zmin" / z1, &
209  "zmax" / z2, &
210 ! "use_fargo" / 1, &
211 ! "fargo" / 2, &
212  "gparam" / gpar)
213 
214  ! boundary conditions
215  boundary => dict( &
216  "western" / bc(west), &
217  "eastern" / bc(east), &
218  "southern" / bc(south), &
219  "northern" / bc(north), &
220  "bottomer" / bc(bottom), &
221  "topper" / bc(top))
222 
223 
224  ! physics settings
225  physics => dict( &
226  "problem" / euler_isotherm, &
227  "units" / geometrical, &
228  "cs" / csiso)
229 
230  ! flux calculation and reconstruction method
231  fluxes => dict( &
232  "order" / linear, &
233  "fluxtype" / kt, &
234  "variables" / primitive, &
235  "limiter" / vanleer, &
236  "theta" / 1.2)
237 
238  ! source term due to viscosity
239  vis => dict( &
240  "stype" / viscosity, &
241  "vismodel" / vistype, &
242  "dynconst" / (1./re), &
243  "exponent" / pl_exp, &
244  "cvis" / cvis)
245 
246  ! source term due to a point mass
247  grav => dict( &
248  "stype" / gravity, &
249  "pmass/gtype"/ pointmass, & ! grav. accel. of a point mass
250  "pmass/mass" / 1.0, & ! non-dim. mass of the accreting object
251  "pmass/outbound" / 0) ! disable accretion
252 
253  ! combine all source terms
254  sources => dict( &
255  "viscosity" / vis, &
256  "grav" / grav)
257 
258  ! time discretization settings
259  timedisc => dict( &
260  "method" / ssprk, &
261  "stoptime" / (tsim * tvis), &
262  "cfl" / 0.4, &
263  "dtlimit" / (1e-10 * tvis), &
264 ! "rhstype" / 1, &
265  "maxiter" / 100000000, &
266  "tol_rel" / 0.01, &
267  "tol_abs" / (/1e-40,1e-5,1e-8,1e-5/))
268 
269 ! ! enable output of analytical solution
270 ! #ifdef HAVE_FGSL
271 ! IF (MGEO.EQ.CYLINDRICAL) &
272 ! CALL SetAttr(timedisc,"output/solution",1)
273 ! #endif
274 
275  ! initialize data input/output
276  datafile => dict( &
277  "fileformat"/ vtk, &
278  "filename" / (trim(odir) // trim(ofname)), &
279  "count" / onum)
280 
281  config => dict( &
282  "mesh" / mesh, &
283  "physics" / physics, &
284  "boundary" / boundary, &
285  "fluxes" / fluxes, &
286  "sources" / sources, &
287  "timedisc" / timedisc, &
288  "datafile" / datafile)
289  END SUBROUTINE makeconfig
290 
291  SUBROUTINE initdata(Sim,Mesh,Physics,Timedisc)
292  USE geometry_generic_mod
293  IMPLICIT NONE
294  !------------------------------------------------------------------------!
295  CLASS(fosite), INTENT(INOUT) :: Sim
296  CLASS(physics_base), INTENT(IN) :: Physics
297  CLASS(mesh_base), INTENT(IN) :: Mesh
298  CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
299  !------------------------------------------------------------------------!
300  ! Local variable declaration
301  CLASS(sources_base), POINTER :: sp
302  CLASS(sources_gravity), POINTER :: gp
303  CLASS(geometry_base), ALLOCATABLE :: geo_cyl
304  TYPE(Dict_TYP), POINTER :: geo_config
305  INTEGER :: i,j,k
306 #ifdef PARALLEL
307  INTEGER :: ierror
308 #endif
309  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3) &
310  :: bccyl
311  REAL :: r,z,vr,vphi,sden,mu,height,mass
312  CHARACTER(LEN=64) :: value
313  !------------------------------------------------------------------------!
314  geo_config => dict("geometry" / cylindrical)
315  CALL new_geometry(geo_cyl,geo_config)
316  CALL deletedict(geo_config)
317 
318  ! compute cylindrical coordinates for each cell bary center
319  CALL geo_cyl%Convert2Curvilinear(mesh%bccart,bccyl)
320 
321  ! custom boundary conditions if requested
322  SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
323  CLASS IS (boundary_custom)
324  CALL bwest%SetCustomBoundaries(mesh,physics, &
325  (/custom_nograd,custom_outflow,custom_kepler,custom_nograd/))
326  END SELECT
327  SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
328  CLASS IS (boundary_custom)
329  CALL beast%SetCustomBoundaries(mesh,physics, &
330  (/custom_nograd,custom_outflow,custom_kepler,custom_nograd/))
331  END SELECT
332 
333  ! get gravitational acceleration
334  sp => sim%Sources
335  DO
336  IF (ASSOCIATED(sp).EQV..false.) RETURN
337  SELECT TYPE(sp)
338  CLASS IS(sources_gravity)
339  gp => sp
340  EXIT
341  END SELECT
342  sp => sp%next
343  END DO
344  IF (.NOT.ASSOCIATED(sp)) CALL physics%Error("pringle::InitData","no gravity term initialized")
345 
346  SELECT CASE(vistype)
347  CASE(beta)
348  mu=1.0
349  CASE(powerlaw)
350  mu=pl_exp
351  CASE DEFAULT
352  CALL timedisc%Error("pringle::InitData","only pringle and beta viscosity possible")
353  END SELECT
354 
355  ! initial condition
356  SELECT TYPE(pvar => timedisc%pvar)
357  TYPE IS(statevector_eulerisotherm)
358  DO k=mesh%KGMIN,mesh%KGMAX
359  DO j=mesh%JGMIN,mesh%JGMAX
360  DO i=mesh%IMIN,mesh%IGMAX
361  ! distance to axis (r-coordinate in cylindrical coordinates)
362  r = bccyl(i,j,k,1)
363  ! z-coordinate
364  z = bccyl(i,j,k,3)
365  ! Keplerian velocity
366  vphi = r/(r**2+z**2)**0.75
367  ! compute pressure scale height
368  height = csiso * sqrt(r**3) ! h = cs / Omega = cs / SQRT(GM/r**3), with GM=1
369 #ifdef HAVE_FGSL
370  ! use exact solution at time TINIT
371  CALL viscousring_analytic(mu,tinit,r,sden,vr)
372 #else
373  ! use approximation for TINIT << 1
374  ! no need for modified Bessel functions
375  CALL viscousring_approx(mu,tinit,r,sden,vr)
376 #endif
377  pvar%density%data3d(i,j,k) = rhomin + sden / (sqrt(2*pi) * height) &
378  * exp(-0.5*(z/height)**2)
379 ! pvar%velocity%data4d(i,j,k,1) = vr
380 ! pvar%velocity%data4d(i,j,k,2) = vphi
381 ! pvar%velocity%data4d(i,j,k,3) = 0.0
382  END DO
383  END DO
384  END DO
385  ! determine disk mass
386  mass = sum(mesh%volume%data1d(:)*pvar%density%data1d(:),mask=mesh%without_ghost_zones%mask1d)
387 #ifdef PARALLEL
388  CALL mpi_allreduce(mpi_in_place,mass,1,default_mpi_real,mpi_sum, &
389  mesh%comm_cart,ierror)
390 #endif
391  ! rescale disk mass
392  pvar%density%data1d(:) = pvar%density%data1d(:) * mdisk / mass
393  pvar%velocity%data1d(:) = 0.0 ! initialize all velocities with zero
394  CALL gp%UpdateGravity(mesh,sim%Physics,sim%Fluxes,pvar,0.0,0.0)
395  pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
396  timedisc%GetCentrifugalVelocity(mesh,sim%Physics,sim%Fluxes,sim%Sources,(/0.,0.,1./), &
397  gp%accel%data4d)
398  CLASS DEFAULT
399  CALL timedisc%Error("pringle::InitData","only isothermal physics possible")
400  END SELECT
401 
402  IF (ASSOCIATED(timedisc%solution)) THEN
403  ! store initial condition in exact solution array
404  timedisc%solution = timedisc%pvar
405  END IF
406 
407  ! transform to conservative variables
408  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
409 
410  CALL mesh%Info(" DATA-----> initial condition: " // "pringle disk")
411  WRITE(value,"(ES9.3)") tvis
412  CALL mesh%Info(" " // "viscous timescale: " //trim(value))
413  WRITE(value,"(ES9.3)") re
414  CALL mesh%Info(" " // "Reynolds number: " //trim(value))
415  WRITE(value,"(ES9.3)") ma
416  CALL mesh%Info(" " // "Mach number: " //trim(value))
417  WRITE(value,"(ES9.3)") mdisk
418  CALL mesh%Info(" " // "Disk mass: " //trim(value))
419 
420  CALL geo_cyl%Finalize()
421  DEALLOCATE(geo_cyl)
422  END SUBROUTINE initdata
423 
424 
425  SUBROUTINE viscousring_approx(mu,tau,r,Sigma,vr)
426  IMPLICIT NONE
427  !--------------------------------------------------------------------!
428  REAL, INTENT(IN) :: mu,tau,r
429  REAL, INTENT(OUT) :: Sigma,vr
430  !--------------------------------------------------------------------!
431  REAL :: lambda,x,r4l,invr4l
432  !--------------------------------------------------------------------!
433  lambda = 1./(4.0-mu)
434  r4l = r**(0.25/lambda)
435  x = 2*lambda**2 / tau * r**(0.25/lambda)
436  invr4l = 1./r4l
437  sigma = mdisk/sqrt((4*pi)**3 * tau) * r4l**(1.5-9*lambda) &
438  * exp(-lambda**2/tau*(1.0-r4l)**2)
439  vr = -1.5/re * ((x*tau)/(2*lambda*lambda))**(4*lambda-2.) &
440  * (1.0 - 0.5*x/lambda * (x*tau/(2*lambda*lambda) - 1.0))
441  END SUBROUTINE viscousring_approx
442 
443 #ifdef HAVE_FGSL
444  SUBROUTINE viscousring_analytic(mu,tau,r,Sigma,vr)
445  IMPLICIT NONE
446  !--------------------------------------------------------------------!
447  REAL, INTENT(IN) :: mu,tau,r
448  REAL, INTENT(OUT) :: Sigma,vr
449  !--------------------------------------------------------------------!
450  REAL :: lambda,r4l,x
451  !--------------------------------------------------------------------!
452  lambda = 1./(4.0-mu)
453  r4l = r**(0.25/lambda)
454  x = 2*lambda**2 / tau * r4l
455  sigma = mdisk/(4*pi) * (lambda/tau) * r4l**2 * r**(-2.25) &
456  * exp(-lambda**2/tau*(1.0+r4l**2)) &
457  * modified_bessel_inu(lambda,x)
458  vr = -1.5/re * ((x*tau)/(2*lambda*lambda))**(4*lambda-2.) &
459  * (1.0 - 0.5*x/lambda * (x*tau/(2*lambda*lambda) &
460  - modified_bessel_inu(1+lambda,x)/modified_bessel_inu(lambda,x)))
461  END SUBROUTINE viscousring_analytic
462 
463  FUNCTION modified_bessel_inu(nu,x) RESULT(Inu)
464  USE fgsl
465  IMPLICIT NONE
466  !--------------------------------------------------------------------!
467  REAL, INTENT(IN) :: nu,x
468  REAL :: Inu
469  !--------------------------------------------------------------------!
470  INTEGER :: n
471  !--------------------------------------------------------------------!
472  IF (aint(nu).EQ.nu) THEN
473  ! nu is an integer -> use modified bessel functions of integer order
474  n = int(nu)
475  inu = fgsl_sf_bessel_icn(n,x)
476  ELSE
477  ! nu is not an integer -> use representation of modified bessel functions
478  ! through hypergeometric functions 0F1
479  inu = (0.5*x)**nu / fgsl_sf_gamma(1.0+nu) * fgsl_sf_hyperg_0f1(1.0+nu,(0.5*x)**2)
480  END IF
481  END FUNCTION modified_bessel_inu
482 #endif
483 
484 END PROGRAM pringle3d
Generated by   doxygen 1.8.14