mmsn.f90

Program and data initialization for circumbinary accretion disk as MMSN

Author
Jannes Klee
Tobias Illenseer
Anna Feiler
Roman Avramenko
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: bindisk.f90 #
5 !# #
6 !# Copyright (C) 2010-2017 #
7 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8 !# Anna Feiler <afeiler@astrophysik.uni-kiel.de> #
9 !# Roman Avramenko <ravramenko@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 !----------------------------------------------------------------------------!
37 !----------------------------------------------------------------------------!
38 PROGRAM init
39  USE fosite_mod
40  !--------------------------------------------------------------------------!
41  ! general constants
42  REAL, PARAMETER :: GN = 1.0 ! Newtons grav. constant !
43  REAL, PARAMETER :: AU = 1.0 ! astronomical unit !
44  REAL, PARAMETER :: MSUN = 1.0 ! solar mass [kg] !
45  ! simulation parameters
46  REAL, PARAMETER :: TSIM = 10 ! simulation time [binary orbits]
47  REAL, PARAMETER :: VALPHA = 1e-3 ! alpha viscosity parameter !
48  ! 1. binary system
49 ! REAL, PARAMETER :: MBH1 = 0.690*MSUN
50 ! REAL, PARAMETER :: MBH2 = 0.203*MSUN
51  REAL, PARAMETER :: MBH1 = 0.4465*msun
52  REAL, PARAMETER :: MBH2 = 0.4465*msun
53 ! REAL, PARAMETER :: EXCENT = 0.159 ! excentricity !
54  REAL, PARAMETER :: EXCENT = 0.0 ! excentricity !
55  REAL, PARAMETER :: SEMMA = 0.224*au ! semi mayor axis !
56  ! 2. disk
57  REAL, PARAMETER :: Rgap = 2.5*semma ! estimated gap size !
58  REAL :: XSCALE = 1.0 ! scaling factor
59  REAL :: HRATIO = 0.05 ! scaling factor
60  ! gas paramter
61  REAL, PARAMETER :: MU = 2.35e-3 ! mean molecular mass [kg/mol] !
62  REAL, PARAMETER :: RG = 8.31447 ! molar gas constant !
63  ! mesh settings
64  REAL, PARAMETER :: GPAR = au ! geometry scaling paramete !
65  REAL, PARAMETER :: RMIN = 1.5*semma ! inner radius of the disk !
66  REAL, PARAMETER :: RMAX = 5.0*au ! outer radius of the grid !
67  INTEGER, PARAMETER :: XRES = 128 ! x-resolution !
68  INTEGER, PARAMETER :: YRES = 128 ! y-resolution !
69  INTEGER, PARAMETER :: ZRES = 1 ! y-resolution !
70  ! output file parameter
71  INTEGER, PARAMETER :: ONUM = 100 ! number of output time st !
72  CHARACTER(LEN=256), PARAMETER &
73  :: ODIR = "./"
74  CHARACTER(LEN=256), PARAMETER & ! output data file name !
75  :: OFNAME = 'mmsn'
76  !--------------------------------------------------------------------------!
77  CLASS(fosite), ALLOCATABLE :: Sim
78  REAL :: OMEGA,PERIOD
79  !--------------------------------------------------------------------------!
80 
81 ALLOCATE(sim)
82 CALL sim%InitFosite()
83 CALL makeconfig(sim, sim%config)
84 CALL sim%Setup()
85 CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, sim%Fluxes, sim%Sources)
86 CALL sim%Run()
87 CALL sim%Finalize()
88 DEALLOCATE(sim)
89 
90 CONTAINS
91 
92  SUBROUTINE makeconfig(Sim, config)
93  IMPLICIT NONE
94  !------------------------------------------------------------------------!
95  CLASS(fosite), INTENT(INOUT) :: Sim
96  TYPE(Dict_TYP),POINTER :: config
97  TYPE(Dict_TYP),POINTER :: mesh, physics, fluxes, boundary,grav, &
98  sources, binary, vis, timedisc, datafile, &
99  rotframe
100  !------------------------------------------------------------------------!
101  ! some derived simulation parameters
102  omega = sqrt(gn*(mbh1+mbh2)/semma)/semma
103  period = 2.*pi / omega
104 ! OMEGA = 0.0
105 
106  ! mesh settings
107  ! stellar orbits must be inside the central hole of the mesh
108  mesh => dict( &
109  "meshtype" / midpoint, &
110  "geometry" / logcylindrical, &
111 ! "geometry" / CYLINDRICAL, &
112  "inum" / xres, &
113  "jnum" / yres, &
114  "knum" / zres, &
115  "xmin" / log(rmin/gpar), &
116  "xmax" / log(rmax/gpar), &
117 ! "xmin" / (RMIN/GPAR), &
118 ! "xmax" / (RMAX/GPAR), &
119  "ymin" / 0.0, &
120  "ymax" / (2.0*pi), &
121  "zmin" / 0.0, &
122  "zmax" / 0.0, &
123  "use_fargo" / 1, &
124  "fargo" / 1, &
125  "decomposition" / (/ -1, 1, 1/), &
126  "gparam" / gpar)
127  IF (omega.GT.0.0) CALL setattr(mesh,"omega",omega)
128 
129  ! physics settings
130  physics => dict( &
131  "problem" / euler_isotherm, &
132  "units" / geometrical,&
133  "output/fcsound" / 1, &
134  "output/bccsound" / 1)
135 
136 
137  ! boundary conditions
138  boundary => dict( &
139  "western" / custom,&
140  "eastern" / custom,&
141  "southern" / periodic, &
142  "northern" / periodic, &
143  "bottomer" / reflecting,&
144  "topper" / reflecting &
145  )
146 
147 
148  ! numerical fluxes and reconstruction method
149  fluxes => dict( &
150  "order" / linear, &
151  "fluxtype" / kt, &
152  "variables" / primitive, &
153  "limiter" / vanleer, &
154  "theta" / 1.2)
155 
156 
157  ! viscosity source term
158  vis => dict( &
159  "stype" / viscosity, &
160  "vismodel" / alpha, &
161  "cvis" / 0.1,&
162  "dynconst" / valpha, &
163  "output/dynvis" / 1)
164 
165  rotframe => dict( &
166  "stype" / rotating_frame)
167 
168  ! gravitational acceleration due to binary system
169  binary => dict( &
170  "gtype" / pointmass_binary, &
171  "mass1" / mbh1, &
172  "mass2" / mbh2, &
173  "mesh_rot" / 1, &
174 ! "excentricity" / EXCENT, &
175  "output/binpos" / 1, &
176  "output/omega" / 1, &
177  "semimayoraxis" / semma)
178 
179  ! source term due to all gravity terms
180  grav => dict( &
181  "stype" / gravity, &
182  "binary" / binary,&
183 ! "self/gtype" / SPECTRAL, &
184 ! "self/green" / 1, &
185  "output/height" / 1, &
186  "output/potential"/ 1, &
187  "output/accel" / 1)
188 
189 
190  sources => dict( &
191  "vis" / vis, &
192  "grav" / grav)
193 ! IF (OMEGA.GT.0.0) CALL SetAttr(sources,"rotframe",rotframe)
194 
195  ! time discretization settings
196  timedisc => dict( &
197  "method" / ssprk,&
198  "cfl" / 0.3, &
199  "stoptime" / (period*tsim), &
200  "dtlimit" / 1.0e-40,&
201  "tol_rel" / 1.0e-3, &
202  "tol_abs" / (/ 1.0e-08, 1., 1.0e-08 /), &
203  "rhstype" / 1, &
204  "maxiter" / 100000000)
205 
206  ! initialize data input/output
207  datafile => dict( &
208  "fileformat" / vtk, &
209  "filename" / (trim(odir) // trim(ofname)), &
210  "count" / onum)
211 
212  config => dict( &
213  "physics" / physics, &
214  "fluxes" / fluxes, &
215  "mesh" / mesh, &
216  "boundary" / boundary, &
217  "sources" / sources, &
218  "timedisc" / timedisc, &
219  "datafile" / datafile)
220  END SUBROUTINE makeconfig
221 
222 
223  SUBROUTINE initdata(Mesh,Physics,Timedisc,Fluxes,Sources)
224  USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
225  IMPLICIT NONE
226  !------------------------------------------------------------------------!
227  CLASS(mesh_base), INTENT(IN) :: Mesh
228  CLASS(physics_base), INTENT(INOUT) :: Physics
229  CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
230  CLASS(fluxes_base), INTENT(INOUT) :: Fluxes
231  CLASS(sources_base), POINTER :: Sources
232  !------------------------------------------------------------------------!
233  ! Local variable declaration
234  CLASS(sources_base), POINTER :: sp
235  CLASS(sources_gravity), POINTER :: gp
236  INTEGER :: i
237 #ifdef PARALLEL
238  INTEGER :: ierror
239 #endif
240  REAL :: Sigma0
241  REAL, DIMENSION(:,:,:), POINTER :: r, Sigma
242  REAL, DIMENSION(:,:,:,:), POINTER :: r_faces
243  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: bccsound
244  REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,6) :: fcsound
245  CHARACTER(LEN=32) :: info_str
246  !------------------------------------------------------------------------!
247  ! set some pointers for convenience
248  r => mesh%RemapBounds(mesh%radius%bcenter)
249  r_faces => mesh%RemapBounds(mesh%radius%faces)
250 
251  ! Set sound speed
252  bccsound = hratio*sqrt((mbh1+mbh2)*physics%constants%GN/r(:,:,:))
253  DO i =1,6
254  fcsound(:,:,:,i) = hratio*sqrt((mbh1+mbh2)*physics%constants%GN/r_faces(:,:,:,i))
255  END DO
256 
257  ! set isothermal sound speeds
258  SELECT TYPE (phys => physics)
259  CLASS IS(physics_eulerisotherm)
260  CALL phys%SetSoundSpeeds(mesh,bccsound)
261  CALL phys%SetSoundSpeeds(mesh,fcsound)
262  CLASS DEFAULT
263  ! abort
264  CALL phys%Error("InitData","physics not supported")
265  END SELECT
266 
267  ! boundary conditions
268  ! custom boundary conditions at western boundary if requested
269  SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
270  CLASS IS (boundary_custom)
271  CALL bwest%SetCustomBoundaries(mesh,physics, &
272  (/custom_nograd,custom_outflow,custom_kepler/))
273  END SELECT
274  SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
275  CLASS IS (boundary_custom)
276  CALL beast%SetCustomBoundaries(mesh,physics, &
277  (/custom_nograd,custom_outflow,custom_kepler/))
278  END SELECT
279 
280  ! get gravitational acceleration
281  sp => sources
282  DO
283  IF (ASSOCIATED(sp).EQV..false.) RETURN
284  SELECT TYPE(sp)
285  CLASS IS(sources_gravity)
286  gp => sp
287  EXIT
288  END SELECT
289  sp => sp%next
290  END DO
291  IF (.NOT.ASSOCIATED(sp)) CALL physics%Error("mmsn::InitData","no gravity term initialized")
292 
293  ! choose Sigma0 in a way that 2% of the mass of the binary is contained in 30 au
294 ! Sigma0 = 0.0008*(MBH1+MBH2)*SQRT(30.0/AU)/(4.*PI)
295 ! Sigma0 = 0.02*(MBH1+MBH2)/(4.*PI*(-(30.0/GPAR)**(-0.5)+(RMIN/GPAR)**(-0.5)))
296 ! Sigma0 = 0.02*(MBH1+MBH2)*SQRT(30.0/AU)/(4.*PI)
297  sigma0 = 0.00311
298  SELECT TYPE (pvar => timedisc%pvar)
299  CLASS IS(statevector_eulerisotherm)
300  ! 1. set surface density and initialize gravitational acceleration
301  pvar%density%data3d(:,:,:) = fgap(r(:,:,:),rgap)*sigma0*0.1*xscale*r(:,:,:)**(-1.5)
302 ! pvar%velocity%data4d(:,:,:,:) = 0.0
303 ! IF (Mesh%FARGO.GT.0) Timedisc%w(:,:) = 0.0
304 ! ! get conservative variables
305 ! CALL Physics%Convert2Conservative(Timedisc%pvar,Timedisc%cvar)
306  ! 2. azimuthal velocity: balance initial radial gravitational acceleration
307  ! with centrifugal acceleration
308  CALL gp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
309  pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
310  timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),gp%accel%data4d)
311  ! transform velocities to rotating frame
312  ! ATTENTION: this works only if the second velocity is the azimuthal velocity
313  pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) - mesh%OMEGA*r(:,:,:)
314  CLASS DEFAULT
315  CALL physics%Error("mmsn::InitData","unsupported state vector")
316  END SELECT
317 
318  ! get conservative variables
319  CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
320 
321  IF (mesh%FARGO.EQ.2) &
322  timedisc%w(:,:) = sqrt(physics%constants%GN*(mbh1+mbh2)/mesh%radius%bcenter(:,mesh%JMIN,:))
323 
324  ! print some information on stdout
325  CALL mesh%Info(" DATA-----> initial condition: " // "MMSN - Setup")
326  WRITE (info_str, '(ES8.2)') xscale
327  CALL mesh%Info(" disk mass: " // trim(info_str) // " MMSN")
328 
329  END SUBROUTINE initdata
330 
331  ELEMENTAL REAL FUNCTION fgap(R, Rgap)
332  IMPLICIT NONE
333  REAL, INTENT(IN) :: R, Rgap
334 
335  fgap = 1./(1. + exp(-(r-rgap)/(0.1*rgap)))
336 
337  END FUNCTION
338 
339 END PROGRAM init
Generated by   doxygen 1.8.14