sgdisk.f90

Program and data initialization a self-gravitating disk (sgdisk)

Author
Jannes Klee
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: sgdisk.f90 #
5 !# #
6 !# Copyright (C) 2010-2019 #
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 !----------------------------------------------------------------------------!
33 !----------------------------------------------------------------------------!
34 PROGRAM sgdisk
35  USE fosite_mod
36 #ifdef NECSXAURORA
37  USE asl_unified
38 #endif
39 #ifdef HAVE_KROME
40  USE krome_main
41  USE krome_user
42 #endif
43 #ifdef PARALLEL
44 #ifdef HAVE_MPI_MOD
45  USE mpi
46 #endif
47 #endif
48  IMPLICIT NONE
49 #ifdef PARALLEL
50 #ifdef HAVE_MPIF_H
51  include 'mpif.h'
52 #endif
53 #endif
54 
55  !--------------------------------------------------------------------------!
56  ! general constants
57  REAL, PARAMETER :: GN = 6.6742d-8 ! grav. constant (cgs) !
58  REAL, PARAMETER :: RG = 8.31447d+7 ! molar gas constant (cgs) !
59  REAL, PARAMETER :: MSUN = 1.989d+33 ! solar mass [g] !
60  REAL, PARAMETER :: AU = 1.49597870691e+13 ! ast. unit [cm] !
61  REAL, PARAMETER :: YEAR = 3.15576e+7 ! Julian year [s] !
62  REAL, PARAMETER :: MSCALE = msun ! mass scaling param. !
63  ! simulation parameters
64  REAL, PARAMETER :: TSIM = 6.0 ! simulation time [ORP] !
65  REAL, PARAMETER :: MU = 2.35d+0 ! mean mol. mass [g/mol] !
66  REAL, PARAMETER :: VALPHA = 0.05 ! alpha visc. parameter !
67  ! disk & central object
68  REAL, PARAMETER :: MBH = 1.0*mscale ! central mass !
69  REAL, PARAMETER :: MDISK = 0.1*mscale ! initial disk mass !
70  REAL, PARAMETER :: T_INIT = 120.0 ! initial temperatur !
71  REAL, PARAMETER :: T_MIN = 10.0 ! minimal temperatur !
72  ! cooling parameters
73  REAL, PARAMETER :: BETA_C = 10.0 ! cooling parameter !
74  REAL, PARAMETER :: GAMMA = 1.6 ! adiabatic index !
75  ! mesh settings
76  INTEGER, PARAMETER :: MGEO = logcylindrical ! geometry of the mesh !
77  REAL, PARAMETER :: GPAR = au ! geom. scal. parameter !
78  REAL, PARAMETER :: RMIN = 1.0*gpar ! inner radius of the disk !
79  REAL, PARAMETER :: RMAX = 25.0*gpar ! outer radius of the grid !
80  INTEGER, PARAMETER :: XRES = 256 ! x-resolution !
81  INTEGER, PARAMETER :: YRES = xres*3 ! y-resolution !
82  INTEGER, PARAMETER :: ZRES = 1 ! z-resolution !
83  ! mestel
84  LOGICAL :: DAMP ! wave damping !
85  REAL, PARAMETER :: NOISE = 0.3 ! initial noise level !
86  ! output file parameter
87  INTEGER, PARAMETER :: ONUM = 100 ! number of outputs !
88  CHARACTER(LEN=256), PARAMETER &
89  :: ODIR = "./" ! output directory !
90  CHARACTER(LEN=256) :: OFNAME = 'sgdisk' ! output name !
91  !--------------------------------------------------------------------------!
92  CLASS(fosite), ALLOCATABLE :: Sim
93  REAL :: CSISO,OMEGA,ORP
94  !--------------------------------------------------------------------------!
95 
96 ALLOCATE(sim)
97 
98 #ifdef NECSXAURORA
99 CALL asl_library_initialize()
100 #endif
101 
102 CALL sim%InitFosite()
103 CALL makeconfig(sim, sim%config)
104 CALL sim%Setup()
105 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, sim%Fluxes, sim%Sources, &
106  sim%Timedisc%pvar, sim%Timedisc%cvar)
107 CALL sim%Run()
108 CALL sim%Finalize()
109 
110 #ifdef NECSXAURORA
111  CALL asl_library_finalize()
112 #endif
113 
114 DEALLOCATE(sim)
115 
116 CONTAINS
117 
118  SUBROUTINE makeconfig(Sim, config)
119  IMPLICIT NONE
120  !------------------------------------------------------------------------!
121  CLASS(fosite), INTENT(INOUT) :: Sim
122  TYPE(Dict_TYP), POINTER :: config
123  TYPE(Dict_TYP), POINTER :: mesh, physics, fluxes, boundary, grav,&
124  sources, timedisc, datafile, cooling, &
125  rotframe, pmass, self, damping
126  !------------------------------------------------------------------------!
127  ! some derived simulation parameters
128  csiso = sqrt(rg/mu*t_init)
129  omega = 0.0
130  orp = 2.*pi/(sqrt(gn*mbh/rmax**3.))
131  damp = .false.
132 
133  ! mesh settings
134  ! stellar orbits must be inside the central hole of the mesh
135  mesh => dict(&
136  "meshtype" / midpoint, &
137  "geometry" / mgeo, &
138  "inum" / xres, &
139  "jnum" / yres, &
140  "knum" / zres, &
141  "xmin" / log(rmin/gpar), &
142  "xmax" / log(rmax/gpar), &
143  "ymin" / 0.0, &
144  "ymax" / (2.0*pi), &
145  "zmin" / 0.0, &
146  "zmax" / 0.0, &
147  "use_fargo" / 1, &
148  "fargo" / 2, &
149  "decomposition" / (/ -1, 1, 1/), &
150 ! "output/radius" / 1, &
151 ! "output/dl" / 1, &
152  "output/position_vector" / 1, &
153  "gparam" / gpar)
154 
155 
156  ! physics settings
157  physics => dict(&
158  "problem" / euler, &
159  "mu" / mu, &
160  "gamma" / gamma, &
161  "units" / cgs)
162 ! "output/bccsound" / 1)
163 
164  ! boundary conditions
165  boundary => dict(&
166  "western" / custom,&
167  "eastern" / custom,&
168  "southern" / periodic, &
169  "northern" / periodic, &
170  "bottomer" / reflecting, &
171  "topper" / reflecting)
172 
173 
174  ! numerical fluxes and reconstruction method
175  fluxes => dict(&
176  "order" / linear, &
177  "fluxtype" / kt, &
178  "variables" / primitive, &
179  "passive_limiting" / .false., &
180  "limiter" / vanleer, &
181  "theta" / 1.2)
182 
183 
184  ! viscosity source term
185  rotframe => dict(&
186  "stype" / rotating_frame)
187 
188 
189  ! gravitational acceleration due to binary system
190  pmass => dict(&
191  "gtype" / pointmass, &
192  "mass" / mbh, &
193  "output/position" / 1)
194 
195  ! cooling function
196  cooling => dict(&
197  "stype" / disk_cooling, &
198  "output/Qcool" / 1, &
199  "rhomin" / 1.0e-30, &
200  "method" / gammie, &
201  "b_cool" / beta_c)
202 
203  ! self-gravity
204  self => dict(&
205  "gtype" / spectral, &
206  "output/potential" / 1, &
207  "self/green" / 1)
208 
209  ! collect all gravity-source terms
210  grav => dict(&
211  "stype" / gravity, &
212  "pmass" / pmass, &
213  "self" / self, &
214  "output/potential" / 1, &
215  "energy" / 0, &
216  "output/height" / 1, &
217  "output/accel" / 1)
218 
219  ! collect all source terms
220  sources => dict(&
221  "cooling" / cooling, &
222  "grav" / grav)
223 
224  ! add wave damping if requested
225  IF (damp) &
226  CALL setattr(sources,"damping",damping)
227 
228  ! time discretization settings
229  timedisc => dict(&
230  "method" / modified_euler, &
231  "cfl" / 0.3, &
232  "dtlimit" / 1.0e-50, &
233  "tol_rel" / 1.0e-3, &
234  "rhstype" / 1, &
235  "tol_abs" / (/ 1.0e-16, 1.0, 1.0e-16, 1.0 /), &
236  "output/energy" / 1, &
237  "output/rhs" / 1, &
238  "stoptime" / (orp*tsim), &
239  "pmin" / 1e-14, &
240  "tmin" / t_min, &
241  "checkdata" / ior(check_nothing,check_tmin), &
242  "maxiter" / 100000000)
243 
244  ! initialize data input/output
245  datafile => dict(&
246  "fileformat" / xdmf, &
247  "filename" / (trim(odir) // trim(ofname)), &
248  "count" / onum)
249 
250  config => dict(&
251  "physics" / physics, &
252  "fluxes" / fluxes, &
253  "mesh" / mesh, &
254  "boundary" / boundary, &
255  "sources" / sources, &
256  "timedisc" / timedisc, &
257  "datafile" / datafile)
258 
259  END SUBROUTINE makeconfig
260 
261 
262  SUBROUTINE initdata(Mesh,Physics,Timedisc,Fluxes,Sources,pvar,cvar)
263  USE physics_euler_mod, ONLY : physics_euler, statevector_euler
264  IMPLICIT NONE
265  !------------------------------------------------------------------------!
266  CLASS(mesh_base), INTENT(IN) :: Mesh
267  CLASS(physics_base), INTENT(INOUT) :: Physics
268  CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
269  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
270  CLASS(sources_base), POINTER :: Sources
271  CLASS(marray_compound),INTENT(INOUT):: pvar,cvar
272  !------------------------------------------------------------------------!
273  ! Local variable declaration
274  CLASS(sources_base), POINTER :: sp
275  CLASS(sources_gravity), POINTER :: gp
276  INTEGER :: i,j,k
277 #ifdef PARALLEL
278  INTEGER :: ierror
279 #endif
280  REAL :: mass
281  REAL, DIMENSION(:,:,:), POINTER :: r, Sigma
282  REAL, DIMENSION(:,:,:,:), POINTER :: r_vec
283  CHARACTER(LEN=32) :: mdisk_str
284  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
285  :: rands
286 #ifdef NECSXAURORA
287  INTEGER :: rng, n
288 #endif
289  !------------------------------------------------------------------------!
290  ! MESTEL'S DISK
291  ! distance from origin to cell bary centers and position vector
292  r => mesh%RemapBounds(mesh%radius%bcenter(:,:,:))
293  r_vec => mesh%RemapBounds(mesh%posvec%bcenter(:,:,:,:))
294  ! pointer to density array
295  sigma => mesh%RemapBounds(pvar%data4d(:,:,:,physics%DENSITY))
296 
297  ! set surface density using radial power law (1/r) with a little noise
298 #ifndef NECSXAURORA
299  CALL initrandseed(timedisc)
300  CALL random_number(rands)
301 #else
302  CALL asl_random_create(rng, asl_randommethod_mt19937_64)
303  CALL asl_random_distribute_uniform(rng)
304  n = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)
305  CALL asl_random_generate_d(rng, n, rands)
306 #endif
307  rands = rands * noise * 2.0 + (1.0 - noise)
308  sigma = rands*(rmin/r(:,:,:))
309 
310 #ifdef NECSXAURORA
311  CALL asl_random_destroy(rng)
312 #endif
313 
314  ! determine disk mass
315  mass = sum(mesh%volume%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) * &
316  sigma(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
317 #ifdef PARALLEL
318  CALL mpi_allreduce(mpi_in_place,mass,1,default_mpi_real,mpi_sum, &
319  mesh%comm_cart,ierror)
320 #endif
321  ! rescale disk mass
322  sigma(:,:,:) = sigma(:,:,:) * mdisk / mass
323 
324  ! 2. azimuthal velocity: balance initial radial acceleration with centrifugal acceleration
325  SELECT TYPE(phys => physics)
326  TYPE IS(physics_euler)
327  pvar%data4d(:,:,:,physics%PRESSURE) = &
328  csiso*csiso*pvar%data4d(:,:,:,physics%DENSITY)/gamma
329  ! necessary in order to calculte height below in UpdateGravity
330  SELECT TYPE(pvarious => pvar)
331  CLASS IS(statevector_euler)
332  CALL phys%UpdateSoundSpeed(pvarious)
333  END SELECT
334  END SELECT
335 
336  sp => sources
337  DO
338  IF (ASSOCIATED(sp).EQV..false.) RETURN
339  SELECT TYPE(sp)
340  CLASS IS(sources_gravity)
341  gp => sp
342  EXIT
343  END SELECT
344  sp => sp%next
345  END DO
346 
347  CALL gp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
348 
349  ! reset velocities
350  pvar%data2d(:,physics%XVELOCITY:physics%YVELOCITY) = 0.0
351  ! ATTENTION: Don't use GetCentrifugalVelocity without the optional acceleration array!
352  ! This would yield undefined data, because GetCentrifugalVelocity calls ComputeRHS
353  ! which calls CenterBoundary. Since the FARFIELD boundary conditions are not
354  ! initialized at this stage (see below), the result is undefined.
355  pvar%data4d(:,:,:,physics%XVELOCITY:physics%XVELOCITY+physics%VDIM-1) = &
356  timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),gp%accel%data4d)
357 
358  IF (mesh%FARGO.EQ.2) &
359  timedisc%w(:,:) = sqrt(physics%constants%GN*mbh/r(:,mesh%JMIN,:))-mesh%OMEGA*r(:,mesh%JMIN,:)
360 
361  ! transform velocities to rotating frame
362  pvar%data4d(:,:,:,physics%YVELOCITY) = pvar%data4d(:,:,:,physics%YVELOCITY) &
363  - mesh%OMEGA*r(:,:,:)
364 
365  ! get conservative variables
366  CALL physics%Convert2Conservative(pvar,cvar)
367 
368  ! setting for custom boundary conditions (western boundary)
369  SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
370  CLASS IS (boundary_custom)
371  CALL bwest%SetCustomBoundaries(mesh,physics, &
372  (/custom_nograd,custom_outflow,custom_kepler,custom_nograd/))
373  END SELECT
374 
375  ! setting for custom boundary conditions (eastern boundary)
376  SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
377  CLASS IS (boundary_custom)
378  CALL beast%SetCustomBoundaries(mesh,physics, &
379  (/custom_nograd,custom_outflow,custom_kepler,custom_nograd/))
380  END SELECT
381 
382  CALL physics%Convert2Conservative(pvar,cvar)
383  ! print some information
384  WRITE (mdisk_str, '(ES8.2)') mdisk
385  CALL mesh%Info(" DATA-----> initial condition: Mestel's disk")
386  CALL mesh%Info(" disk mass: " // trim(mdisk_str))
387 
388  END SUBROUTINE initdata
389 
390 
391  SUBROUTINE initrandseed(Timedisc)
392  IMPLICIT NONE
393  !------------------------------------------------------------------------!
394  CLASS(timedisc_base),INTENT(IN) :: Timedisc
395  INTEGER :: i, n, clock
396  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
397  !------------------------------------------------------------------------!
398  ! Initialize random number generator with a seed based on the systems time
399  ! source: http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
400  CALL random_seed(size = n)
401  ALLOCATE(seed(n))
402  CALL system_clock(count=clock)
403  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
404 #ifdef PARALLEL
405  seed = seed + timedisc%GetRank()
406 #endif
407  CALL random_seed(put = seed)
408  DEALLOCATE(seed)
409  END SUBROUTINE initrandseed
410 
411 
412 END PROGRAM sgdisk
Generated by   doxygen 1.8.14