shearingsheet.f90

Program and data initialization for a shearing-box simulation with standard initialization

Author
Jannes Klee

This program runs the self-gravitating shearingsheet simulation first done by Gammie (2001) [gammie2001] . For initialization a constant background density \( \Sigma \) is used and the pressure is set in a way to fullfill the Toomre-criterion \( Q = 1 \). The initial velocities \( v_x, v_y \) are randomly distributed with sub-sonic values. Whether you use a fast cooling with \( \beta = 2.0 \) or a slow cooling \( \beta = 10.0 \), you see fragmentation or a settling into the gravitoturbulent state.

Simulation parameters
cooling parameter \( \beta \) \( 10 \text{ (slow)}, 2 \text{ (fast)} \)
resolution \( N_x \times N_y \) \( 1024 \times 1024 \)
flux solver \( \mathtt{KT} \)
time-discretization \( \mathtt{DORMAND-PRINCE} \)
limiter \( \mathtt{VANLEER} \)
boundaries (north, south) shearing
boundaries (east, west) periodic
heat capacity ratio \( \gamma \) \( 2.0 \)
Initial condition
density \( \Sigma \) \( 1.0 \)
pressure \( P \) \( \frac{2.5^2 \pi^2 G \Sigma^3}{\gamma \Omega^2} \)
velocity \( v_x \) \( q \Omega y + \delta v_x \)
velocity \( v_y \) \( \delta v_y \)
random velocities \( \delta v_x, \delta v_y \) \( < c_{\mathrm{s}} \)
slow cooling
fastcooling

You can find some movies showing the temporal evolution of the column density in the gallery.


Attention
Long runtime.

References:

  • [gammie2001] Charles F. Gammie (2001). "Nonlinear outcome of gravitational instability in cooling, gaseous disks" The Astrophysical Journal 553: 174-183
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: shearingsheet.f90 #
5 !# #
6 !# Copyright (C) 2015-2018 #
7 !# Jannes Klee <jklee@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 !----------------------------------------------------------------------------!
78 !----------------------------------------------------------------------------!
79 PROGRAM shearingsheet
80  USE fosite_mod
81 #ifdef NECSXAURORA
82  USE asl_unified
83 #endif
84 #ifdef PARALLEL
85 #ifdef HAVE_MPI_MOD
86  USE mpi
87 #endif
88 #endif
89  IMPLICIT NONE
90 #ifdef PARALLEL
91 #ifdef HAVE_MPI_H
92  include 'mpif.h'
93 #endif
94 #endif
95  !--------------------------------------------------------------------------!
96  ! general constants
97  REAL, PARAMETER :: GN = 1.0 ! grav. constant [GEOM] !
98  INTEGER, PARAMETER :: UNITS = geometrical ! needs to be consistent !
99  ! simulation parameter
100  REAL, PARAMETER :: OMEGA = 1.0 ! rotation at fid. point !
101  REAL, PARAMETER :: SIGMA0 = 1.0 ! mean surf.dens. !
102  REAL, PARAMETER :: TSIM = 100./omega ! simulation time !
103  REAL, PARAMETER :: GAMMA = 2.0 ! dep. on vert. struct. !
104  REAL, PARAMETER :: BETA_C = 10.0 ! cooling parameter !
105 ! REAL, PARAMETER :: BETA_C = 2.0 ! 2 -> collapse !
106  REAL, PARAMETER :: Q = 1.5 ! shearing parameter !
107  ! mesh settings
108  INTEGER, PARAMETER :: MGEO = cartesian
109  INTEGER, PARAMETER :: SHEAR_DIRECTION = 1 ! enables shearingsheet with direcion !
110  INTEGER, PARAMETER :: XRES = 1024 ! cells in x-direction !
111  INTEGER, PARAMETER :: YRES = 1024 ! cells in y-direction !
112  INTEGER, PARAMETER :: ZRES = 1 ! cells in z-direction !
113  REAL :: DOMAINX = 320.0 ! domain size [GEOM] !
114  REAL :: DOMAINY = 320.0 ! domain size [GEOM] !
115  ! number of output time steps
116  INTEGER, PARAMETER :: ONUM = 100
117  ! output directory and output name
118  CHARACTER(LEN=256), PARAMETER :: ODIR = "./"
119  CHARACTER(LEN=256), PARAMETER :: OFNAME = "shearingsheet"
120  !--------------------------------------------------------------------------!
121  CLASS(fosite), ALLOCATABLE :: Sim
122  !--------------------------------------------------------------------------!
123 
124 ALLOCATE(sim)
125 
126 #ifdef NECSXAURORA
127 CALL asl_library_initialize()
128 #endif
129 
130 CALL sim%InitFosite()
131 CALL makeconfig(sim, sim%config)
132 CALL sim%Setup()
133 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc%pvar, sim%Timedisc%cvar)
134 CALL sim%Run()
135 CALL sim%Finalize()
136 
137 
138 #ifdef NECSXAURORA
139 CALL asl_library_finalize()
140 #endif
141 
142 DEALLOCATE(sim)
143 
144 CONTAINS
145  SUBROUTINE makeconfig(Sim,config)
146  IMPLICIT NONE
147  !--------------------------------------------------------------------------!
148  CLASS(fosite) :: Sim
149  TYPE(Dict_TYP), POINTER :: config
150  !--------------------------------------------------------------------------!
151  ! local variable declaration
152  TYPE(Dict_TYP), POINTER :: mesh,physics,fluxes,grav,cooling, &
153  shearingbox,sources,timedisc,datafile
154  REAL :: XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX, SOUNDSPEED
155  !--------------------------------------------------------------------------!
156  domainx = domainx*gn*sigma0/(omega*omega)
157  domainy = domainy*gn*sigma0/(omega*omega)
158  xmin = -0.5*domainx
159  xmax = +0.5*domainx
160  ymin = -0.5*domainy
161  ymax = +0.5*domainy
162  zmin = 0.0
163  zmax = 0.0
164  soundspeed = pi*gn*sigma0/omega ! Toomre-criterion
165 
166  ! physics settings
167  physics => dict(&
168  "problem" / euler, &
169  "gamma" / gamma, &
170  "units" / geometrical &
171  )
172 
173  ! mesh settings
174  mesh => dict(&
175  "meshtype" / midpoint, &
176  "geometry" / mgeo, &
177  "omega" / omega, &
178  "shear_dir" / shear_direction, &
179  "decomposition" / (/ 1, -1, 1/), &
180  "inum" / xres, &
181  "jnum" / yres, &
182  "knum" / zres, &
183  "xmin" / xmin, &
184  "xmax" / xmax, &
185  "ymin" / ymin, &
186  "ymax" / ymax, &
187  "zmin" / zmin, &
188  "zmax" / zmax, &
189  "Q" / q &
190  )
191 
192  ! fluxes settings
193  fluxes => dict(&
194  "order" / linear, &
195  "fluxtype" / kt, &
196  "variables" / primitive, &
197  "limiter" / vanleer &
198  )
199 
200  ! gravity settings (source term)
201  grav => dict(&
202  "stype" / gravity, &
203  "self/gtype" / sboxspectral, &
204  "self/output/phi" / 0, &
205  "self/output/accel_x" / 0, &
206  "self/output/Fmass2D" / 0, &
207  "self/output/accel_y" / 0 &
208  )
209 
210  ! parametrized cooling from Gammie (2001)
211  cooling => dict(&
212  "stype" / disk_cooling, &
213  "method" / gammie_sb, &
214  "output/Qcool" / 1, &
215  "b_cool" / beta_c &
216  )
217 
218  ! shearing box fictious forces
219  shearingbox => dict(&
220  "stype" / shearbox &
221  )
222 
223  ! sources settings (contains source terms)
224  sources => dict(&
225  "cooling" / cooling, &
226  "shearing" / shearingbox, &
227  "grav" / grav &
228  )
229 
230  ! time discretization settings
231  timedisc => dict(&
232  "method" / dormand_prince, &
233  "cfl" / 0.4, &
234  "stoptime" / tsim, &
235  "dtlimit" / 1e-40, &
236  "output/external_sources" / 0, &
237  "output/density_cvar" / 0, &
238  "output/energy" / 0, &
239  "output/xmomentum" / 0, &
240  "output/ymomentum" / 0, &
241  "output/rhs" / 0, &
242  "output/geometrical_sources" / 0, &
243  "maxiter" / 100000000, &
244  "tol_rel" / 1.0e-3 &
245  )
246 
247  ! data i/o settings
248  datafile => dict(&
249  "fileformat" / xdmf, &
250  "filepath" / trim(odir), &
251  "filename" / trim(ofname), &
252  "count" / onum &
253  )
254 
255  ! overall config settings
256  config => dict(&
257  "mesh" / mesh, &
258  "physics" / physics, &
259  "fluxes" / fluxes, &
260  "sources" / sources, &
261  "timedisc" / timedisc, &
262  "datafile" / datafile &
263  )
264  END SUBROUTINE makeconfig
265 
266  SUBROUTINE initdata(Mesh,Physics, pvar, cvar)
267  IMPLICIT NONE
268  !------------------------------------------------------------------------!
269  CLASS(mesh_base), INTENT(IN) :: Mesh
270  CLASS(physics_base), INTENT(IN) :: Physics
271  CLASS(marray_compound), POINTER, INTENT(INOUT) :: pvar,cvar
272  !------------------------------------------------------------------------!
273  REAL :: SOUNDSPEED
274  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) &
275  :: rands
276 #ifdef NECSXAURORA
277  INTEGER :: rng, n
278 #endif
279  !------------------------ standard run ----------------------------------!
280  SELECT TYPE(p => pvar)
281  TYPE IS(statevector_euler)
282  ! constant initial density
283  p%density%data3d(:,:,:) = sigma0
284 
285  ! constant initial pressure determined by Q = 1
286  soundspeed = pi*physics%Constants%GN*sigma0/omega ! Toomre-criterion
287  p%pressure%data3d(:,:,:) = 2.5**2*pi*pi* &
288  physics%Constants%GN**2.*sigma0**3./(gamma*omega**2)
289 
290  ! create random numbers for setup of initial velocities
291 #ifndef NECSXAURORA
292  CALL initrandseed(physics)
293 #else
294  CALL asl_random_create(rng, asl_randommethod_mt19937_64)
295  CALL asl_random_distribute_uniform(rng)
296  n = (mesh%IGMAX-mesh%IGMIN+1)*(mesh%JGMAX-mesh%JGMIN+1)*(mesh%KGMAX-mesh%KGMIN+1)
297 #endif
298  IF (mesh%shear_dir.EQ.2) THEN
299 #ifndef NECSXAURORA
300  CALL random_number(rands)
301 #else
302  CALL asl_random_generate_d(rng, n, rands)
303 #endif
304  rands = (rands-0.5)*0.1*soundspeed
305  p%velocity%data4d(:,:,:,1) = rands(:,:,:)
306 
307 #ifndef NECSXAURORA
308  CALL random_number(rands)
309 #else
310  CALL asl_random_generate_d(rng, n, rands)
311 #endif
312  rands = (rands-0.5)*0.1*soundspeed
313  p%velocity%data4d(:,:,:,2) = -q*omega*mesh%bcenter(:,:,:,1) + rands(:,:,:)
314  ELSE IF (mesh%shear_dir.EQ.1) THEN
315 #ifndef NECSXAURORA
316  CALL random_number(rands)
317 #else
318  CALL asl_random_generate_d(rng, n, rands)
319 #endif
320  rands = (rands-0.5)*0.1*soundspeed
321  p%velocity%data4d(:,:,:,1) = q*omega*mesh%bcenter(:,:,:,2) + rands(:,:,:)
322 
323 #ifndef NECSXAURORA
324  CALL random_number(rands)
325 #else
326  CALL asl_random_generate_d(rng, n, rands)
327 #endif
328  rands = (rands-0.5)*0.1*soundspeed
329  p%velocity%data4d(:,:,:,2) = rands(:,:,:)
330  END IF
331  CLASS DEFAULT
332  CALL physics%Error("shear::InitData","only non-isothermal HD supported")
333  END SELECT
334 #ifdef NECSXAURORA
335  CALL asl_random_destroy(rng)
336 #endif
337  !------------------------------------------------------------------------!
338  CALL physics%Convert2Conservative(pvar,cvar)
339  CALL mesh%Info(" DATA-----> initial condition: " // &
340  "Standard run shearingsheet")
341  END SUBROUTINE initdata
342 
344  SUBROUTINE initrandseed(Physics)
345  IMPLICIT NONE
346  !------------------------------------------------------------------------!
347  CLASS(physics_base),INTENT(IN) :: Physics
348  INTEGER :: i, n, clock
349  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
350  !------------------------------------------------------------------------!
351  ! Initialize random number generator with a seed based on the systems time
352  ! source: http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
353  CALL random_seed(size = n)
354  ALLOCATE(seed(n))
355  CALL system_clock(count=clock)
356  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
357 #ifdef PARALLEL
358  seed = seed + physics%GetRank()
359 #endif
360  CALL random_seed(put = seed)
361  DEALLOCATE(seed)
362  END SUBROUTINE initrandseed
363 END PROGRAM shearingsheet
Generated by   doxygen 1.8.14