DISPATCH
particle_list_mod.f90
1 !===============================================================================
2 !> Particle list, extending a doubly-linked list. Each particle maintains
3 !> arrays with previous positions, velocities, and times, which are brought
4 !> along if/when it changes owner patch or rank. The memory footprint of a
5 !> particle is 12 words for position, 12 words for velocity, and 8 words for
6 !> time, all together 32 words = 128 bytes (plus a few words for id and weight).
7 !> This could be reduced to half, by keeping only two previous positions and
8 !> velocities in the data type.
9 !===============================================================================
11  USE io_unit_mod
12  USE dll_mod
13  USE omp_timer_mod
14  USE particle_mod
15  USE patch_mod
16  implicit none
17  private
18  !-----------------------------------------------------------------------------
19  type, public, extends(dll_t):: particle_list_t
20  integer:: verbose=0
21  integer:: nt=4
22  contains
23  procedure:: init
24  procedure:: force
25  procedure:: force_field
26  procedure:: test
27  procedure:: print
28  end type
29  !-----------------------------------------------------------------------------
30 CONTAINS
31 
32 !===============================================================================
33 !> Initialize a particle list
34 !===============================================================================
35 SUBROUTINE init (self, name)
36  class(particle_list_t):: self
37  character(len=*), optional:: name
38  !-----------------------------------------------------------------------------
39  call self%dll_t%init (name)
40  call particle%init
41 END SUBROUTINE init
42 
43 !===============================================================================
44 !> Initialize a particle list
45 !===============================================================================
46 SUBROUTINE append (self, new)
47  class(particle_list_t):: self
48  class(particle_t), pointer:: new
49  class(dll_node_t), pointer:: new_dll
50  !-----------------------------------------------------------------------------
51  new_dll => new
52  call self%dll_t%append (new_dll)
53 END SUBROUTINE append
54 
55 !===============================================================================
56 !> Compute force from other particles by direct summation
57 !===============================================================================
58 FUNCTION force (self, p)
59  class(particle_list_t):: self
60  class(particle_t), pointer:: p
61  real(8):: force(3)
62  !.............................................................................
63  class(dll_node_t), pointer:: o
64  integer:: i, i0, i1, j
65  real(8):: rp(3), pt, qt, r(3), r1, r2, r3
66  !-----------------------------------------------------------------------------
67  force = 0.0_8
68  rp = p%r(:,p%it)
69  o => self%head
70  do while (associated(o))
71  if (.not.associated (o, p)) then
72  select type (o)
73  class is (particle_t)
74  do i=2,self%nt-1
75  i1 = o%iit(i)
76  if (o%t(i1) > p%time) exit
77  end do
78  i0 = o%iit(i-1)
79  pt = (p%time-o%t(i0))/max(o%t(i1)-o%t(i0),tiny(1d0))
80  qt = 1d0-pt
81  r2 = tiny(1d0)
82  do j=1,3
83  r(j) = (qt*o%r(j,i0) + pt*o%r(j,i1)) - rp(j)
84  r2 = r2 + r(j)**2
85  end do
86  r1 = 1.0/sqrt(r2)
87  r3 = r1**3
88  force = force - o%mass*r*r3
89  end select
90  end if
91  o => o%next
92  end do
93 END FUNCTION force
94 
95 !===============================================================================
96 !===============================================================================
97 SUBROUTINE force_field (self, patch)
98  class(particle_list_t):: self
99  class(patch_t):: patch
100  !.............................................................................
101  class(dll_node_t), pointer:: p
102  real(8), allocatable:: rp(:,:), m(:)
103  real(8):: force(3), pt, qt, rc(3), r(3), r1, r3
104  integer:: i, i0, i1, ip, ix, iy, iz
105  !-----------------------------------------------------------------------------
106  ! Cache all particle position into array, in order to be able to loop over
107  ! particles
108  !-----------------------------------------------------------------------------
109  allocate (rp(self%n,3), m(self%n))
110  ip = 0
111  p => self%head
112  do while (associated(p))
113  ip = ip+1
114  select type (p)
115  class is (particle_t)
116  do i=2,self%nt-1
117  i0 = p%iit(i-1)
118  i1 = p%iit(i)
119  if (p%t(i1) > patch%time) exit
120  end do
121  pt = (patch%time - p%t(i0))/max(p%t(i1)-p%t(i0),tiny(1d0))
122  qt = 1d0-pt
123  rp(:,ip) = qt*p%r(:,i0) + pt*p%r(:,i1)
124  m(ip) = p%mass
125  end select
126  p => p%next
127  end do
128  !-----------------------------------------------------------------------------
129  ! For each cell, loop over particles and accumulate force in double precision
130  !-----------------------------------------------------------------------------
131  associate(m1=>patch%mesh(1), m2=>patch%mesh(2), m3=>patch%mesh(3))
132  do iz=m3%li,m3%ui
133  rc(3) = m3%p + m3%r(iz)
134  do iy=m2%li,m2%ui
135  rc(2) = m2%p + m2%r(iy)
136  do ix=m1%li,m1%ui
137  rc(1) = m1%p + m1%r(ix)
138  force = 0d0
139  do ip=1,self%n
140  r = rp(:,ip) - rc
141  r1 = 1d0/sqrt(r(1)**2 + r(2)**2 + r(3)**2 + tiny(1d0))
142  r3 = r1*r1*r1
143  force = force - m(ip)*r*r3
144  end do
145  patch%force_per_unit_mass(ix,iy,iz,:) = &
146  patch%force_per_unit_mass(ix,iy,iz,:) + force
147  end do
148  end do
149  end do
150  end associate
151  deallocate (rp)
152 END SUBROUTINE force_field
153 
154 !===============================================================================
155 !> Test particle list functionality
156 !===============================================================================
157 SUBROUTINE test (self)
158  class(particle_list_t):: self
159  integer:: id, n
160  class(particle_t), pointer:: p
161  class(dll_node_t), pointer :: node, next, new
162  real(8)::used
163  logical:: flip
164  !-----------------------------------------------------------------------------
165  call self%init
166  !-----------------------------------------------------------------------------
167  ! Make a list with 100000 particles
168  !-----------------------------------------------------------------------------
169  n = 100000
170  used = wallclock()
171  do id=1,n
172  allocate (p)
173  call p%init
174  p%r = 0.5_8
175  p%ds = 1.0_8
176  p%v = 1.0_8
177  node => p
178  call self%append (node)
179  end do
180  used = wallclock()-used
181  print *, 'particle list n =', self%n, used/self%n
182  !-----------------------------------------------------------------------------
183  ! Remove every 2nd particle
184  !-----------------------------------------------------------------------------
185  flip = .false.
186  used = wallclock()
187  node => self%head
188  do while (associated(node))
189  next => node%next
190  select type (node)
191  type is (particle_t)
192  p => node
193  end select
194  if (flip) then
195  call self%remove (node)
196  deallocate (node)
197  end if
198  flip = .not.flip
199  node => next
200  end do
201  used = wallclock()-used
202  print *, 'particle list n =', self%n, used/self%n
203  !-----------------------------------------------------------------------------
204  ! Add one new particle for every 2nd existing one
205  !-----------------------------------------------------------------------------
206  flip = .false.
207  used = wallclock()
208  node => self%head
209  do while (associated(node))
210  next => node%next
211  allocate (p)
212  new => p
213  if (flip) then
214  call self%insert_before (node, new)
215  end if
216  flip = .not.flip
217  node => next
218  end do
219  used = wallclock()-used
220  print *, 'particle list n =', self%n, used/self%n
221 END SUBROUTINE test
222 
223 !===============================================================================
224 !> Print a particle list
225 !===============================================================================
226 SUBROUTINE print (self)
227  class(particle_list_t):: self
228  class(dll_node_t), pointer:: p
229  p => self%head
230  do while (associated(p))
231  select type (p)
232  class is (particle_t)
233  print *, p%id, p%r(:,p%it)
234  end select
235  p => p%next
236  end do
237 END SUBROUTINE print
238 
239 END MODULE particle_list_mod
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
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
Particle list, extending a doubly-linked list. Each particle maintains arrays with previous positions...