62 integer,
parameter:: nt=4
63 integer,
parameter:: particle_kind=8
65 real:: per_cell, per_mass
67 integer:: id, it, time_order, n_trade
69 real,
allocatable:: v(:,:,:,:)
70 real(8),
allocatable:: wt(:)
81 procedure:: post_update
83 procedure,
private:: mpi_append_particles
84 procedure,
private:: mpi_read_particles
92 SUBROUTINE init (self, name)
94 character(len=*),
optional:: name
96 call trace%begin (
'trace_particles_t%init')
103 SUBROUTINE dealloc (self)
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')
113 END SUBROUTINE dealloc
118 SUBROUTINE add (self, link)
120 class(
link_t),
pointer:: link
122 real,
pointer:: d(:,:,:), pv(:,:,:,:)
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
132 call trace%begin (
'trace_particles_t%add')
133 self%patch => task2patch(link%task)
134 associate(patch => self%patch)
139 read (io%input, trace_particle_params, iostat=iostat)
141 write (io%output, trace_particle_params)
149 self%per_cell = per_cell
150 self%time_order = time_order
151 self%verbose = verbose
159 call self%particle_list_t%init
167 time_order = min(time_order,patch%nt)
168 allocate (self%wt(time_order))
169 select case (time_order)
171 self%wt(1:2) = [1.5, -0.5]
173 self%wt(1:3) = [23./12., -16./12., 5./12.]
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')
184 self%v(:,:,:,i) = pv(:,:,:,i)/d(:,:,:)
189 call self%particle_list_t%init (
'trace')
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)
199 if (random%ran1() < (self%per_cell-np))
then 207 call io%bits_mem (-storage_size(self%v), &
208 product(shape(self%v)),
'-tracep%v')
215 ran = random%ran1() - 0.5
223 associate(patch => self%patch)
227 p%r(:,it) = [ix+ran(),iy+ran(),iz+ran()]
228 p%v(:,it) = self%v(ix,iy,iz,:)*patch%dtime
232 call self%append (node)
244 SUBROUTINE update (self)
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
256 if (.not. self%on)
return 257 call trace%begin (
'trace_particles_t%update', itimer=itimer)
258 associate(patch => self%patch)
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))
268 self%v(:,:,:,i) = pv(:,:,:,i)/d(:,:,:)
273 if (self%import%n > 0)
then 274 node => self%import%head
275 do while (
associated(node))
277 call self%import%remove (node)
278 call self%append (node)
281 if (self%verbose > 0) &
282 write (io%output,*) self%id,
' has', self%n,
' particles' 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')
301 do while (
associated(node))
303 class is (particle_t)
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)
315 p%v(:,o1)*self%wt(1) + &
316 p%v(:,o2)*self%wt(2) + &
321 p%v(:,o1)*self%wt(1) + &
331 if (any(r < self%lf .or. r > self%uf))
then 332 call self%remove (node)
333 call self%export%append (node)
338 if (self%verbose > 0 .and. nrem>0) &
339 write (io%output,*) self%id,
' has', self%n,
' particles', nrem
340 call trace%end (itimer)
342 END SUBROUTINE update
350 SUBROUTINE interpolate (self)
354 integer:: sz(4), ip(3), i, j
356 class(dll_node_t),
pointer:: node
357 class(particle_t),
pointer:: part
358 integer,
save:: itimer=0
365 call trace%begin (
'trace_particles_t%interpolate', itimer=itimer)
368 do while (
associated(node))
374 p = 1.0 + part%r(:,it)
376 ip = max(1,min(ip,sz(1:3)))
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)))
391 call trace%end (itimer)
392 END SUBROUTINE interpolate
397 SUBROUTINE trade (self, nbor)
400 class(dll_node_t),
pointer:: node, next
401 class(particle_t),
pointer:: part
402 real(kind=particle_kind):: r(3)
405 if (self%export%n == 0)
return 406 node => self%export%head
407 do while (
associated(node))
410 class is (particle_t)
413 r = part%r(:,part%it) + nbor%offset
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
433 SUBROUTINE post_update (self)
436 class(link_t),
pointer:: nbor
437 class(task_t),
pointer:: nbtask
438 class(*),
pointer:: nbpart
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))
456 nbpart => nbtask%connect%trace_particles
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)
468 END SUBROUTINE post_update
473 SUBROUTINE pack (self)
475 integer,
pointer:: buffer(:)
476 real,
allocatable:: r(:,:,:)
477 class(dll_node_t),
pointer:: part
480 call trace%begin (
'trace_particles_t%pack')
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)
496 SUBROUTINE unpack (self)
498 integer,
pointer:: buffer(:)
499 real,
allocatable:: r(:,:,:)
500 class(dll_node_t),
pointer:: part
503 call trace%begin (
'trace_particles_t%pack')
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)
514 END SUBROUTINE unpack
519 SUBROUTINE mpi_append_particles (self, particles)
521 class(particle_list_t):: particles
522 real,
allocatable:: r(:,:,:)
523 class(dll_node_t),
pointer:: part
526 call trace%begin (
'mpi_buffer%append_particles')
530 allocate (r(3,nt,self%n))
533 do while (
associated(part))
535 class is (particle_t)
544 call self%mpi_buffer%append (r)
547 END SUBROUTINE mpi_append_particles
552 SUBROUTINE mpi_read_particles (self, particles)
554 class(particle_list_t):: particles
555 real,
allocatable:: r(:,:,:)
556 class(dll_node_t),
pointer:: part
559 call trace%begin (
'mpi_buffer%read_particles')
563 allocate (r(3,nt,self%n))
564 call self%mpi_buffer%read (r)
567 do while (
associated(part))
569 class is (particle_t)
580 END SUBROUTINE mpi_read_particles
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...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Doubly linked list (DLL), carrying anything, as simply as possible.
Module with list handling for generic class task_t objects.
Particle list, extending a doubly-linked list. Each particle maintains arrays with previous positions...
Template module for tasks.