33 real(8):: mass=0d0, mom(3)
34 real(8):: acceleration(3)=0d0
35 real(8),
dimension(:,:) ,
allocatable:: p, v, a
36 real(4),
dimension(:,:,:),
allocatable:: r
43 procedure:: init_task_list
50 procedure:: check_refine
51 procedure,
private:: refine_check
55 logical:: do_translate=.true.
56 logical:: force_refine=.true.
57 logical:: only_high_levels=.true.
59 integer:: center_star=0
61 integer:: max_sinks=10000
62 real:: accretion_rate=0.01
63 real:: accretion_efficiency=0.5
64 real:: accretion_fraction=1e-6
65 real:: accretion_radius=4.
66 real:: exclusion_radius=32.
67 real:: jeans_resolution=2.
69 real:: rho_limit_factor=8.
70 real:: rho_fractionr=1e-6
72 real:: softening_length=1.0
73 real(8):: out_time=0.1
74 character(len=64):: file =
'sink.dat' 76 integer,
save:: n_sinks=0
84 SUBROUTINE init (self)
87 logical,
save:: first_time=.true.
89 namelist /sink_patch_params/ &
90 on, do_translate, force_refine, only_high_levels, center_star, verbose, &
91 accretion_rate, accretion_efficiency, accretion_fraction, accretion_radius, &
92 exclusion_radius, jeans_resolution, rho_limit, rho_limit_factor, rho_cut, &
93 max_sinks, softening_length, out_time, file
95 call trace%begin (
'sink_patch_t%init')
100 read (io_unit%input, sink_patch_params, iostat=iostat)
101 if (io_unit%master)
write (io_unit%output, sink_patch_params)
105 self%particle_solver => particle_solver
112 SUBROUTINE dealloc (self)
115 call trace%begin (
'sink_patch_t%dealloc')
116 if (
allocated(self%r))
then 117 call io%bits_mem (-storage_size(self%r), product(shape(self%r)),
'-sink%r')
120 if (
allocated(self%p))
then 121 call io%bits_mem (-storage_size(self%p), product(shape(self%p)),
'-sink%p')
124 if (
allocated(self%v))
then 125 call io%bits_mem (-storage_size(self%v), product(shape(self%v)),
'-sink%v')
128 if (
allocated(self%a))
then 129 call io%bits_mem (-storage_size(self%a), product(shape(self%a)),
'-sink%a')
133 END SUBROUTINE dealloc
138 SUBROUTINE init_task_list (self, task_list)
140 class(
list_t),
pointer:: task_list
142 self%task_list => task_list
143 END SUBROUTINE init_task_list
148 FUNCTION check (self, patch)
154 if (n_sinks >= max_sinks)
return 155 call trace%begin (
'sink_patch_t%check')
157 check = check .and. patch%fmaxval(patch%idx%d) > rho_limit
159 print *,
'sink_patch_t%check: dmax in', &
160 maxloc(patch%mem(:,:,:,patch%idx%d,patch%it,1))
163 print *,
'n_sinks =', n_sinks
171 SUBROUTINE init_task (self, patch)
175 class(
link_t),
pointer:: nbor, prev
176 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: d
177 real(8),
dimension(:),
pointer:: r
178 integer:: i, l(3), u(3), m(3), ix, iy, iz
181 call trace%begin (
'sink_patch_t%init_task')
182 allocate (self%p(3,self%nt))
183 allocate (self%v(3,self%nt))
184 allocate (self%a(3,self%nt))
185 call io%bits_mem (storage_size(self%p), product(shape(self%p)),
'sink%p')
186 call io%bits_mem (storage_size(self%v), product(shape(self%v)),
'sink%v')
187 call io%bits_mem (storage_size(self%a), product(shape(self%a)),
'sink%a')
190 d => patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),patch%idx%d,patch%it,1)
194 self%position(i) = patch%position(i) + r(m(i))
195 self%ds(i) = patch%mesh(i)%d
197 self%n = 2*ceiling(accretion_radius)
205 self%no_mans_land = .false.
206 self%n = (self%n/2)*2 + 1
207 self%size = self%ds*(self%n-1)
212 call self%patch_t%setup (self%position, self%size, self%n, self%ng, self%nt, self%nv, self%nw)
213 if (self%id==io%id_debug)
then 214 print *,
'sink_patch_t%init_task: id, patch%id, m, dmax', self%id, patch%id, m, maxval(d)
215 print *,
'sink_patch_t%init_task: (pos-patch%pos)/size', (self%position-patch%position)/patch%size
216 print *,
'sink_patch_t%init_task: patch%pos, d', patch%myposition(m), patch%mem(m(1),m(2),m(3),patch%idx%d,patch%it,1)
217 print *,
'sink_patch_t%init_task: self%n , patch%n ', self%n , patch%n
220 self%level = patch%level
221 self%time = patch%time
223 self%mesh(i)%h = patch%mesh(i)%h
227 self%unsigned(self%idx%d) = .true.
228 call self%set (bits%internal)
232 allocate (self%mask(self%mesh(1)%gn,self%mesh(2)%gn,self%mesh(3)%gn))
233 allocate (self%r (self%mesh(1)%gn,self%mesh(2)%gn,self%mesh(3)%gn))
234 call io%bits_mem (storage_size(self%r ), product(shape(self%r )),
'sink%r')
235 call io%bits_mem (storage_size(self%mask), product(shape(self%mask)),
'sink%mask')
236 do iz=1,self%mesh(3)%gn
237 do iy=1,self%mesh(2)%gn
238 do ix=1,self%mesh(1)%gn
239 r2 = (ix-self%mesh(1)%o)**2 &
240 + (iy-self%mesh(2)%o)**2 &
241 + (iz-self%mesh(3)%o)**2
242 self%r(ix,iy,iz) = max(sqrt(r2),1e-3)
243 self%mask(ix,iy,iz) = merge(.true., .false., r2<accretion_radius**2)
248 END SUBROUTINE init_task
255 SUBROUTINE dnload (self, only)
257 integer,
optional:: only
259 call trace%begin (
'sink_patch_t%dnload')
260 call download%download_link (self%link, all_cells=.true., only=self%idx%d)
261 call download%download_link (self%link, all_cells=.true., only=self%idx%px)
262 call download%download_link (self%link, all_cells=.true., only=self%idx%py)
263 call download%download_link (self%link, all_cells=.true., only=self%idx%pz)
265 END SUBROUTINE dnload
269 SUBROUTINE update (self)
271 call trace%begin (
'sink_patch_t%update')
272 call self%particle_solver%force_field (self)
275 self%mem(:,:,:,:,self%new,1) = self%mem(:,:,:,:,self%it,1)
277 END SUBROUTINE update
281 SUBROUTINE accrete (self)
284 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: d, px, py, pz
285 real(kind=KindScalarVar),
dimension(:,:,:),
allocatable:: u2
291 call trace%begin (
'sink_patch_t%accrete')
292 allocate (u2(self%mesh(1)%gn,self%mesh(2)%gn,self%mesh(3)%gn))
296 d => self%mem(:,:,:,self%idx%d ,self%it,1)
297 px => self%mem(:,:,:,self%idx%px,self%it,1)
298 py => self%mem(:,:,:,self%idx%py,self%it,1)
299 pz => self%mem(:,:,:,self%idx%pz,self%it,1)
303 self%u_max = sqrt(self%fmaxval(u2))
305 call self%courant_condition
307 dv = product(self%ds)
310 if (self%mass==0d0)
then 311 do iz=self%mesh(3)%li,self%mesh(3)%ui
312 do iy=self%mesh(2)%li,self%mesh(2)%ui
313 do ix=self%mesh(1)%li,self%mesh(1)%ui
314 q = min(rho_limit/d(ix,iy,iz),1d0)
315 q = merge(q, 1d0, self%mask(ix,iy,iz))
317 if (self%mask(ix,iy,iz))
then 318 dm = dm + d(ix,iy,iz)*dv*p
320 print
'(3i4,3f12.6,l4)', ix, iy, iz, d(ix,iy,iz), self%r(ix,iy,iz), p, self%mask(ix,iy,iz)
321 self%mom(1) = self%mom(1) + px(ix,iy,iz)*dv*p
322 self%mom(2) = self%mom(2) + py(ix,iy,iz)*dv*p
323 self%mom(3) = self%mom(3) + pz(ix,iy,iz)*dv*p
324 d(ix,iy,iz) = d(ix,iy,iz)*q
325 px(ix,iy,iz) = px(ix,iy,iz)*q
326 py(ix,iy,iz) = py(ix,iy,iz)*q
327 pz(ix,iy,iz) = pz(ix,iy,iz)*q
338 do iz=1,self%mesh(3)%gn
339 do iy=1,self%mesh(2)%gn
340 do ix=1,self%mesh(1)%gn
341 omega = sqrt(scaling%grav*self%mass/self%r(ix,iy,iz))/self%r(ix,iy,iz)
342 p = merge(accretion_rate*self%dtime*omega, 0d0, self%mask(ix,iy,iz))
344 dm = dm + d(ix,iy,iz)*dv*p
345 self%mom(1) = self%mom(1) + px(ix,iy,iz)*dv*p
346 self%mom(2) = self%mom(2) + py(ix,iy,iz)*dv*p
347 self%mom(3) = self%mom(3) + pz(ix,iy,iz)*dv*p
348 d(ix,iy,iz) = d(ix,iy,iz)*q
349 px(ix,iy,iz) = px(ix,iy,iz)*q
350 py(ix,iy,iz) = py(ix,iy,iz)*q
351 pz(ix,iy,iz) = pz(ix,iy,iz)*q
357 self%mass = self%mass + dm
358 self%velocity = self%mom/self%mass
360 print
'(4(a,1p,e12.4,2x),2(a,3e12.4,2x))', &
361 'sink: time=', self%time,
'mass=', self%mass, &
362 'dmin=',self%fminval(self%idx%d),
'dmax=',self%fmaxval(self%idx%d), &
363 'vel=',self%velocity,
'pos=',self%position
365 END SUBROUTINE accrete
370 SUBROUTINE move (self)
372 real(8):: acceleration(3)
373 real,
dimension(:,:,:),
pointer:: phi
374 integer:: ix, iy, iz, ii(3)
376 call trace%begin (
'sink_patch_t%move')
381 ii = floor(self%mesh%o)
382 ix = ii(1); iy=ii(2); iz=ii(3)
383 self%acceleration = 0.0
387 self%velocity = self%velocity + 0.5d0*self%acceleration*self%dtime
391 self%time = self%time + self%dtime
392 self%position = self%position + self%dtime*self%velocity
397 self%acceleration = 0.0
401 self%velocity = self%velocity + 0.5d0*self%acceleration*self%dtime
407 INTEGER FUNCTION check_refine (self, patch)
412 call self%refine_check (patch%link, check_refine)
419 SUBROUTINE refine_check (self, link, refine)
421 class(
link_t),
pointer:: link
424 logical:: check_sinkparticle
425 type(
list_t),
pointer:: task_list
426 class(
task_t),
pointer:: task, patch
428 class(
link_t),
pointer:: new_link, nbor
431 real(kind=KindScalarVar),
dimension(:,:,:),
pointer:: d
435 check_sinkparticle = .false.
436 if (.not.sink_patch%on)
return 437 task_list => self%task_list
444 call task_list%refresh_nbors(link)
447 check_sinkparticle = sink_patch%check (gpatch)
448 if (check_sinkparticle)
then 449 do_trace = io%do_trace
457 new_link%parent => link
458 new_link%task => task
464 task%link => new_link
465 task%particle_solver => self%particle_solver
466 task%particle => particle
467 call self%particle_solver%append (particle)
468 call task%init_task (gpatch)
472 call task_list%init_nbors (new_link)
473 call task_list%append_link (new_link)
474 nbor => new_link%nbor
475 do while (
associated(nbor))
476 print *,
'nbor:', nbor%task%id,
real((nbor%task%position-task%position)/task%ds)
483 d => task%mem(:,:,:,task%idx%d,task%it,1)
484 d = max(d,task%fminval(d))
485 if (task%id==io%id_debug)
then 486 m = nint(task%offset)
487 print *,
'sink_patch_t%create: task%pos, d', task%myposition(m), &
488 task%mem(m(1),m(2),m(3),task%idx%d,task%it,1)
489 print *,
'MM', minval(d), task%fminval(task%idx%d), &
490 maxval(d), task%fmaxval(task%idx%d)
501 do while (
associated(nbor))
502 call task_list%init_nbors (nbor%link)
505 call task_list%init_nbors (link)
506 io%do_trace = do_trace
509 END SUBROUTINE refine_check
download_link: takes care of downloads to linktask same: called for patches on the same level differ:...
A communications hub between patches. Three methods are used: download_same, download_higher and down...
The gpath_t layer now essentially only handles restarts.
Module with list handling for generic class task_t objects.
Define code units, in terms of (here) CGS units.
Particle data type, extends a dll_node data type, so it can become part of a particle_list_t data typ...
Doubly linked list (DLL), carrying anything, as simply as possible.
A sink particle is represented as primarily a particle task, whose position is updated with a particl...
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Module with list handling for generic class task_t objects.
Kick-Drift-Kick (KDK) N-body solver.
Template module for tasks.