DISPATCH
sink_patch_mod.f90
1 !===============================================================================
2 !> A sink particle is represented as primarily a particle task, whose position
3 !> is updated with a particle solver. The patch is usually much smaller than an
4 !> MHD patch, and its values are downloaded from contributing MHD patches.
5 !>
6 !> A particular quality level indicates to the download_mod that all cells from
7 !> the sink particle patch should be downloaded to overlapping patches, thus
8 !> communicating back the removal of accreted mass and momentum.
9 !>
10 !> The procedure is called from refine_t%extras%check_refine which handles
11 !> creating a task list link, adding an nbor list, and adding the task to the
12 !> task list.
13 !===============================================================================
15  USE io_mod
16  USE trace_mod
17  USE gpatch_mod
18  USE index_mod
19  USE download_mod
20  USE data_hub_mod
21  USE link_mod
22  USE list_mod
23  USE bits_mod
24  USE scaling_mod
25  USE kinds_mod
26  USE task_mod
27  USE dll_mod
28  USE particle_mod
30  implicit none
31  private
32  type, public, extends (gpatch_t):: sink_patch_t
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
37  class(particle_solver_t), pointer:: particle_solver
38  class(dll_node_t), pointer:: particle
39  logical:: on
40  contains
41  procedure:: init
42  procedure:: dealloc
43  procedure:: init_task_list
44  procedure:: check
45  procedure:: init_task
46  procedure:: dnload
47  procedure:: update
48  procedure:: accrete
49  procedure:: move
50  procedure:: check_refine
51  procedure, private:: refine_check
52  end type
53  !-----------------------------------------------------------------------------
54  logical:: on=.false. ! turn feature on/off
55  logical:: do_translate=.true. ! make a transformation to the velocity of the sink
56  logical:: force_refine=.true. ! insist on refining down to levelmax at star
57  logical:: only_high_levels=.true. ! only accrete at high levels
58  integer:: id=0 ! sink id
59  integer:: center_star=0 ! index of star to zoom in to
60  integer:: verbose=2 ! verbosity level
61  integer:: max_sinks=10000 ! for testing
62  real:: accretion_rate=0.01 ! The rate of accretion, measured in Omega_Kepler
63  real:: accretion_efficiency=0.5 ! The accretion radius, measured in cells
64  real:: accretion_fraction=1e-6 ! The accretion radius, measured in cells
65  real:: accretion_radius=4. ! The accretion radius, measured in cells
66  real:: exclusion_radius=32. ! The radius, measured in cells, inside which no other sink are born
67  real:: jeans_resolution=2. ! Minimum, umber of cells per Jeans' length
68  real:: rho_limit=1e6 ! The limit above which sink_particles are created
69  real:: rho_limit_factor=8. ! Express rho_limit relative to Jeans' ladder top
70  real:: rho_fractionr=1e-6 ! The fration of rho_limit above which accretion happens
71  real:: rho_cut=1.0 ! ??
72  real:: softening_length=1.0 ! Softening of gravitational potential
73  real(8):: out_time=0.1 ! cadence of output into sink.dat
74  character(len=64):: file ='sink.dat'! name of sink data file
75  !-----------------------------------------------------------------------------
76  integer, save:: n_sinks=0 ! counter
77  type(sink_patch_t), public:: sink_patch
78 CONTAINS
79 
80 !===============================================================================
81 !> Initialize an instance of the sink_patch_t data type, with parameters for
82 !> sink creation read from a sink_patch_params namelist.
83 !===============================================================================
84 SUBROUTINE init (self)
85  class(sink_patch_t):: self
86  !.............................................................................
87  logical, save:: first_time=.true.
88  integer:: iostat
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
94  !-----------------------------------------------------------------------------
95  call trace%begin ('sink_patch_t%init')
96  !$omp critical (input_cr)
97  if (first_time) then
98  first_time = .false.
99  rewind(io_unit%input)
100  read (io_unit%input, sink_patch_params, iostat=iostat)
101  if (io_unit%master) write (io_unit%output, sink_patch_params)
102  end if
103  !$omp end critical (input_cr)
104  self%on = on
105  self%particle_solver => particle_solver
106  call trace%end()
107 END SUBROUTINE init
108 
109 !===============================================================================
110 !> Deallocate permanent allocatables
111 !===============================================================================
112 SUBROUTINE dealloc (self)
113  class(sink_patch_t):: self
114  !-----------------------------------------------------------------------------
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')
118  deallocate (self%r)
119  end if
120  if (allocated(self%p)) then
121  call io%bits_mem (-storage_size(self%p), product(shape(self%p)), '-sink%p')
122  deallocate (self%p)
123  end if
124  if (allocated(self%v)) then
125  call io%bits_mem (-storage_size(self%v), product(shape(self%v)), '-sink%v')
126  deallocate (self%v)
127  end if
128  if (allocated(self%a)) then
129  call io%bits_mem (-storage_size(self%a), product(shape(self%a)), '-sink%a')
130  deallocate (self%a)
131  end if
132  call trace%end()
133 END SUBROUTINE dealloc
134 
135 !===============================================================================
136 !> Allows storing a link to a task list in the sink_patch_t data type
137 !===============================================================================
138 SUBROUTINE init_task_list (self, task_list)
139  class(sink_patch_t):: self
140  class(list_t), pointer:: task_list
141  !.............................................................................
142  self%task_list => task_list
143 END SUBROUTINE init_task_list
144 
145 !===============================================================================
146 !> Test if a sink particle should be created in a given patch
147 !===============================================================================
148 FUNCTION check (self, patch)
149  class(sink_patch_t):: self
150  class(gpatch_t), target:: patch
151  logical:: check
152  !-----------------------------------------------------------------------------
153  check = .false.
154  if (n_sinks >= max_sinks) return
155  call trace%begin ('sink_patch_t%check')
156  check = .true.
157  check = check .and. patch%fmaxval(patch%idx%d) > rho_limit
158  if (check) then
159  print *, 'sink_patch_t%check: dmax in', &
160  maxloc(patch%mem(:,:,:,patch%idx%d,patch%it,1))
161  !$omp atomic
162  n_sinks = n_sinks+1
163  print *,'n_sinks =', n_sinks
164  end if
165  call trace%end()
166 END FUNCTION check
167 
168 !===============================================================================
169 !> Create a sink partícle
170 !===============================================================================
171 SUBROUTINE init_task (self, patch)
172  class(sink_patch_t):: self
173  class(gpatch_t), target:: patch
174  !.............................................................................
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
179  real:: r2, dmin
180  !-----------------------------------------------------------------------------
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')
188  l = patch%mesh%li
189  u = patch%mesh%ui
190  d => patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),patch%idx%d,patch%it,1)
191  m = maxloc(d)+l-1
192  do i=1,3
193  r => patch%mesh(i)%r
194  self%position(i) = patch%position(i) + r(m(i))
195  self%ds(i) = patch%mesh(i)%d
196  end do
197  self%n = 2*ceiling(accretion_radius)
198  self%ng = 2
199  self%nv = patch%nv
200  self%nt = patch%nt
201  self%nw = 1
202  !-----------------------------------------------------------------------------
203  ! We should use no_mans_land=.false., and an odd number for self%n
204  !-----------------------------------------------------------------------------
205  self%no_mans_land = .false.
206  self%n = (self%n/2)*2 + 1
207  self%size = self%ds*(self%n-1)
208  !-----------------------------------------------------------------------------
209  ! Set a high level on the sink particle patch, to force downloading of interior
210  ! data
211  !-----------------------------------------------------------------------------
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
218  end if
219  self%quality = 1024
220  self%level = patch%level
221  self%time = patch%time
222  do i=1,3
223  self%mesh(i)%h = patch%mesh(i)%h
224  end do
225  self%mass = 0d0
226  self%velocity = 0d0
227  self%unsigned(self%idx%d) = .true.
228  call self%set (bits%internal)
229  !-----------------------------------------------------------------------------
230  ! Pre-load radius and mask array for use in accretion calculations
231  !-----------------------------------------------------------------------------
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)
244  end do
245  end do
246  end do
247  call trace%end()
248 END SUBROUTINE init_task
249 
250 !===============================================================================
251 !> Download values to all the sinkparticle%patch cells, but only to the
252 !> density and momentum slots. This relies on that the nbor list is
253 !> comprehensive, including also the parent patch
254 !===============================================================================
255 SUBROUTINE dnload (self, only)
256  class(sink_patch_t):: self
257  integer, optional:: only
258  !-----------------------------------------------------------------------------
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)
264  call trace%end()
265 END SUBROUTINE dnload
266 
267 !===============================================================================
268 !===============================================================================
269 SUBROUTINE update (self)
270  class(sink_patch_t):: self
271  call trace%begin ('sink_patch_t%update')
272  call self%particle_solver%force_field (self)
273  call self%accrete
274  call self%move
275  self%mem(:,:,:,:,self%new,1) = self%mem(:,:,:,:,self%it,1)
276  call trace%end()
277 END SUBROUTINE update
278 
279 !===============================================================================
280 !===============================================================================
281 SUBROUTINE accrete (self)
282  class(sink_patch_t):: self
283  !.............................................................................
284  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d, px, py, pz
285  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: u2
286  real:: v2_kepler
287  integer:: ix, iy, iz
288  real:: omega, dv
289  real(8):: p, q, dm
290  !-----------------------------------------------------------------------------
291  call trace%begin ('sink_patch_t%accrete')
292  allocate (u2(self%mesh(1)%gn,self%mesh(2)%gn,self%mesh(3)%gn))
293  !-----------------------------------------------------------------------------
294  ! Short hand
295  !-----------------------------------------------------------------------------
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)
300  u2 = (px/d)**2
301  u2 = (py/d)**2 + u2
302  u2 = (px/d)**2 + u2
303  self%u_max = sqrt(self%fmaxval(u2))
304  deallocate (u2)
305  call self%courant_condition
306  dm = 0d0
307  dv = product(self%ds)
308  !print *, 'XX0', minval(d), maxval(d), self%fmaxval(d)
309  !print *, 'XX1', dv, self%mesh%n, self%mesh%gn
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))
316  p = 1d0-q
317  if (self%mask(ix,iy,iz)) then
318  dm = dm + d(ix,iy,iz)*dv*p
319  if (verbose>1) &
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
328  end if
329  end do
330  end do
331  end do
332  !print *, '1 id,dm:', self%id, dm
333  !-----------------------------------------------------------------------------
334  ! Remove a fraction p of per-volume quantities and keep q=(1-p)
335  !-----------------------------------------------------------------------------
336  else
337  !print *, 'XX', scaling%grav, self%mass, minval(self%r), maxval(self%r)
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))
343  q = 1d0-p
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
352  end do
353  end do
354  end do
355  !print *, '2 id,dm:', self%id, dm
356  end if
357  self%mass = self%mass + dm
358  self%velocity = self%mom/self%mass
359  if (verbose > 0) &
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
364  call trace%end()
365 END SUBROUTINE accrete
366 
367 !===============================================================================
368 !> Move the sink particle one time step self%dt, using kidk-drift-kick
369 !===============================================================================
370 SUBROUTINE move (self)
371  class(sink_patch_t):: self
372  real(8):: acceleration(3)
373  real, dimension(:,:,:), pointer:: phi
374  integer:: ix, iy, iz, ii(3)
375  !-----------------------------------------------------------------------------
376  call trace%begin ('sink_patch_t%move')
377  !-----------------------------------------------------------------------------
378  ! Kick. The sink is assumed to be in the exact grid position floor(mesh%o)
379  !-----------------------------------------------------------------------------
380  !phi = self%mem(:,:,:,self%idx%phi,self%it,1)
381  ii = floor(self%mesh%o)
382  ix = ii(1); iy=ii(2); iz=ii(3)
383  self%acceleration = 0.0
384  !self%acceleration = [(phi(ix+1,iy ,iz )-phi(ix-1,iy ,iz ))/(2.*self%ds(1)), &
385  ! (phi(ix ,iy+1,iz )-phi(ix ,iy-1,iz ))/(2.*self%ds(2)), &
386  ! (phi(ix ,iy ,iz+1)-phi(ix ,iy ,iz-1))/(2.*self%ds(3))]
387  self%velocity = self%velocity + 0.5d0*self%acceleration*self%dtime
388  !-----------------------------------------------------------------------------
389  ! Drift + update phi to the new space-time position
390  !-----------------------------------------------------------------------------
391  self%time = self%time + self%dtime
392  self%position = self%position + self%dtime*self%velocity
393  !call download%download_link (self%link, all_cells=.true., only=self%idx%phi)
394  !-----------------------------------------------------------------------------
395  ! Kick
396  !-----------------------------------------------------------------------------
397  self%acceleration = 0.0
398  !self%acceleration = [(phi(ix+1,iy ,iz )-phi(ix-1,iy ,iz ))/(2.*self%ds(1)), &
399  ! (phi(ix ,iy+1,iz )-phi(ix ,iy-1,iz ))/(2.*self%ds(2)), &
400  ! (phi(ix ,iy ,iz+1)-phi(ix ,iy ,iz-1))/(2.*self%ds(3))]
401  self%velocity = self%velocity + 0.5d0*self%acceleration*self%dtime
402  call trace%end()
403 END SUBROUTINE move
404 
405 !===============================================================================
406 !===============================================================================
407 INTEGER FUNCTION check_refine (self, patch)
408  class(sink_patch_t):: self
409  class(gpatch_t), pointer:: patch
410  !.............................................................................
411  check_refine = -1
412  call self%refine_check (patch%link, check_refine)
413 END FUNCTION
414 
415 !===============================================================================
416 !> Test for the need to create a new sink particle, and handle the task list
417 !> aspect of that; creating an nbor list, and adding the link to the task list.
418 !===============================================================================
419 SUBROUTINE refine_check (self, link, refine)
420  class(sink_patch_t):: self
421  class(link_t), pointer:: link
422  integer:: refine
423  !.............................................................................
424  logical:: check_sinkparticle
425  type(list_t), pointer:: task_list
426  class(task_t), pointer:: task, patch
427  class(gpatch_t), pointer:: gpatch
428  class(link_t), pointer:: new_link, nbor
429  class(sink_patch_t), pointer:: new_sink
430  class(dll_node_t), pointer:: particle
431  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d
432  integer:: m(3)
433  logical:: do_trace
434  !-----------------------------------------------------------------------------
435  check_sinkparticle = .false.
436  if (.not.sink_patch%on) return
437  task_list => self%task_list
438  patch => link%task
439  select type (patch)
440  class is (sink_patch_t)
441  !---------------------------------------------------------------------------
442  ! Refresh nbor list
443  !---------------------------------------------------------------------------
444  call task_list%refresh_nbors(link)
445  class is (gpatch_t)
446  gpatch => patch
447  check_sinkparticle = sink_patch%check (gpatch)
448  if (check_sinkparticle) then
449  do_trace = io%do_trace
450  !io%do_trace = .true.
451  !-------------------------------------------------------------------------
452  ! Allocate a new sink_patch_t instance, as the task of a new link
453  !-------------------------------------------------------------------------
454  allocate (new_link)
455  allocate (sink_patch_t:: task)
456  allocate (particle_t:: particle)
457  new_link%parent => link
458  new_link%task => task
459  !-------------------------------------------------------------------------
460  ! Create the new sinkparticle, using data from the MHD patch
461  !-------------------------------------------------------------------------
462  select type (task)
463  class is (sink_patch_t)
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)
469  !-------------------------------------------------------------------------
470  ! Initialize its neighbor list, and add it to the task list
471  !-------------------------------------------------------------------------
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)
477  nbor => nbor%next
478  end do
479  !------------------------------------------------------------------------
480  ! Download variables and accrete
481  !------------------------------------------------------------------------
482  call task%dnload
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)
491  end if
492  call task%accrete
493  end select
494  !-------------------------------------------------------------------------
495  ! Refresh the nbor lists of the parent and its nbors; the parent should
496  ! be done after the nbors, to avoid putting the child link on its nbor
497  ! list until the other ones are updated. All nbors that happen to have
498  ! overlap with the sink particle patch will have it as an nbor
499  !-------------------------------------------------------------------------
500  nbor => link%nbor
501  do while (associated(nbor))
502  call task_list%init_nbors (nbor%link)
503  nbor => nbor%next
504  end do
505  call task_list%init_nbors (link)
506  io%do_trace = do_trace
507  end if
508  end select
509 END SUBROUTINE refine_check
510 
511 END MODULE sink_patch_mod
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...
Definition: data_hub_mod.f90:5
The gpath_t layer now essentially only handles restarts.
Definition: gpatch_mod.f90:4
Module with list handling for generic class task_t objects.
Definition: list_mod.f90:4
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Particle data type, extends a dll_node data type, so it can become part of a particle_list_t data typ...
Definition: particle_mod.f90:5
Doubly linked list (DLL), carrying anything, as simply as possible.
Definition: dll_mod.f90:4
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...
Definition: index_mod.f90:7
Definition: io_mod.f90:4
Kick-Drift-Kick (KDK) N-body solver.
Template module for tasks.
Definition: task_mod.f90:4