42 REAL,
PARAMETER :: gn = 1.0
43 REAL,
PARAMETER :: au = 1.0
44 REAL,
PARAMETER :: msun = 1.0
46 REAL,
PARAMETER :: tsim = 10
47 REAL,
PARAMETER :: valpha = 1e-3
51 REAL,
PARAMETER :: mbh1 = 0.4465*msun
52 REAL,
PARAMETER :: mbh2 = 0.4465*msun
54 REAL,
PARAMETER :: excent = 0.0
55 REAL,
PARAMETER :: semma = 0.224*au
57 REAL,
PARAMETER :: rgap = 2.5*semma
61 REAL,
PARAMETER :: mu = 2.35e-3
62 REAL,
PARAMETER :: rg = 8.31447
64 REAL,
PARAMETER :: gpar = au
65 REAL,
PARAMETER :: rmin = 1.5*semma
66 REAL,
PARAMETER :: rmax = 5.0*au
67 INTEGER,
PARAMETER :: xres = 128
68 INTEGER,
PARAMETER :: yres = 128
69 INTEGER,
PARAMETER :: zres = 1
71 INTEGER,
PARAMETER :: onum = 100
72 CHARACTER(LEN=256),
PARAMETER &
74 CHARACTER(LEN=256),
PARAMETER &
77 CLASS(
fosite),
ALLOCATABLE :: sim
85CALL initdata(sim%Mesh, sim%Physics, sim%Timedisc, sim%Fluxes, sim%Sources)
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, &
102 omega = sqrt(gn*(mbh1+mbh2)/semma)/semma
103 period = 2.*pi /
omega
109 "meshtype" / midpoint, &
110 "geometry" / logcylindrical, &
115 "xmin" / log(rmin/gpar), &
116 "xmax" / log(rmax/gpar), &
125 "decomposition" / (/ -1, 1, 1/), &
127 IF (
omega.GT.0.0)
CALL setattr(mesh,
"omega",
omega)
131 "problem" / euler_isotherm, &
132 "units" / geometrical,&
133 "output/fcsound" / 1, &
134 "output/bccsound" / 1)
141 "southern" / periodic, &
142 "northern" / periodic, &
143 "bottomer" / reflecting,&
144 "topper" / reflecting &
152 "variables" / primitive, &
153 "limiter" / vanleer, &
159 "stype" / viscosity, &
160 "vismodel" / alpha, &
162 "dynconst" / valpha, &
166 "stype" / rotating_frame)
170 "gtype" / pointmass_binary, &
175 "output/binpos" / 1, &
176 "output/omega" / 1, &
177 "semimajoraxis" / semma)
185 "output/height" / 1, &
186 "output/potential"/ 1, &
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 /), &
204 "maxiter" / 100000000)
208 "fileformat" / vtk, &
209 "filename" / (trim(odir) // trim(ofname)), &
213 "physics" / physics, &
216 "boundary" / boundary, &
217 "sources" / sources, &
218 "timedisc" / timedisc, &
219 "datafile" / datafile)
223 SUBROUTINE initdata(Mesh,Physics,Timedisc,Fluxes,Sources)
229 CLASS(mesh_base),
INTENT(IN) :: Mesh
230 CLASS(physics_base),
INTENT(INOUT) :: Physics
231 CLASS(timedisc_base),
INTENT(INOUT) :: Timedisc
232 CLASS(fluxes_base),
INTENT(INOUT) :: Fluxes
233 CLASS(sources_list),
ALLOCATABLE,
INTENT(INOUT) :: Sources
243 REAL,
DIMENSION(:,:,:),
POINTER :: r, Sigma
244 REAL,
DIMENSION(:,:,:,:),
POINTER :: r_faces
245 REAL,
DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: bccsound
246 REAL,
DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,6) :: fcsound
247 CHARACTER(LEN=32) :: info_str
250 IF (
ALLOCATED(sources)) &
251 sp => sources%GetSourcesPointer(gravity)
252 IF (.NOT.
ASSOCIATED(sp)) &
253 CALL physics%Error(
"mmsn::InitData",
"no gravity term initialized")
256 r => mesh%RemapBounds(mesh%radius%bcenter)
257 r_faces => mesh%RemapBounds(mesh%radius%faces)
260 bccsound = hratio*sqrt((mbh1+mbh2)*physics%constants%GN/r(:,:,:))
262 fcsound(:,:,:,i) = hratio*sqrt((mbh1+mbh2)*physics%constants%GN/r_faces(:,:,:,i))
266 SELECT TYPE (phys => physics)
268 CALL phys%SetSoundSpeeds(mesh,bccsound)
269 CALL phys%SetSoundSpeeds(mesh,fcsound)
272 CALL phys%Error(
"InitData",
"physics not supported")
277 SELECT TYPE(bwest => timedisc%Boundary%boundary(west)%p)
278 CLASS IS (boundary_custom)
279 CALL bwest%SetCustomBoundaries(mesh,physics, &
280 (/custom_nograd,custom_outflow,custom_kepler/))
282 SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
283 CLASS IS (boundary_custom)
284 CALL beast%SetCustomBoundaries(mesh,physics, &
285 (/custom_nograd,custom_outflow,custom_kepler/))
293 SELECT TYPE (pvar => timedisc%pvar)
296 pvar%density%data3d(:,:,:) =
fgap(r(:,:,:),rgap)*sigma0*0.1*xscale*r(:,:,:)**(-1.5)
303 CALL gp%UpdateGravity(mesh,physics,fluxes,pvar,0.0,0.0)
304 pvar%velocity%data4d(:,:,:,1:physics%VDIM) = &
305 timedisc%GetCentrifugalVelocity(mesh,physics,fluxes,sources,(/0.,0.,1./),gp%accel%data4d)
308 pvar%velocity%data4d(:,:,:,2) = pvar%velocity%data4d(:,:,:,2) - mesh%OMEGA*r(:,:,:)
310 CALL physics%Error(
"mmsn::InitData",
"unsupported state vector")
314 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
316 IF (mesh%fargo%GetType().EQ.2) &
317 timedisc%w(:,:) = sqrt(physics%constants%GN*(mbh1+mbh2)/mesh%radius%bcenter(:,mesh%JMIN,:))
320 CALL mesh%Info(
" DATA-----> initial condition: " //
"MMSN - Setup")
321 WRITE (info_str,
'(ES10.2)') xscale
322 CALL mesh%Info(
" disk mass: " // trim(info_str) //
" MMSN")
326 ELEMENTAL REAL FUNCTION fgap(R, Rgap)
328 REAL,
INTENT(IN) :: r, rgap
330 fgap = 1./(1. + exp(-(r-rgap)/(0.1*rgap)))
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
subroutine makeconfig(Sim, config)
elemental real function fgap(R, Rgap)
physics module for 1D,2D and 3D isothermal Euler equations
generic source terms module providing functionaly common to all source terms
generic gravity terms module providing functionaly common to all gravity terms
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)