DISPATCH
trace_particles_mod.f90
1 !===============================================================================
2 !> Trace particle module. Two classes of trace particles will be supported
3 !> (this is not yet fully implemented):
4 !>
5 !> 1) The bulk of the particles are used to accurate trace the motion of mass,
6 !> with a fair number of particles per cell, constantly renormalized to accurately
7 !> represent the mass in the patch.
8 !>
9 !> 2) A small fraction of trace particles, intended for graphics, are flagged
10 !> with a bit that results in a higher cadence of saving to disk, with particles
11 !> moving only a few cells between snapshots.
12 !>
13 !> Disk files should be organized with one file per snapshot for the bulk
14 !> particles, while the graphics particles should be accessible in one file
15 !> per rank, since it needs to be written too frequently for collective MPI
16 !> calls. When reading and plotting, one can open these rank files, and make
17 !> a selection, based on for example their first location. Subsequently, as
18 !> particles transition to other ranks, this has to be handled by the reader.
19 !>
20 !> One example: For superonic turbulence, one may choose to save only two
21 !> snapshots per dynamic time L/M, while that would correspond to at least N/2
22 !> cell transitions. Hence, the number of particles needs to be of order N
23 !> smaller in the graphics files, if disks files should have about the same
24 !> size. In principle, only positions and particle IDs need to be save in the
25 !> bulk snapshots, which means of order 16 bytes per particle, and one can then
26 !> allow a few particles per cell, w/o a major increase in data volume.
27 !>
28 !> Conversely, one can allow of order N times fewer particles in the graphics
29 !> files, or about one particle per 1000 cells, or one particle per radius 10.
30 !> One could choose to make the graphics particles have a constant fraction of
31 !> the mass, especially if one ensures that their paths accurately reflect the
32 !> the path of the mass. For such a small fraction it does not make sense to
33 !> even consider mass conservation, and one should instead extra precision in
34 !> their motion.
35 !>
36 !> The particle ID and a few flags may be encoded in a 64 bit integer, with the
37 !> upper 4 bits reserved for flags, being masked out when considering IDs.
38 !>
39 !> In the case of the bulk particles, we are interested in two partially
40 !> conflicting requirements: 1) they should be renomalized frequently, to keep
41 !> the weight of each particle from becoming too extreme (small or large)
42 !> compared with the amount of mass in the cell (or in the patch). 2) particle
43 !> heritage should be accessible, so one can identify exactly where a particle
44 !> comes from. This could be achieved by reserving the lowermost No bits for
45 !> their "original" ID.
46 !===============================================================================
48  USE io_mod
49  USE trace_mod
50  USE dll_mod
51  USE particle_mod
53  USE omp_timer_mod
54  USE random_mod
55  USE link_mod
56  USE patch_mod
57  USE task_mod
58  USE link_mod
59  USE mpi_buffer_mod
60  implicit none
61  private
62  integer, parameter:: nt=4
63  integer, parameter:: particle_kind=8
64  type, public, extends(particle_list_t):: trace_particles_t
65  real:: per_cell, per_mass
66  real:: lf(3), uf(3)
67  integer:: id, it, time_order, n_trade
68  integer:: offset(3)
69  real, allocatable:: v(:,:,:,:)
70  real(8), allocatable:: wt(:)
71  real(8):: position(3)
72  logical:: on
73  type(particle_list_t):: export, import
74  type(mpi_buffer_t):: mpi_buffer
75  class(patch_t), pointer:: patch
76  contains
77  procedure:: init
78  procedure:: dealloc
79  procedure:: add
80  procedure:: update
81  procedure:: post_update
82  procedure:: pack
83  procedure, private:: mpi_append_particles
84  procedure, private:: mpi_read_particles
85  end type
86  type(random_t):: random
87 CONTAINS
88 
89 !===============================================================================
90 !> Initialize trace particles
91 !===============================================================================
92 SUBROUTINE init (self, name)
93  class(trace_particles_t):: self
94  character(len=*), optional:: name
95  !-----------------------------------------------------------------------------
96  call trace%begin ('trace_particles_t%init')
97  call trace%end
98 END SUBROUTINE init
99 
100 !===============================================================================
101 !> Deallocate
102 !===============================================================================
103 SUBROUTINE dealloc (self)
104  class(trace_particles_t):: self
105  !-----------------------------------------------------------------------------
106  call trace%begin ('trace_particles_t%dealloc')
107  if (allocated(self%v)) then
108  call io%bits_mem (-storage_size(self%v), &
109  product(shape(self%v)), '-tracep%v')
110  deallocate (self%v)
111  end if
112  call trace%end
113 END SUBROUTINE dealloc
114 
115 !===============================================================================
116 !> Add trace particles, if requested by the namelist
117 !===============================================================================
118 SUBROUTINE add (self, link)
119  class(trace_particles_t):: self
120  class(link_t), pointer:: link
121  !.............................................................................
122  real, pointer:: d(:,:,:), pv(:,:,:,:)
123  integer:: iostat
124  integer:: i, ip, np, ix, iy, iz, no
125  integer, save:: time_order=3, verbose=0
126  real, save:: per_cell=0.0
127  real, save:: per_mass=1.0
128  logical, save:: on=.false.
129  logical, save:: first_time=.true.
130  namelist /trace_particle_params/ on, verbose, per_cell, per_mass, time_order
131  !-----------------------------------------------------------------------------
132  call trace%begin ('trace_particles_t%add')
133  self%patch => task2patch(link%task)
134  associate(patch => self%patch)
135  !$omp critical (input_cr)
136  if (first_time) then
137  first_time = .false.
138  rewind(io%input)
139  read (io%input, trace_particle_params, iostat=iostat)
140  if (io%master) &
141  write (io%output, trace_particle_params)
142  call random%init
143  end if
144  !$omp end critical (input_cr)
145  select type (patch)
146  class is (patch_t)
147  self%on = on
148  self%nt = patch%nt
149  self%per_cell = per_cell
150  self%time_order = time_order
151  self%verbose = verbose
152  if (.not.on) then
153  call trace%end ()
154  return
155  end if
156  !----------------------------------------------------------------------------
157  ! Allocate and initialize
158  !----------------------------------------------------------------------------
159  call self%particle_list_t%init
160  !-----------------------------------------------------------------------------
161  ! Initialize Adams-Bashforth integration; the concept of symplectic integration
162  ! does not apply when particles are just folling the gas, so there is no point
163  ! in doing kick-drift-kick for trace particles. One could possibly improve
164  ! accuracy by using the fact that one has both location and velocity history,
165  ! but a higher order cell velocity interpolation may be more important.
166  !-----------------------------------------------------------------------------
167  time_order = min(time_order,patch%nt)
168  allocate (self%wt(time_order))
169  select case (time_order)
170  case (2)
171  self%wt(1:2) = [1.5, -0.5]
172  case (3)
173  self%wt(1:3) = [23./12., -16./12., 5./12.]
174  case default
175  end select
176  !-----------------------------------------------------------------------------
177  !-----------------------------------------------------------------------------
178  d => patch%mem(:,:,:,patch%idx%d ,patch%it,1)
179  pv => patch%mem(:,:,:,patch%idx%px:patch%idx%pz,patch%it,1)
180  allocate (self%v(patch%gn(1),patch%gn(2),patch%gn(3),3))
181  call io%bits_mem (storage_size(self%v), &
182  product(shape(self%v)), 'tracep%v')
183  do i=1,3
184  self%v(:,:,:,i) = pv(:,:,:,i)/d(:,:,:)
185  end do
186  !-----------------------------------------------------------------------------
187  ! Add particles in proportion to the density, relative to the average
188  !-----------------------------------------------------------------------------
189  call self%particle_list_t%init ('trace')
190  self%id = patch%id
191  if (self%per_cell > 0.0) then
192  np = ifix(self%per_cell)
193  do iz=patch%mesh(3)%li,patch%mesh(3)%ui
194  do iy=patch%mesh(2)%li,patch%mesh(2)%ui
195  do ix=patch%mesh(3)%li,patch%mesh(1)%ui
196  do ip=1,ifix(self%per_cell)
197  call add1
198  end do
199  if (random%ran1() < (self%per_cell-np)) then
200  call add1
201  end if
202  end do
203  end do
204  end do
205  end if
206  end select
207  call io%bits_mem (-storage_size(self%v), &
208  product(shape(self%v)), '-tracep%v')
209  deallocate (self%v)
210  call trace%end ()
211  end associate
212 contains
213  !-----------------------------------------------------------------------------
214  real function ran()
215  ran = random%ran1() - 0.5
216  end function
217  !-----------------------------------------------------------------------------
218  subroutine add1
219  class(particle_t), pointer:: p
220  class(dll_node_t), pointer:: node
221  integer:: it
222  !---------------------------------------------------------------------------
223  associate(patch => self%patch)
224  allocate (p)
225  call p%init
226  do it=1,nt
227  p%r(:,it) = [ix+ran(),iy+ran(),iz+ran()]
228  p%v(:,it) = self%v(ix,iy,iz,:)*patch%dtime
229  p%t(it) = patch%time
230  end do
231  node => p
232  call self%append (node)
233  end associate
234  end subroutine add1
235  !-----------------------------------------------------------------------------
236 END SUBROUTINE add
237 
238 !===============================================================================
239 !> Update particle positions. To optimize, particle positions are already
240 !> normalized to floating point array coordinates. This allows using single
241 !> precision, also on disk, where the coordinates just need to be complemented
242 !> with the double precision patch position.
243 !===============================================================================
244 SUBROUTINE update (self)
245  class(trace_particles_t):: self
246  !.............................................................................
247  real(8):: time, dtime, dt
248  real, pointer:: d(:,:,:), pv(:,:,:,:)
249  integer:: it, o1, o2, o3, i ,j, nrem
250  class(dll_node_t), pointer:: node, next
251  class(particle_t), pointer:: p
252  class(*), pointer:: car
253  real(kind=particle_kind), pointer:: r(:)
254  integer, save:: itimer=0
255  !-----------------------------------------------------------------------------
256  if (.not. self%on) return
257  call trace%begin ('trace_particles_t%update', itimer=itimer)
258  associate(patch => self%patch)
259  time = patch%time
260  dtime = patch%dtime
261  !----------------------------------------------------------------------------
262  ! Update velocity
263  !----------------------------------------------------------------------------
264  d => patch%mem(:,:,:,patch%idx%d ,patch%it,1)
265  pv => patch%mem(:,:,:,patch%idx%px:patch%idx%pz,patch%it,1)
266  allocate (self%v(patch%gn(1),patch%gn(2),patch%gn(3),3))
267  do i=1,3
268  self%v(:,:,:,i) = pv(:,:,:,i)/d(:,:,:)
269  end do
270  !-----------------------------------------------------------------------------
271  ! Import any particles exported by neighbors
272  !-----------------------------------------------------------------------------
273  if (self%import%n > 0) then
274  node => self%import%head
275  do while (associated(node))
276  next => node%next
277  call self%import%remove (node)
278  call self%append (node)
279  node => next
280  end do
281  if (self%verbose > 0) &
282  write (io%output,*) self%id, ' has', self%n, ' particles'
283  end if
284  self%it = o1
285  !-----------------------------------------------------------------------------
286  ! Export and prepare import
287  !-----------------------------------------------------------------------------
288  self%lf = patch%mesh%li - 0.5
289  self%uf = patch%mesh%ui + 0.5
290  call self%export%init ('export')
291  call self%import%init ('import')
292  call interpolate (self)
293  call io%bits_mem (-storage_size(self%v), &
294  product(shape(self%v)), '-tracep%v')
295  deallocate (self%v)
296  !-----------------------------------------------------------------------------
297  ! Direct, non-vectorized update
298  !-----------------------------------------------------------------------------
299  node => self%head
300  nrem = 0
301  do while (associated(node))
302  select type (node)
303  class is (particle_t)
304  p => node
305  end select
306  it = 1+mod(p%it,nt)
307  o1 = 1+modulo(it-2,nt)
308  o2 = 1+modulo(it-3,nt)
309  dt = min(dtime,time-p%t(o1))
310  p%v(:,o1) = p%v(:,o1)*dt
311  if (self%time_order == 3) then
312  o3 = 1+modulo(it-4,nt)
313  p%r(:,it) = &
314  p%r(:,o1) + &
315  p%v(:,o1)*self%wt(1) + &
316  p%v(:,o2)*self%wt(2) + &
317  p%v(:,o3)*self%wt(3)
318  else
319  p%r(:,it) = &
320  p%r(:,o1) + &
321  p%v(:,o1)*self%wt(1) + &
322  p%v(:,o2)*self%wt(2)
323  end if
324  p%it = it
325  p%t(it) = time
326  !---------------------------------------------------------------------------
327  ! Particles outside the patch interior are appended to an export list and
328  ! removed from the patch particle list
329  !---------------------------------------------------------------------------
330  r => p%r(:,it)
331  if (any(r < self%lf .or. r > self%uf)) then
332  call self%remove (node)
333  call self%export%append (node)
334  nrem = nrem+1
335  end if
336  node => node%next
337  end do
338  if (self%verbose > 0 .and. nrem>0) &
339  write (io%output,*) self%id, ' has', self%n, ' particles', nrem
340  call trace%end (itimer)
341  end associate
342 END SUBROUTINE update
343 
344 !===============================================================================
345 !> Tri-linear interpolate in a velocity array. The particle positions within
346 !> a patch are maintained in units of mesh points, both to simplify interpolation
347 !> and updates, and to improve precision, while staying with singe precision
348 !> in coordinates and velocities.
349 !===============================================================================
350 SUBROUTINE interpolate (self)
351  class(trace_particles_t):: self
352  integer:: it
353  !.............................................................................
354  integer:: sz(4), ip(3), i, j
355  real:: p(3), q(3)
356  class(dll_node_t), pointer:: node
357  class(particle_t), pointer:: part
358  integer, save:: itimer=0
359  !-----------------------------------------------------------------------------
360  ! Straight loops. This could be optimized further, by using arrays instead
361  ! of lists, and accpting "holes" in the arrays up to a certain fraction. But
362  ! this is already fast enough, given the disk space contsraints on the number
363  ! of particles.
364  !-----------------------------------------------------------------------------
365  call trace%begin ('trace_particles_t%interpolate', itimer=itimer)
366  sz = size(self%v)
367  node => self%head
368  do while (associated(node))
369  select type (node)
370  type is (particle_t)
371  part => node
372  it = part%it
373  end select
374  p = 1.0 + part%r(:,it)
375  ip = floor(p)
376  ip = max(1,min(ip,sz(1:3)))
377  p = p-ip
378  do j=1,3
379  part%v(j,it) = &
380  q(3)*(q(2)*(q(1)*self%v(ip(1) ,ip(2) ,ip(3) ,j) &
381  + p(1)*self%v(ip(1)+1,ip(2) ,ip(3) ,j)) &
382  + p(2)*(q(1)*self%v(ip(1) ,ip(2)+1,ip(3) ,j) &
383  + p(1)*self%v(ip(1)+1,ip(2)+1,ip(3) ,j))) &
384  + p(3)*(q(2)*(q(1)*self%v(ip(1) ,ip(2) ,ip(3)+1,j) &
385  + p(1)*self%v(ip(1)+1,ip(2) ,ip(3)+1,j)) &
386  +p(2)*(q(1)*self%v(ip(1) ,ip(2)+1,ip(3)+1,j) &
387  + p(1)*self%v(ip(1)+1,ip(2)+1,ip(3)+1,j)))
388  end do
389  node => node%next
390  end do
391  call trace%end (itimer)
392 END SUBROUTINE interpolate
393 
394 !===============================================================================
395 !> Trade particles between patches.
396 !===============================================================================
397 SUBROUTINE trade (self, nbor)
398  class(trace_particles_t):: self
399  class(trace_particles_t):: nbor
400  class(dll_node_t), pointer:: node, next
401  class(particle_t), pointer:: part
402  real(kind=particle_kind):: r(3)
403  !-----------------------------------------------------------------------------
404  self%n_trade = 0
405  if (self%export%n == 0) return
406  node => self%export%head
407  do while (associated(node))
408  next => node%next
409  select type (node)
410  class is (particle_t)
411  part => node
412  end select
413  r = part%r(:,part%it) + nbor%offset
414  !---------------------------------------------------------------------------
415  ! Particles in the interior of an nbor node are removed from the export list
416  ! of this patch, and put on the import list of that patch, with adjusted r
417  !---------------------------------------------------------------------------
418  if (all(r >= self%lf .and. r <= self%uf)) then
419  if (self%verbose > 1) &
420  write (io%output,*) nbor%id, self%id, part%id, part%r(:,part%it), nbor%offset
421  part%r(:,part%it) = r
422  call self%export%remove (node)
423  call nbor%import%append (node)
424  self%n_trade = self%n_trade + 1
425  end if
426  node => node%next
427  end do
428 END SUBROUTINE trade
429 
430 !===============================================================================
431 !> Trade particles between patches.
432 !===============================================================================
433 SUBROUTINE post_update (self)
434  class(trace_particles_t):: self
435  integer:: i
436  class(link_t), pointer:: nbor
437  class(task_t), pointer:: nbtask
438  class(*), pointer:: nbpart
439  !----------------------------------------------------------------------------
440  ! If particles have been accumlated in the "export" particle list, loop over
441  ! nbors and transfer particles to their "import" particle lists. This is to
442  ! minimize the need for locking or critical regions. It should be enough to
443  ! protect the import lists with OMP locks, while they are being updated, and
444  ! while they are being read.
445  !----------------------------------------------------------------------------
446  associate(patch => self%patch)
447  if (self%export%n > 0) then
448  if (self%verbose > 0) &
449  write (io%output,*) &
450  patch%id, self%export%n, ' particles to export'
451  nbor => patch%link%nbor
452  do while (associated(nbor))
453  nbtask => nbor%task
454  select type (nbtask)
455  class is (patch_t)
456  nbpart => nbtask%connect%trace_particles
457  select type (nbpart)
458  class is (trace_particles_t)
459  nbpart%offset = nint((mod(patch%position-nbtask%position + &
460  1.5*patch%box,patch%box)-0.5*patch%box)/patch%ds)
461  call trade (self, nbpart)
462  end select
463  end select
464  nbor => nbor%next
465  end do
466  end if
467  end associate
468 END SUBROUTINE post_update
469 
470 !===============================================================================
471 !> Prepare an MPI buffer
472 !===============================================================================
473 SUBROUTINE pack (self)
474  class(trace_particles_t):: self
475  integer, pointer:: buffer(:)
476  real, allocatable:: r(:,:,:)
477  class(dll_node_t), pointer:: part
478  integer:: np
479  !-----------------------------------------------------------------------------
480  call trace%begin ('trace_particles_t%pack')
481  !-----------------------------------------------------------------------------
482  ! Prepare an MPI buffer
483  !-----------------------------------------------------------------------------
484  associate(patch => self%patch)
485  np = self%n + self%export%n
486  call self%mpi_buffer%init (np + 10, 'trace_particles_t')
487  call self%mpi_append_particles (self)
488  call self%mpi_append_particles (self%export)
489  call trace%end ()
490  end associate
491 END SUBROUTINE pack
492 
493 !===============================================================================
494 !> Prepare an MPI buffer
495 !===============================================================================
496 SUBROUTINE unpack (self)
497  class(trace_particles_t):: self
498  integer, pointer:: buffer(:)
499  real, allocatable:: r(:,:,:)
500  class(dll_node_t), pointer:: part
501  integer:: np
502  !-----------------------------------------------------------------------------
503  call trace%begin ('trace_particles_t%pack')
504  !-----------------------------------------------------------------------------
505  ! Prepare an MPI buffer
506  !-----------------------------------------------------------------------------
507  associate(patch => self%patch)
508  np = self%n + self%export%n
509  call self%mpi_buffer%init (np + 10, 'trace_particles_t')
510  call self%mpi_read_particles (self)
511  call self%mpi_read_particles (self%export)
512  call trace%end ()
513  end associate
514 END SUBROUTINE unpack
515 
516 !===============================================================================
517 !> Append particle data to an mpi_buffer
518 !===============================================================================
519 SUBROUTINE mpi_append_particles (self, particles)
520  class(trace_particles_t):: self
521  class(particle_list_t):: particles
522  real, allocatable:: r(:,:,:)
523  class(dll_node_t), pointer:: part
524  integer:: i
525  !-----------------------------------------------------------------------------
526  call trace%begin ('mpi_buffer%append_particles')
527  !-----------------------------------------------------------------------------
528  ! Run through the particle list and collect the data in an array
529  !-----------------------------------------------------------------------------
530  allocate (r(3,nt,self%n))
531  i = 0
532  part => self%head
533  do while (associated(part))
534  select type (part)
535  class is (particle_t)
536  i = i+1
537  r(:,:,i) = part%r
538  end select
539  part => part%next
540  end do
541  !-----------------------------------------------------------------------------
542  ! Append the buffer and deallocate
543  !-----------------------------------------------------------------------------
544  call self%mpi_buffer%append (r)
545  deallocate (r)
546  call trace%end ()
547 END SUBROUTINE mpi_append_particles
548 
549 !===============================================================================
550 !> Append particle data to an mpi_buffer
551 !===============================================================================
552 SUBROUTINE mpi_read_particles (self, particles)
553  class(trace_particles_t):: self
554  class(particle_list_t):: particles
555  real, allocatable:: r(:,:,:)
556  class(dll_node_t), pointer:: part
557  integer:: i
558  !-----------------------------------------------------------------------------
559  call trace%begin ('mpi_buffer%read_particles')
560  !-----------------------------------------------------------------------------
561  ! Run through the particle list and collect the data in an array
562  !-----------------------------------------------------------------------------
563  allocate (r(3,nt,self%n))
564  call self%mpi_buffer%read (r)
565  i = 0
566  part => self%head
567  do while (associated(part))
568  select type (part)
569  class is (particle_t)
570  i = i+1
571  part%r = r(:,:,i)
572  end select
573  part => part%next
574  end do
575  !-----------------------------------------------------------------------------
576  ! Append the buffer and deallocate
577  !-----------------------------------------------------------------------------
578  deallocate (r)
579  call trace%end ()
580 END SUBROUTINE mpi_read_particles
581 
582 END MODULE trace_particles_mod
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Trace particle module. Two classes of trace particles will be supported (this is not yet fully implem...
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
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Doubly linked list (DLL), carrying anything, as simply as possible.
Definition: dll_mod.f90:4
Definition: io_mod.f90:4
Particle list, extending a doubly-linked list. Each particle maintains arrays with previous positions...
Template module for tasks.
Definition: task_mod.f90:4