File: N:\mfix\model\des\des_cluster_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: DES_CLUSTER                                            !
4     !                                                                      !
5     !  Purpose: Common elements for MFIX-DEM cluster identification        !
6     !  condition.                                                          !
7     !                                                                      !
8     !  Author: J.Galvin, J.Musser                         Date:  Nov-12    !
9     !  Modified: S. Benyahia                              Date:  Dec-12    !
10     !                                                                      !
11     !  Comments: Info on clusters such as average void fraction, Re and    !
12     !  Comments: cluster size can now be printed from this file            !
13     !                                                                      !
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
15           MODULE DES_CLUSTER
16     
17     
18     !-----------------------------------------------
19     ! Modules
20     !-----------------------------------------------
21           USE param
22           USE param1
23           USE compar
24           USE fldvar
25           USE physprop
26           USE discretelement
27     
28           use mpi_utility
29           use desmpi_wrapper
30           use desmpi
31           use compar
32           use desmpi
33           use parallel
34           use sendrecv
35     
36     
37           IMPLICIT NONE
38     !-----------------------------------------------
39     
40     
41     ! Number of clusters.
42           INTEGER ClusterCount
43     
44     ! Number of base particles held in current search list
45           INTEGER pSearchHistoryCount
46     
47     !----------------------------------------------->>>
48           TYPE PARTICLE_TYPE
49     ! Global index of particle.
50              INTEGER :: ID
51     ! Successor in particle list
52              TYPE (PARTICLE_TYPE), POINTER :: next_particle
53           END TYPE PARTICLE_TYPE
54     !-----------------------------------------------<<<
55     
56     ! Basic particle link-list.
57           TYPE(PARTICLE_TYPE), POINTER :: PSEARCH_HISTORY_LL
58     
59     !----------------------------------------------->>>
60           TYPE CLUSTER_TYPE
61     ! Cluster index.
62              INTEGER :: ID
63     ! Number of particles in cluster.
64              INTEGER :: ParticleCount
65     ! Successor in cluster link-list.
66              TYPE (CLUSTER_TYPE), POINTER :: next_cluster
67     ! Particle link-list
68              TYPE (PARTICLE_TYPE), POINTER :: PARTICLE_LL
69           END TYPE CLUSTER_TYPE
70     !-----------------------------------------------<<<
71     
72     ! Cluster link-list.
73           TYPE(CLUSTER_TYPE), POINTER :: CLUSTER_LL
74     
75     !<<<-------------------  DMP Related Variables  -------------------->>>!
76     
77     ! Process that does global cluster calculations.
78           integer, parameter :: clusterPE = 0
79     
80           TYPE pType
81              INTEGER :: map
82              INTEGER :: id
83     ! Successor in particle list
84              TYPE (pType), POINTER :: next
85           END TYPE pType
86     
87     
88           TYPE cType
89     ! Number of particles in cluster.
90              INTEGER :: size
91     ! Particle link-list
92              TYPE (pType), POINTER :: particle
93           END TYPE cType
94     
95     ! Cluster array.
96           TYPE(cType), dimension(:), allocatable, target :: clusters
97           INTEGER clusterCount_all
98     
99     ! Send/Recv information:
100     !----------------------->>
101     
102     ! Number of data entrires being sent from the local process to Root
103     ! during a call to cluster_gather.
104           INTEGER send_cnt
105     
106     ! Number of data entries being received by root from each process
107     ! during a call to cluster_gather.
108           INTEGER, dimension(:), allocatable :: recv_cnt
109     
110     ! Displacement of data from all processes during a cluster_gather.
111           INTEGER, dimension(:), allocatable :: recv_dsp
112     
113     ! Size of array needed by Root to accept all data during a gather call.
114           INTEGER recv_sum
115     
116     
117     ! Module interfaces:
118     !----------------------->>
119     
120           interface getClusterParticleData
121              module procedure getClusterParticleData_1i
122              module procedure getClusterParticleData_2i
123              module procedure getClusterParticleData_1d
124              module procedure getClusterParticleData_2d
125           end interface
126     
127           interface getClusterFieldData
128              module procedure getClusterFieldData_1d
129              module procedure getClusterFieldData_3d
130           end interface
131     
132           interface sendClusterData
133              module procedure sendClusterData_1d
134           end interface
135     
136     
137           contains
138     
139     !......................................................................!
140     !                                                                      !
141     !......................................................................!
142           SUBROUTINE CREATE_CLUSTER(cluster)
143     
144           TYPE(CLUSTER_TYPE), INTENT(OUT), POINTER :: cluster
145     
146           ALLOCATE(cluster)
147     
148           NULLIFY(cluster%next_cluster)
149           NULLIFY(cluster%PARTICLE_LL)
150     
151           if(ClusterCount == 0) then
152     ! no clusters have been created/identified
153              if(associated(CLUSTER_LL)) then
154                 print*, ' Error - cluster pointer already associated!'
155                 CALL MFIX_EXIT(myPE)
156              else
157                 ClusterCount = 1
158                 Cluster%ParticleCount = 0
159                 cluster%ID = ClusterCount
160     ! create first cluster of linked list of clusters. with cluster_ll
161     ! always being the 'first' in the list
162                 CLUSTER_LL => cluster
163              endif
164           else
165              if(.NOT.associated(CLUSTER_LL)) then
166                 print*, ' Error - cluster pointer is not associated!'
167                 CALL MFIX_EXIT(myPE)
168              else
169                 ClusterCount = ClusterCount + 1
170                 Cluster%ParticleCount = 0
171                 cluster%ID = ClusterCount
172     ! establish the link between the new cluster and existing cluster list
173                 cluster%next_cluster => CLUSTER_LL
174     ! reassign/point cluster_ll to be 'first' in the linked list
175                 CLUSTER_LL => cluster
176     ! as a result the linked list of clusters is created from 'bottom-up'
177     ! with the new cluster always being inserted before any existing
178     ! clusters
179              endif
180           ENDIF
181     
182           END SUBROUTINE CREATE_CLUSTER
183     
184     
185     !......................................................................!
186     !                                                                      !
187     !......................................................................!
188           SUBROUTINE GetTopCluster(cluster)
189     
190           TYPE(CLUSTER_TYPE), INTENT(INOUT), POINTER :: cluster
191     
192           NULLIFY(cluster)
193     
194     ! first check that a list has been formed
195           if(.NOT.associated(CLUSTER_LL) .AND. ClusterCount == 0) then
196              write(*,"(//,3x,A,//)") 'No clusters to delete!'
197           elseif(.NOT.associated(CLUSTER_LL) .AND. ClusterCount /= 0) then
198              write(*,"(//,3x,A,//)") ' Error with ClusterCount and pointer - Er:1'
199           elseif(associated(CLUSTER_LL) .AND. ClusterCount == 0) then
200              write(*,"(//,3x,A,//)") ' Error with ClusterCount and pointer - Er:2'
201           else
202     ! now identify top cluster
203              CALL getNextCluster(cluster)
204           endif
205           END SUBROUTINE GetTopCluster
206     
207     
208     !......................................................................!
209     !                                                                      !
210     !......................................................................!
211           SUBROUTINE DeleteTopCluster(cluster)
212     
213           TYPE(CLUSTER_TYPE), INTENT(INOUT), POINTER :: cluster
214     
215     ! dbg
216     !      write(*,"(/A,I7)")' Deleting top cluster: ',cluster%ID
217           CALL DELETE_PARTICLES_IN_CLUSTER(cluster)
218     
219           if(associated(cluster%next_cluster)) then
220              CLUSTER_LL => cluster%next_cluster
221           else
222              NULLIFY(CLUSTER_LL)
223           endif
224     
225           ClusterCount = ClusterCount - 1
226           NULLIFY(cluster%next_cluster)
227           NULLIFY(cluster%PARTICLE_LL)
228           DEALLOCATE(cluster)
229           NULLIFY(cluster)
230     
231           if(ClusterCount <0) then
232              write(*,"()")' ClusterCount < 0'
233              CALL MFIX_EXIT(myPE)
234           endif
235     
236           END SUBROUTINE DeleteTopCluster
237     
238     
239     !......................................................................!
240     !                                                                      !
241     !......................................................................!
242           SUBROUTINE DELETE_CLUSTERS()
243     
244           TYPE(CLUSTER_TYPE), POINTER :: cluster
245           INTEGER cL, cDeleted
246     
247           NULLIFY(cluster)
248           cDeleted = 0
249     
250           if(.NOT.associated(CLUSTER_LL) .AND. ClusterCount == 0) then
251     !         write(*,"(//,3x,A,//)") 'No clusters to delete!'
252           elseif(.NOT.associated(CLUSTER_LL) .AND. ClusterCount /= 0) then
253              write(*,"(//,3x,A,//)") ' Error with ClusterCount and pointer - Er:1'
254           elseif(associated(CLUSTER_LL) .AND. ClusterCount == 0) then
255              write(*,"(//,3x,A,//)") ' Error with ClusterCount and pointer - Er:2'
256           else
257              do cL =1, ClusterCount
258                 CALL getNextCluster(cluster)
259     ! dbg
260     !            write(*,"(/A,I7)")' Deleting cluster: ',cluster%ID
261                 CALL DELETE_PARTICLES_IN_CLUSTER(cluster)
262     
263                 if(associated(cluster%next_cluster)) then
264                    CLUSTER_LL => cluster%next_cluster
265                 else
266                    NULLIFY(CLUSTER_LL)
267                 endif
268                 cDeleted = cDeleted + 1
269                 NULLIFY(cluster%next_cluster)
270                 NULLIFY(cluster%PARTICLE_LL)
271                 DEALLOCATE(cluster)
272                 NULLIFY(cluster)
273              enddo
274     
275              if(cDeleted == ClusterCount) then
276     ! dbg
277     !            write(*,"(4X,A,I7)")'Number of clusters deleted: ',cDeleted
278                 ClusterCount = 0
279              else
280                 write(*,"()")' cDeleted /= ClusterCount'
281                 CALL MFIX_EXIT(myPE)
282              endif
283     
284           endif
285     
286           END SUBROUTINE DELETE_CLUSTERS
287     
288     
289     
290     !......................................................................!
291     !                                                                      !
292     !......................................................................!
293           SUBROUTINE DELETE_PARTICLES_IN_CLUSTER(cluster)
294     
295           TYPE(CLUSTER_TYPE), INTENT(INOUT), POINTER :: cluster
296     
297           TYPE(PARTICLE_TYPE), POINTER :: particle
298           INTEGER pL, pDeleted
299     
300           NULLIFY(particle)
301           pDeleted = 0
302     
303           if(.NOT.associated(cluster%PARTICLE_LL) .AND. cluster%ParticleCount == 0) then
304              write(*,"(//,3x,A,//)") 'No particles to delete!'
305           elseif(.NOT.associated(cluster%PARTICLE_LL) .AND. cluster%ParticleCount /= 0) then
306              write(*,"(//,3x,A,//)") ' Error with ParticleCount and pointer - Er:1'
307           elseif(associated(cluster%PARTICLE_LL) .AND. cluster%ParticleCount == 0) then
308              write(*,"(//,3x,A,//)") ' Error with ParticleCount and pointer - Er:2'
309           else
310              do pL =1, cluster%ParticleCount
311                 CALL GetNextParticle(cluster, particle)
312     ! dbg
313     !            write(*,"(6X,A,I7)")' Deleting particle: ',particle%ID
314     
315                 if(associated(particle%next_particle)) then
316                    cluster%PARTICLE_LL => particle%next_particle
317                 else
318                    NULLIFY(cluster%PARTICLE_LL)
319                 endif
320                 pDeleted = pDeleted + 1
321                 NULLIFY(particle%next_particle)
322                 DEALLOCATE(particle)
323                 NULLIFY(particle)
324              enddo
325     
326              if(pDeleted == cluster%ParticleCount) then
327     ! dbg
328     !            if(pDeleted >1 ) &
329     !               write(*,"(4X,A,I7)")'Number of deleted particles: ',&
330     !               pDeleted
331                 cluster%ParticleCount = 0
332              else
333                 write(*,"()")' pDeleted /= cluster%ParticleCount'
334                 CALL MFIX_EXIT(myPE)
335              endif
336     
337           endif
338     
339     
340           END SUBROUTINE DELETE_PARTICLES_IN_CLUSTER
341     
342     
343     
344     !......................................................................!
345     !                                                                      !
346     !......................................................................!
347           SUBROUTINE ADD_PARTICLE_TO_CLUSTER(cluster, pID)
348     
349           TYPE(CLUSTER_TYPE), INTENT(INOUT), POINTER :: cluster
350           INTEGER, INTENT(IN) :: pID
351     
352           TYPE (PARTICLE_TYPE), POINTER :: particle
353     
354     ! Create a particle instance.
355           ALLOCATE(particle)
356           NULLIFY(particle%next_particle)
357     ! Set the particle id.
358           particle%ID = pID
359     ! Set the cluster pointer to particle.
360           IF(cluster%ParticleCount == 0) THEN
361              cluster%ParticleCount = 1
362              cluster%PARTICLE_LL => particle
363           ELSE
364              cluster%ParticleCount = cluster%ParticleCount + 1
365              particle%next_particle => cluster%PARTICLE_LL
366              cluster%PARTICLE_LL => particle
367           ENDIF
368     
369     
370           END SUBROUTINE ADD_PARTICLE_TO_CLUSTER
371     
372     
373     !......................................................................!
374     !                                                                      !
375     !......................................................................!
376           SUBROUTINE getNextCluster(cluster)
377     
378           TYPE(CLUSTER_TYPE), INTENT(INOUT), POINTER :: cluster
379     
380           if(.NOT.associated(cluster)) then
381     ! point cluster to first cluster in linked list of clusters
382              cluster => CLUSTER_LL
383           elseif(associated(cluster%next_cluster)) then
384              cluster => cluster%next_cluster
385           else
386              print*,' You are looking for a cluster that does not exist'
387              CALL MFIX_EXIT(myPE)
388           endif
389           END SUBROUTINE getNextCluster
390     
391     
392     
393     !......................................................................!
394     !                                                                      !
395     !......................................................................!
396           SUBROUTINE GetNextParticle(cluster, particle)
397     
398           TYPE(CLUSTER_TYPE), INTENT(IN), POINTER :: cluster
399           TYPE(PARTICLE_TYPE), INTENT(INOUT), POINTER :: particle
400     
401           if(.NOT.associated(particle)) then
402              particle => cluster%PARTICLE_LL
403           elseif(associated(particle%next_particle)) then
404              particle => particle%next_particle
405           else
406              print*,' You are looking for a particle that does not exist'
407              CALL MFIX_EXIT(myPE)
408           endif
409           END SUBROUTINE GetNextParticle
410     
411     
412     
413     !......................................................................!
414     !                                                                      !
415     !......................................................................!
416           SUBROUTINE ADD_PARTICLE_TO_PSEARCHHISTORY(particle, pID)
417     
418           TYPE(PARTICLE_TYPE), INTENT(OUT), POINTER :: particle
419           INTEGER, INTENT(IN) :: pID
420     
421           ALLOCATE(particle)
422     
423           NULLIFY(particle%next_particle)
424     
425     !dbg
426     !      print *, '----------> addparticle_to_psearchhistory'
427     
428           if(pSearchHistoryCount == 0) then
429     ! no particles in search list have been created/identified
430              if(associated(PSEARCH_HISTORY_LL)) then
431                      print*, ' Error - particle history pointer already &
432                         & associated!'
433                 CALL MFIX_EXIT(myPE)
434              else
435                 pSearchHistoryCount = 1
436                 particle%ID = pID
437     ! create first particle of linked list of particles. with
438     ! psearch_history_ll always pointing to the 'first' in the list
439                 PSEARCH_HISTORY_LL => particle
440     ! dbg
441     !            write(*,"(/A,I7)") 'Create history w/ particle ', pID
442              endif
443           else
444              if(.NOT.associated(PSEARCH_HISTORY_LL)) then
445                 print*, ' Error - particle history  pointer is not &
446                    & associated!'
447                 CALL MFIX_EXIT(myPE)
448              else
449                 pSearchHistoryCount = pSearchHistoryCount + 1
450                 particle%ID = pID
451     ! establish the link between the new particle and existing particle list
452                 particle%next_particle => PSEARCH_HISTORY_LL
453     ! reassign/point psearch_history_ll to be 'first' in the linked list
454                 PSEARCH_HISTORY_LL => particle
455     ! as a result the linked list of particles is created from 'bottom-up'
456     ! with the new particle  always being inserted before any existing
457     ! particles
458     ! dbg
459     !            print *, 'Add particle in history ', pID
460     
461              endif
462           ENDIF
463     
464           END SUBROUTINE ADD_PARTICLE_TO_PSEARCHHISTORY
465     
466     
467     
468     !......................................................................!
469     !                                                                      !
470     !......................................................................!
471           SUBROUTINE GetTopParticle_in_PSearchHistory(particle)
472     
473           TYPE(PARTICLE_TYPE), INTENT(INOUT), POINTER :: particle
474     
475     !dbg
476     !     print *, '----------> gettopparticle_in_psearchhistory'
477     
478           if(.NOT.associated(particle)) then
479     ! point particle to first particle in linked list of particles
480              particle => PSEARCH_HISTORY_LL
481           elseif(associated(particle%next_particle)) then
482     ! this will never happen if we nullify particle at start of this
483     ! routine.
484              particle => particle%next_particle
485           else
486              print*,' You are looking for a particle that does not exist'
487              CALL MFIX_EXIT(myPE)
488           endif
489           END SUBROUTINE GetTopParticle_in_PSearchHistory
490     
491     
492     
493     !......................................................................!
494     !                                                                      !
495     !......................................................................!
496           SUBROUTINE DeleteTopParticle_in_PSearchHistory()
497     
498           TYPE(PARTICLE_TYPE), POINTER :: particle
499     
500           NULLIFY(particle)
501     !dbg
502     !      print *, '----------> deletetopparticle_in_psearchhistory'
503     
504           if(.NOT.associated(PSEARCH_HISTORY_LL) .AND. &
505              pSearchHistoryCount == 0) then
506              write(*,"(//,3x,A,//)") 'No particles in history to delete!'
507           elseif(.NOT.associated(PSEARCH_HISTORY_LL) .AND. &
508              pSearchHistoryCount /= 0) then
509              write(*,"(//,3x,A,//)") &
510                 ' Error with pSearchHistoryCount and pointer - Er:1'
511           elseif(associated(PSEARCH_HISTORY_LL) .AND. &
512              pSearchHistoryCount == 0) then
513              write(*,"(//,3x,A,//)") &
514                 ' Error with pSearchHistoryCount and pointer - Er:2'
515           else
516              CALL GetTopParticle_In_PSearchHistory(particle)
517     ! dbg
518     !         write(*,"(A,I7)")' Deleting particle: ',particle%ID
519     
520              if(associated(particle%next_particle)) then
521                 PSEARCH_HISTORY_LL => particle%next_particle
522              else
523                 NULLIFY(PSEARCH_HISTORY_LL)
524              endif
525              pSearchHistoryCount = pSearchHistoryCount - 1
526              NULLIFY(particle%next_particle)
527              DEALLOCATE(particle)
528              NULLIFY(particle)
529     
530              if(pSearchHistoryCount < 0) then
531                 write(*,"(4X,A)")'pSearchHistoryCount < 0'
532                 CALL MFIX_EXIT(myPE)
533              endif
534     
535           endif
536     
537           END SUBROUTINE DeleteTopParticle_In_PSearchHistory
538     
539     
540     
541     !......................................................................!
542     !                                                                      !
543     !......................................................................!
544           SUBROUTINE DELETE_PSEARCHHISTORY()
545     
546           TYPE(PARTICLE_TYPE), POINTER :: particle
547           INTEGER pL, pDeleted
548     
549           NULLIFY(particle)
550           pDeleted = 0
551     !dbg
552     !      print *, '----------> delete_psearch_history'
553     
554           if(.NOT.associated(PSEARCH_HISTORY_LL) .AND. &
555              pSearchHistoryCount == 0) then
556     ! dbg
557     !         write(*,"(//,3x,A,//)") 'No particles in history to delete!'
558           elseif(.NOT.associated(PSEARCH_HISTORY_LL) .AND. &
559              pSearchHistoryCount /= 0) then
560              write(*,"(//,3x,A,//)") &
561                 ' Error with pSearchHistoryCount and pointer - Er:1'
562           elseif(associated(PSEARCH_HISTORY_LL) .AND. &
563              pSearchHistoryCount == 0) then
564              write(*,"(//,3x,A,//)") &
565                 ' Error with pSearchHistoryCount and pointer - Er:2'
566           else
567              do pL =1, pSearchHistoryCount
568     
569                 CALL GetTopParticle_In_PSearchHistory(particle)
570     ! dbg
571     !            write(*,"(/A,I7)")' Deleting particle: ',particle%ID
572     
573                 if(associated(particle%next_particle)) then
574                    PSEARCH_HISTORY_LL => particle%next_particle
575                 else
576                    NULLIFY(PSEARCH_HISTORY_LL)
577                 endif
578                 pDeleted = pDeleted + 1
579                 NULLIFY(particle%next_particle)
580                 DEALLOCATE(particle)
581                 NULLIFY(particle)
582              enddo
583     
584              if(pDeleted == pSearchHistoryCount) then
585     ! dbg
586     !            write(*,"(4X,A,I7)")'Number of particles deleted: ',pDeleted
587                 pSearchHistoryCount = 0
588              else
589                 write(*,"()")' pDeleted /= PSearchHistoryCount'
590                 CALL MFIX_EXIT(myPE)
591              endif
592     
593           endif
594     
595           END SUBROUTINE DELETE_PSEARCHHISTORY
596     
597     
598     
599     !......................................................................!
600     !                                                                      !
601     !......................................................................!
602           SUBROUTINE PRINT_CLUSTERS
603     
604           use run
605     
606           implicit none
607     
608     ! Generic loop counters
609           integer lc1, lc2, lc3, lc4
610     ! Particle index (mapped to cluster send/recv)
611           integer L
612     
613     ! Flags used for debugging.
614     !   0: No messages
615     !   1: Screen only messages
616     !   2: Detailed Files are written by each process
617           INTEGER, parameter :: dbg_level = 0
618     
619     ! Particle properties (for cluster bound particles only)
620           double precision, dimension(:),   allocatable :: lRad
621           double precision, dimension(:,:), allocatable :: lPos
622           double precision, dimension(:,:), allocatable :: lVel_s
623     
624     ! File properticles (for cluster bound particles only)
625           double precision, dimension(:),   allocatable :: lEpg
626           double precision, dimension(:),   allocatable :: lRog
627           double precision, dimension(:),   allocatable :: lMug
628           double precision, dimension(:,:), allocatable :: lVel_g
629     
630     ! An array for post process data storage.
631           double precision, dimension(:), allocatable :: lPost
632     
633     ! Cluster analysis variables:
634           integer lMin(1:3), lMax(1:3)
635     
636           double precision avgEpg, avgRe, avgSlip, lSlip
637           double precision avgVel_s(1:3), avgVel_g(1:3)
638     
639           double precision posMin(1:3), posMax(1:3)
640           double precision clSize(1:3), clDiameter
641     
642           double precision cSize
643     
644           Type(cType), pointer :: cThis
645           Type(pType), pointer :: pThis
646     
647     ! Merge cluster data from all processes and set up send/recv datav.
648           call init_print_clusters(dbg_level)
649     
650     ! Get particle data from all processes. Only collect what is needed.
651           call getClusterParticleData(DES_RADIUS,  lRad)   ! Radius
652           call getClusterParticleData(DES_POS_NEW, lPos)   ! Position
653           call getClusterParticleData(DES_VEL_NEW, lVel_s) ! Velocity
654     
655     ! Get field variable data from all processes.
656           call getClusterFieldData(Ep_g, lEpg)
657           call getClusterFieldData(RO_G, lRog)
658           call getClusterFieldData(Mu_g, lMug)
659           call getClusterFieldData(U_g, V_g, W_g, lVel_g)
660     
661     ! Allocate output data variable. (Root process only.) This data must
662     ! be sent back to the individual processes.
663           allocate(lPost( recv_sum)); lPost = zero
664     ! Initialize the global storage variable.
665           postCluster = zero
666     
667     ! Calculate and report cluster stats.
668           if(myPE == clusterPE) then
669              open (unit=203, file='clusterInfo.dat', &
670                 status='unknown', position='append')
671              if(clusterCount_all == 0) then
672                 write(203,"(/3X,'Time: ',F9.6,3x,'No clusters to print.')")&
673                    time
674              else
675                 write(203,"(/3X,'Time: ',F9.6)") time
676     
677                 nullify(cThis)
678                 nullify(pThis)
679                 lc2 = 0
680                 lp_lc10: do lc1 = 1, clusterCount_all
681     ! Only clusters with more than 3 particles are processed. Keep track of
682     ! the number reported vs. number processed.
683                    if(clusters(lc1)%size < 4) cycle lp_lc10
684                    lc2 = lc2 + 1
685     
686     ! Initialize variables.
687                    avgEpg = zero
688                    avgVel_g = zero
689                    avgVel_s = zero
690                    avgRe = zero
691     
692                    posMin = large_number;  lMin = 0
693                    posMax = zero;          lMax = 0
694     
695     ! Loop over particles comprising this cluster.
696                    lc3 = 0
697                    cThis => clusters(lc1)
698                    if(associated(cThis%particle)) then
699                       pThis => cThis%particle
700                    else
701                       nullify(pThis)
702                    endif
703     ! Although this loop could be based on the cluster size, it is safer
704     ! to loop through the link-list. If there was an allocation error, this
705     ! should prevent the code from crashing.
706                    do while(associated(pThis))
707                       lc3 = lc3 + 1
708     ! Get the particle map index. This index is used to get data from the
709     ! send/recv variables.
710                       L = pThis%map
711     
712     ! Report a problem if viscosity is zero. (sanity check)
713                       if(lMug(L) == zero) then
714                          write(*,"(3x,'Invalid Mu_g. ', &
715                             &'Omitting cluster data: ',I4)") lc1
716                          lc2 = lc2 - 1
717                          cycle lp_lc10
718                       endif
719     
720                       lSlip = zero
721                       do lc4=1, merge(2, 3, NO_K)
722     ! Calculate the cluster's dimensions:
723                          if((lPos(L,lc4)-lRad(L)) < posMin(lc4)) then
724                             posMin(lc4) = lPos(L,lc4)-lRad(L)
725                             lMin(lc4) = L
726                          endif
727                          if((lPos(L,lc4)+lRad(L)) > posMax(lc4)) then
728                             posMax(lc4) = lPos(L,lc4)+lRad(L)
729                             lMax(lc4) = L
730                          endif
731     ! Calculate the slip velocity.
732                          lSlip = lSlip + (lVel_g(L,lc4)-lVel_s(L,lc4))**2
733                       enddo
734                       lSlip = dsqrt(lSlip)
735     
736     ! This is a sanity check on the data. Although a zero slip velocity is
737     ! possible, (hopefully) there is always some small variance.
738                       if(lSlip == zero) then
739                          write(*,"(3x,'Invalid lSlip. ', &
740                             &'Omitting cluster data: ',I4)") lc1
741                          lc2 = lc2 - 1
742                          cycle lp_lc10
743                       endif
744     
745                       avgVel_s = avgVel_s + lVel_s(L,:)
746                       avgVel_g = avgVel_g + lVel_g(L,:)
747                       avgEpg = avgEpg + lEpg(L)
748                       avgRe = avgRe + lRog(L)*lEpg(L)*lSlip*2.0d0*lRad(L)/lMug(L)
749     
750     ! Store the cluster size with the particle. (Output data in vtk files.)
751                       lPost(L) = float(cThis%size)
752     
753     ! Increment the loop.
754                       if(associated(pThis%next)) then
755                          pThis => pThis%next
756                       else
757                          nullify(pThis)
758                       endif
759                    enddo
760     ! This is where the above loop is checked against the previously
761     ! calculated cluster size. They should match.
762                    if(lc3 /= cThis%size) then
763                       write(*,"(3x,'Error processing particles. ', &
764                          &'Omitting cluster data: ',I4)") lc1
765                       lc2 = lc2 - 1
766                       cycle lp_lc10
767                    endif
768     
769                    cSize = dble(cThis%size)
770                    avgEpg   = avgEpg   / cSize
771                    avgVel_g = avgVel_g / cSize
772                    avgVel_s = avgVel_s / cSize
773                    avgRe    = avgRe    / cSize
774     
775                    avgSlip = zero
776                    clSize = zero
777                    do lc4=1, merge(2, 3, NO_K)
778                       avgSlip = avgSlip + (avgVel_g(lc4)-avgVel_s(lc4))**2
779                       clSize(lc4) = posMax(lc4) - posMin(lc4)
780                    enddo
781                    avgSlip = dsqrt(avgSlip)
782     
783     ! Cluster diameter is calculated by using an equivalent circle diameter
784     ! of an ellipse in each plane. Then weighting this area by the slip
785     ! velocity normal to it. This in effect gives a hydrodynamic diameter
786     ! of the cluster.- 1660
787                    if(DO_K) then
788                       clDiameter = &
789                          clSize(2)*clSize(3)*(avgVel_g(1)-avgVel_s(1))**2 &
790                        + clSize(1)*clSize(3)*(avgVel_g(2)-avgVel_s(2))**2 &
791                        + clSize(1)*clSize(2)*(avgVel_g(3)-avgVel_s(3))**2
792     
793                       clDiameter = dsqrt(clDiameter) / avgSlip
794                    else
795                       clDiameter = undefined
796                    endif
797     
798                    write(203,"(4X,I7,4X,I7,4(1X,G13.6))") lc1, &
799                       cThis%size, avgEpg, avgRe, avgSlip, clDiameter
800     
801                 enddo lp_lc10
802     
803                 write(203,"(4X,'Clusters reported: ', I7, &
804                    &4x,'Clusters processed: ', I7)") clusterCount_all, lc2
805              endif
806              close(203)
807           endif
808     
809           call sendClusterData(lPost, postCluster)
810     
811           nullify(cThis)
812           nullify(pThis)
813     
814           if(allocated(lRad)) deallocate(lRad)
815           if(allocated(lPos)) deallocate(lPos)
816           if(allocated(lVel_s)) deallocate(lVel_s)
817     
818           if(allocated(lEpg)) deallocate(lEpg)
819           if(allocated(lRog)) deallocate(lRog)
820           if(allocated(lMug)) deallocate(lMug)
821           if(allocated(lVel_g)) deallocate(lVel_g)
822     
823           if(allocated(lPost)) deallocate(lPost)
824     
825           call finl_print_clusters(dbg_level)
826     
827           END SUBROUTINE PRINT_CLUSTERS
828     
829     
830     !......................................................................!
831     !  Module name: INIT_PRINT_CLUSTERS                                    !
832     !                                                                      !
833     !  Purpose: Collect cluster information from all process on root.      !
834     !   > Number of clusters on each process.                              !
835     !   > Number of particles in each cluster.                             !
836     !                                                                      !
837     !  This information is used to construct a send/receive map for        !
838     !  cluster data.                                                       !
839     !                                                                      !
840     !  Author: J.Musser                                   Date:  Dec-12    !
841     !......................................................................!
842           SUBROUTINE INIT_PRINT_CLUSTERS(dbg_level)
843     
844           use run
845     
846           implicit none
847     
848     ! Flags used for debugging.
849     !   0: No messages
850     !   1: Screen only messages
851     !   2: Detailed Files are written by each process
852           INTEGER, intent(in) :: dbg_level
853     
854     
855     ! Counter of processes
856           INTEGER proc, ierr
857     ! Generic loop counters
858           INTEGER lc1, lc2, lc3, lc4
859     
860           character(LEN=255) :: filename
861     
862     ! Local process data:
863     !----------------------->>
864     
865           TYPE(CLUSTER_TYPE), POINTER :: cluster
866           TYPE(PARTICLE_TYPE), POINTER :: particle
867     
868     ! Number of partcies in each cluster on local process.
869           INTEGER, dimension(:), allocatable :: pCnt
870     ! Number of clusters on local process.
871           INTEGER, dimension(:) :: cCnt(0:numPEs-1)
872     ! Global IDs of particles in clusters (array)
873           INTEGER, dimension(:), allocatable :: gpIDs
874     ! IS_GHOST entries for particles in clusters. (array)
875           LOGICAL, dimension(:), allocatable :: gpGhost
876     
877     ! Root process data:
878     !----------------------->>
879     
880     ! The number of clusters detected on each process.
881           INTEGER, dimension(:) :: cCnt_all(0:numPEs-1)
882     ! Sum of cCnt_all entries
883           INTEGER cCnt_sum
884     
885     ! The number of particles in each cluster.
886           INTEGER, dimension(:), allocatable :: pCnt_all
887     
888     ! Displacements used for send/receive and cluster mapping.
889           INTEGER, dimension(:) :: pCnt_dsp(0:numPEs-1)
890     
891           INTEGER, dimension(:), allocatable :: gp_dsp
892     
893     ! Global IDs of particles in clusters (array)
894           INTEGER, dimension(:), allocatable :: gpIDs_all
895     ! Ghost entries for particles in clusters. (array)
896           LOGICAL, dimension(:), allocatable :: gpGhost_all
897     
898           INTEGER, dimension(:,:), allocatable :: mergeMap
899           LOGICAL merging
900     
901           Type(cType), pointer :: this
902     
903           call des_mpi_barrier()
904     
905           if(myPE == clusterPE .and. dbg_level >=1) &
906              write(*,"(//1x,'Start data dump ',50('-'),'>'/)")
907     
908           if(dbg_level >= 2) then
909              filename=''
910              write(filename,"('dbg_check_',I2.2,'.txt')")myPE
911              open(convert='big_endian',unit=202, file=filename, &
912                 status='replace', position='append')
913              write(202,"(3x,'check A')")
914           endif
915     
916     ! Allocate receive buffer and receive count arrays.
917           allocate( recv_cnt(0:numPEs-1) )
918           allocate( recv_dsp(0:numPEs-1) )
919     
920     ! Set the gather counts array.
921           cCnt = 0; cCnt_all = 0
922           cCnt(mype) = clusterCount
923           call global_sum(cCnt ,cCnt_all)
924           call des_mpi_barrier()
925           if(dbg_level >= 1) call dbg_print_clusters(1)
926           if(dbg_level >= 2) write(202,"(5x,'check 1')")
927     
928     ! Create an array identifying the number of particles in each cluster.
929           allocate(pCnt(cCnt(myPE)))
930           NULLIFY(cluster)
931           do lc1 = 1, cCnt(myPE)
932              CALL getNextCluster(cluster)
933              pCnt(lc1) = cluster%ParticleCount
934           enddo
935           if(dbg_level >= 2) then
936              call dbg_print_clusters(2)
937              write(202,"(5x,'check 2')")
938           endif
939     
940     ! Calculate the total number of clusters in the entire domain.
941           cCnt_sum = 1
942           if(myPE == clusterPE) cCnt_sum = sum(cCnt_all)
943     
944     ! Allocate an array sized to the total number of clusters in the
945     ! entire domain.
946           allocate(pCnt_all(cCnt_sum)); pCnt_all(:) = 1
947           if(myPE == clusterPE) then
948              pCnt_dsp(0) = 0
949     ! Calculate the array displacement for MPI_GATHER
950              do proc = 1,numPEs-1
951                 pCnt_dsp(proc) = pCnt_dsp(proc-1) + cCnt_all(proc-1)
952              enddo
953           endif
954           if(dbg_level >= 2) then
955              call dbg_print_clusters(3)
956              write(202,"(5x,'check 3')")
957           endif
958     
959     ! Populate pCnt_all on clusterPE.
960           call des_mpi_gatherv(pCnt, cCnt(myPE), pCnt_all, cCnt_all,       &
961              pCnt_dsp, clusterPE, ierr)
962           if(dbg_level >= 1) call dbg_print_clusters(4)
963           if(dbg_level >= 2) write(202,"(5x,'check 4')")
964           call des_mpi_barrier()
965     
966     
967           send_cnt = sum(pCnt)
968           if(myPE == clusterPE) then
969              recv_cnt(:) = 0
970              do proc = 0, numPEs-1
971                 if(cCnt_all(proc) > 0) then
972                    do lc1=1, cCnt_all(proc)
973                       recv_cnt(proc) = recv_cnt(proc) + &
974                          pCnt_all(lc1 + pCnt_dsp(proc))
975                    enddo
976                 else
977                    recv_cnt(proc) = 0
978                 endif
979              enddo
980              recv_dsp(:) = 0
981              do proc = 1, numPEs-1
982                 recv_dsp(proc) = recv_dsp(proc-1) + recv_cnt(proc-1)
983              enddo
984              recv_sum = sum(recv_cnt)
985           else
986              recv_sum = 1
987           endif
988           if(dbg_level >= 1) call dbg_print_clusters(5)
989           if(dbg_level >= 2) write(202,"(5x,'check 5')")
990     
991     
992           allocate( gp_dsp(cCnt_sum) ); gp_dsp(:) = 0
993           if(myPE == clusterPE) then
994              do lc1 = 2, cCnt_sum
995                 gp_dsp(lc1) = gp_dsp(lc1-1) + pCnt_all(lc1-1)
996              enddo
997           endif
998           if(dbg_level >= 1) call dbg_print_clusters(6)
999           if(dbg_level >= 2) write(202,"(5x,'check 6')")
1000     
1001     
1002     ! Transfer global particle IDs:
1003           call getClusterParticleData(iglobal_id, gpIDs_all, gpIDs)
1004           if(dbg_level >= 2) then
1005              call dbg_print_clusters(7)
1006              write(202,"(5x,'check 7')")
1007           endif
1008     
1009     
1010     ! Transfer particle ghost array:
1011           call getClusterGhostData(gpGhost_all, gpGhost)
1012           if(dbg_level >= 2) then
1013              call dbg_print_clusters(7)
1014              write(202,"(5x,'check 7')")
1015              call dbg_print_clusters(9)
1016              write(202,"(5x,'check 9')")
1017              call dbg_print_clusters(10)
1018              write(202,"(5x,'check 10')")
1019              call dbg_print_clusters(11)
1020              write(202,"(5x,'check 11')")
1021           elseif(dbg_level >= 1) then
1022              call dbg_print_clusters(11)
1023           endif
1024     
1025     
1026           if(myPE == clusterPE) then
1027     ! If there are two or more clusters, check to see if there is any of
1028     ! the clusters intersect. (Merge clusters that span multiple processes.)
1029              if( cCnt_sum > 1) then
1030                 allocate(mergeMap(cCnt_sum, cCnt_sum)); mergeMap = 0
1031                 do lc1 = 1, cCnt_sum
1032                    mergeMap(lc1,lc1) = 1
1033                 enddo
1034     ! Loop over clusters, comparing the global IDs of the particles that
1035     ! they contain. Any clusters with overlapping particles are later
1036     ! merged as a single cluster.
1037                 lp_lc1: do lc1 = 1, cCnt_sum-1
1038                 lp_lc2: do lc2 = 1+gp_dsp(lc1), gp_dsp(lc1)+pCnt_all(lc1)
1039                 lp_lc3: do lc3 = lc1+1, cCnt_sum
1040                 lp_lc4: do lc4 = 1+gp_dsp(lc3), gp_dsp(lc3)+pCnt_all(lc3)
1041                    if(gpIDs_all(lc2) == gpIDs_all(lc4)) then
1042                       mergeMap(lc1,lc3) = 1
1043                       mergeMap(lc3,lc1) = 1
1044                       if(dbg_level >= 1) &
1045                          write(*,"(3x,'Merge: ',I5,' and ',I5)") lc1, lc3
1046                    endif
1047                 enddo lp_lc4
1048                 enddo lp_lc3
1049                 enddo lp_lc2
1050                 enddo lp_lc1
1051                 if(dbg_level >= 1) then
1052                    write(*,*)''
1053                    call dbg_print_clusters(8, lmsg='before')
1054                 endif
1055     
1056     ! Merge any clusters that share common particles. Common particles are
1057     ! identified by their global particle IDs. Note that if two processes
1058     ! contain a common particle, the particle should be 'real' on one
1059     ! process and a ghost on the other.
1060                 merging = .true.
1061                 do while(merging)
1062                    merging = .false.
1063                    lp_lc11: do lc1=1,cCnt_sum-1
1064                       if(sum(mergeMap(lc1,:)) == 0) cycle lp_lc11
1065                       lp_lc21: do lc2=lc1+1,cCnt_sum
1066                          if(mergeMap(lc1,lc2) /= 1) cycle lp_lc21
1067                          lc3=lc2
1068                          lp_lc41: do lc4 = 1, cCnt_sum
1069                             if(lc3 == lc4) cycle lp_lc41
1070                             if(mergeMap(lc3,lc4) == 1) then
1071                                mergeMap(lc1,lc4) = 1
1072                                mergeMap(lc3,lc4) = 0
1073                                mergeMap(lc3,lc3) = 0
1074                                merging = .true.
1075                             endif
1076                          enddo lp_lc41
1077                       enddo lp_lc21
1078                    enddo lp_lc11
1079                 enddo
1080                 if(dbg_level >= 1) &
1081                    call dbg_print_clusters(8, lmsg='after')
1082     
1083     ! Count the actual number of clusters in the domain.
1084                 clusterCount_all = 0
1085                 do lc1=1,cCnt_sum
1086                    if(sum(mergeMap(lc1,:)) > 0) then
1087                       clusterCount_all = clusterCount_all + 1
1088                    endif
1089                 enddo
1090                 if(dbg_level >= 1) then
1091                    write(*,"(3x,'Number of clusters reported by ', &
1092                       &'all processes: ',I6)") cCnt_sum
1093                    write(*,"(3x,'Actual number of clusters: ',I6)") &
1094                       clusterCount_all
1095                    write(*,*)''
1096                 endif
1097     
1098     
1099     ! Allocate the clusters pointer array. This array is populated with
1100     ! data used to map particles to send/recv data. Ghost particles are not
1101     ! included in the map so that cluster analysis only contains data from
1102     ! real particles.
1103                 allocate(clusters(clusterCount_all))
1104                 do lc1=1,clusterCount_all
1105                    nullify(clusters(lc1)%particle)
1106                    clusters(lc1)%size = 0
1107                 enddo
1108                 lc3 = 0
1109                 lp_lc12: do lc1=1, cCnt_sum
1110                    if(sum(mergeMap(lc1,:)) == 0) cycle lp_lc12
1111                    lc3 = lc3 + 1
1112                    this => clusters(lc3); this%size = 0
1113                    lp_lc22: do lc2=1, cCnt_sum
1114                       if(mergeMap(lc1,lc2) == 0) cycle lp_lc22
1115                       lp_lc42: do lc4=1+gp_dsp(lc2),gp_dsp(lc2)+pCnt_all(lc2)
1116                          if(gpGhost_all(lc4)) cycle lp_lc42
1117                          call addParticle(this, lc4, gpIDs_all(lc4))
1118                       enddo lp_lc42
1119                    enddo lp_lc22
1120                 enddo lp_lc12
1121                 deallocate(mergeMap)
1122                 if(dbg_level >= 1) &
1123                    call dbg_print_clusters(12, dbg=dbg_level)
1124     
1125     ! If there is only one cluster, there is no need to check for
1126     ! intersecting clusters.
1127              elseif( cCnt_sum == 1) then
1128                 if(dbg_level >= 1) &
1129                    write(*,"(3x,'Just one cluster: No merge.',/)")
1130                 clusterCount_all = 1
1131                 allocate(clusters(clusterCount_all))
1132                 nullify(clusters(clusterCount_all)%particle)
1133                 clusters(clusterCount_all)%size = 0
1134                 this => clusters(clusterCount_all); this%size = 0
1135                 do lc1=1, pCnt_all(1)
1136                    if(gpGhost_all(lc1)) write(*,"(3x, &
1137                       &'Error: Ghost particle detected: ',I7)") lc1
1138                    call addParticle(this, lc1, gpIDs_all(lc1))
1139                 enddo
1140                 if(dbg_level >= 1) &
1141                    call dbg_print_clusters(12, dbg=dbg_level)
1142              else
1143                 if(dbg_level >= 1) write(*,"(3x,'No clusters to merge.')")
1144                 clusterCount_all = 0
1145                 allocate(clusters(1))
1146                 nullify(clusters(1)%particle)
1147                 clusters(1)%size = 0
1148              endif
1149           else
1150              clusterCount_all = 0
1151              allocate( clusters(1) )
1152              nullify(clusters(1)%particle)
1153              clusters(1)%size = 0
1154           endif
1155           call des_mpi_barrier()
1156     
1157     ! clean up before returning.
1158           nullify(cluster)
1159           nullify(particle)
1160     
1161           if(allocated(pCnt))      deallocate(pCnt)
1162           if(allocated(gpIDs))     deallocate(gpIDs)
1163           if(allocated(gpGhost))     deallocate(gpGhost)
1164     
1165           if(allocated(pCnt_all))  deallocate(pCnt_all)
1166           if(allocated(gpIDs_all)) deallocate(gpIDs_all)
1167           if(allocated(gpGhost_all)) deallocate(gpGhost_all)
1168     
1169           if(allocated(gp_dsp))    deallocate(gp_dsp)
1170     
1171           return
1172     
1173           contains
1174     
1175     !......................................................................!
1176     !  Module name: dbg_PRINT_CLUSTERS                                     !
1177     !                                                                      !
1178     !  Purpose: A collection of various write statements for debugging.    !
1179     !                                                                      !
1180     !  Author: J.Musser                                   Date:  Dec-12    !
1181     !......................................................................!
1182           SUBROUTINE dbg_print_clusters(lmsgID, lmsg, dbg)
1183           use run
1184     
1185           implicit none
1186     
1187     ! Index for printing a specific message.
1188           INTEGER, INTENT(IN) :: lmsgID
1189     
1190     ! A message to write with data.
1191           CHARACTER(len=*), intent(in), optional :: lmsg
1192     
1193           INTEGER, INTENT(IN), optional :: dbg
1194     
1195     ! Generic loop counters
1196           INTEGER proc, lc1, lc2, lc3, lc4
1197     ! Generic write buffers.
1198           CHARACTER(len=120) wbuff, wbuff2
1199     ! String for generating filenames.
1200           CHARACTER(LEN=255) :: filename
1201     
1202     ! Generic cluster pointer.
1203           Type(cType), pointer :: cThis
1204     ! Generic cluster pointer.
1205           Type(pType), pointer :: pThis
1206     
1207           SELECT CASE(lmsgID)
1208     
1209     !``````````````````````````````````````````````````````````````````````!
1210     ! Root Process Only -                                                  !
1211     !                                                                      !
1212     ! Purpose:: Message to write out the number of clusters detected on    !
1213     ! each process.                                                        !
1214     !                                                                      !
1215     !``````````````````````````````````````````````````````````````````````!
1216           CASE (1); if(myPE /= clusterPE) return
1217     
1218              do proc=0,numPEs-1
1219                 write(*,"(3x,'Process ',I2,' reporting ',I4,' clusters.')") &
1220                    proc, cCnt_all(proc)
1221              enddo
1222              write(*,"(3x,'Total reported clusters: ',I6)") sum(cCnt_all)
1223              write(*,*)''
1224     
1225     
1226     !``````````````````````````````````````````````````````````````````````!
1227     ! All Processes -                                                      !
1228     !                                                                      !
1229     ! Purpose:: Each process creates a file and records the particles      !
1230     ! contained in clusters. Only the data from the current pass is kept   !
1231     ! (file status='replace').                                             !
1232     !                                                                      !
1233     !``````````````````````````````````````````````````````````````````````!
1234           CASE (2) ! All processes
1235     
1236              filename = ''
1237              write(filename,"('dbg_pCnt_',I2.2,'.txt')") myPE
1238              open(convert='big_endian',unit=201, file=filename, status='replace')
1239              write(201,"(//3x,'Time:',F18.6)")Time
1240              write(201,"(3x,'Number of Clusters: ',I4)") cCnt(myPE)
1241     
1242              if(cCnt(myPE) > 0) then
1243                 nullify(cluster)
1244                 do lc2 = 1, clusterCount
1245                    write(201,"(3x,'/Cluster ',I5,' has ',I8,' members.')")
1246                    write(201,"(5x,'Membership:')")
1247                    CALL getNextCluster(cluster)
1248                    NULLIFY(particle)
1249                    do lc1 = 1, cluster%ParticleCount
1250                       CALL GetNextParticle(cluster, particle)
1251                       write(201,"(5x,'Global ID: ',I8)") &
1252                          iglobal_id(particle%ID)
1253                    enddo
1254                 enddo
1255              endif
1256              close(201)
1257     
1258     
1259     !``````````````````````````````````````````````````````````````````````!
1260     ! Root Process Only -                                                  !
1261     !                                                                      !
1262     ! Purpose:: Root process constructs a file that lists all the particle !
1263     ! from all processes that are contained in clusters. Only data from    !
1264     ! the current pass is kept (file status='replace').                    !
1265     !                                                                      !
1266     !``````````````````````````````````````````````````````````````````````!
1267           CASE (3); if(myPE /= clusterPE) return
1268     
1269              filename = 'dbg_pCnt_dsp.txt'
1270              open(convert='big_endian',unit=201, file=filename, status='replace')
1271              write(201,"(/3x,'Time:',F18.6)")Time
1272              write(201,"(3x,'cCnt_sum: ',I4)") cCnt_sum
1273              do lc1 =0, numPEs-1
1274                 write(201,"(5x,'pCnt_dsp(',I4,'): ',I6)")lc1, pCnt_dsp(lc1)
1275              enddo
1276              close(201)
1277     
1278     
1279     !``````````````````````````````````````````````````````````````````````!
1280     ! Root Process Only -                                                  !
1281     !                                                                      !
1282     ! Purpose:: Root process reports the number of clusters on each        !
1283     ! process and the number of particles comprising each cluster.         !
1284     !                                                                      !
1285     !``````````````````````````````````````````````````````````````````````!
1286           CASE (4); if(myPE /= clusterPE) return
1287     
1288              if(cCnt_sum > 0)then
1289                 write(*,"(3x,'Total number of clusters: ',I6)") cCnt_sum
1290                 do proc = 0, numPEs-1
1291                    if(cCnt_all(proc) > 0) then
1292                       write(*,"(5x,'Process ',I2,': cluster count: ',I4)") &
1293                          proc, cCnt_all(proc)
1294                       do lc1 = 1, cCnt_all(proc)
1295                          write(*,"(7x,'Particles in cluster ',I4,': ',I6)")&
1296                             lc1, pCnt_all(lc1 + pCnt_dsp(proc))
1297                       enddo
1298                    else
1299                       write(*,"(3x,'Process ',I2,' reports no clusters.')")&
1300                          proc
1301                    endif
1302                 enddo
1303                 write(*,*)''
1304              endif
1305     
1306     
1307     !``````````````````````````````````````````````````````````````````````!
1308     ! Root Process Only -                                                  !
1309     !                                                                      !
1310     ! Purpose:: Root process reports the information on the send/recv      !
1311     ! process for cluster_gather routines.                                 !
1312     !                                                                      !
1313     !``````````````````````````````````````````````````````````````````````!
1314           CASE (5); if(myPE /= clusterPE) return
1315     
1316              if(recv_sum > 0)then
1317                 write(*,"(3x,'send_cnt: ',I8)") send_cnt
1318                 write(*,"(3x,'recv_sum: ',I8)") recv_sum
1319                 write(*,"(3x,'Generic send/recv setup:')")
1320                 do proc = 0, numPEs-1
1321                    if(recv_cnt(proc) > 0) then
1322                       write(*,"(5x,'Process ',I2,': count: ',I6,&
1323                          &' disp:',I6)") &
1324                          proc, recv_cnt(proc), recv_dsp(proc)
1325                    else
1326                       write(*,"(3x,'Process ',I2,' is empty.')") proc
1327                    endif
1328                 enddo
1329              else
1330                 write(*,"(3x,'No data to send/recv.')")
1331              endif
1332              write(*,*)''
1333     
1334     
1335     !``````````````````````````````````````````````````````````````````````!
1336     ! Root Process Only -                                                  !
1337     !                                                                      !
1338     ! Purpose:: Root process reports the displacements for particle data   !
1339     ! stored in the receive buffer for cluster_gather.                     !
1340     !                                                                      !
1341     !``````````````````````````````````````````````````````````````````````!
1342           CASE (6); if(myPE /= clusterPE) return
1343     
1344              if(cCnt_sum > 0)then
1345                 write(*,"(3x,'Global cluster-to-particle offset:')")
1346                 do lc1 = 1, cCnt_sum
1347                    write(*,"(7x,' Cluster ',I4,': ',I6)") &
1348                       lc1, gp_dsp(lc1)
1349                 enddo
1350                 write(*,*)''
1351              endif
1352     
1353     
1354     !``````````````````````````````````````````````````````````````````````!
1355     ! Root Process Only -                                                  !
1356     !                                                                      !
1357     ! Purpose:: Root process constructs a file containing the global IDs   !
1358     ! of all particles contained in clusters.                              !
1359     !                                                                      !
1360     !``````````````````````````````````````````````````````````````````````!
1361           CASE (7); if(myPE /= clusterPE) return
1362     
1363              open(convert='big_endian',unit=201, file='dbg_gpIDs.txt', status='replace')
1364              do lc1=1,size(gpIDs_all)
1365                 write(201,"(5x,' Global particle ID: ',I8)") gpIDs_all(lc1)
1366              enddo
1367     
1368              close(201)
1369     
1370     
1371     !``````````````````````````````````````````````````````````````````````!
1372     ! Root Process Only -                                                  !
1373     !                                                                      !
1374     ! Purpose:: Root process reports cluster merging data. Output is based !
1375     ! on the array data set in call to this routine.                       !
1376     !                                                                      !
1377     !``````````````````````````````````````````````````````````````````````!
1378           CASE (8); if(myPE /= clusterPE) return
1379     
1380              if(.not.present(lmsg)) then
1381                 write(*,"(/3x,'Invalid use of case debug report.')")
1382                 return
1383              endif
1384     
1385              filename = ''
1386              write(filename,"('dbg_mergeMap_',A,I2.2,'.txt')") trim(lmsg)
1387              open(convert='big_endian',unit=201, file=filename, status='replace')
1388     
1389              write(201,"(/3x,'Time:',F18.6)")Time
1390     
1391              lc1 = size(mergeMap(:,1))         ! rows
1392              lc2 = min(size(mergeMap(1,:)),25) ! columns
1393     
1394              wbuff = ''
1395              write(wbuff,"(7x,'1')")
1396              do lc3=2,lc2
1397                 write(wbuff,"(A,3x,I1)") trim(wbuff), lc3
1398              enddo
1399     
1400              wbuff2 = ''
1401              write(wbuff2,"(5x,'|')")
1402              do lc3=1,lc2
1403                 write(wbuff2,"(2A)") trim(wbuff2),'---|'
1404              enddo
1405     
1406              write(201,*)trim(wbuff)
1407              write(201,*)trim(wbuff2)
1408     
1409              do lc3=1,lc1
1410                 wbuff = ''
1411                 write(wbuff,"(3X,I1,' |')") lc3
1412                 do lc4 = 1, lc2
1413                    write(wbuff,"(A,1x,I1,' |')") &
1414                       trim(wbuff), mergeMap(lc3,lc4)
1415                 enddo
1416                 write(201,*)trim(wbuff)
1417                 write(201,*)trim(wbuff2)
1418              enddo
1419              close(201)
1420     
1421     
1422     !``````````````````````````````````````````````````````````````````````!
1423     ! All Processes -                                                      !
1424     !                                                                      !
1425     ! Purpose:: Each process constructs a file and records ghost particle  !
1426     ! information.                                                         !
1427     !                                                                      !
1428     !``````````````````````````````````````````````````````````````````````!
1429           CASE (9) ! All processes
1430     
1431              filename = ''
1432              write(filename,"('dbg_isghost_',I2.2,'.txt')") myPE
1433              open(convert='big_endian',unit=201, file=filename, status='replace')
1434              write(201,"(/3x,'Time:',F18.6)")Time
1435              do lc1 =1, send_cnt
1436                 write(201,"(5x,'IS_GHOST(',I8,'): ',L2)")lc1, gpGhost(lc1)
1437              enddo
1438              close(201)
1439     
1440     
1441     !``````````````````````````````````````````````````````````````````````!
1442     ! Root Process Only -                                                  !
1443     !                                                                      !
1444     ! Purpose:: Root process constructs a file and records ghost particle  !
1445     ! information for all particles contained in clusters.                 !
1446     !                                                                      !
1447     !``````````````````````````````````````````````````````````````````````!
1448           CASE (10); if(myPE /= clusterPE) return
1449     
1450              open(convert='big_endian',unit=201, file='dbg_isghost_all.txt', status='replace')
1451              write(201,"(/3x,'Time:',F18.6)")Time
1452              do lc1 =1, recv_sum
1453                 write(201,"(5x,'IS_GHOST(',I8,'): ',L2)")lc1, gpGhost_all(lc1)
1454              enddo
1455              close(201)
1456     
1457     
1458     !``````````````````````````````````````````````````````````````````````!
1459     ! Root Process Only -                                                  !
1460     !                                                                      !
1461     ! Purpose:: Root process reports particles in clusters that are        !
1462     ! identified as ghost particles on a process.                          !
1463     !                                                                      !
1464     !``````````````````````````````````````````````````````````````````````!
1465           CASE (11); if(myPE /= clusterPE) return
1466     
1467              if(recv_sum <= 0) return
1468     
1469              lc2 = 0
1470              do proc = 0, numPEs-1
1471              do lc1 = 1, cCnt_all(proc)
1472              lc2 = lc2 + 1
1473              do lc3 = 1, pCnt_all(lc2)
1474                 if(gpGhost_all(gp_dsp(lc2) + lc3)) then
1475                    write(*,"(3x,'Particle ',I8,' in cluster ',I6, &
1476                       &' is a ghost on process ',I2,'.')")        &
1477                       gpIDs_all(gp_dsp(lc2) + lc3), lc2, proc
1478                 endif
1479              enddo
1480              enddo
1481              enddo
1482              write(*,*)''
1483     
1484     
1485     !``````````````````````````````````````````````````````````````````````!
1486     ! Root Process Only -                                                  !
1487     !                                                                      !
1488     ! Purpose:: Root process reports cluster statistics.                   !
1489     !                                                                      !
1490     !``````````````````````````````````````````````````````````````````````!
1491           CASE (12); if(myPE /= clusterPE) return
1492     
1493              if(.not.present(dbg)) return
1494     
1495              if(dbg >= 1) then
1496                 do lc1=1, clusterCount_all
1497                    write(*,"(3x,'Cluster ',I6,' reports ',I6,&
1498                       &' particles.')") lc1, clusters(lc1)%size
1499                 enddo
1500              endif
1501     
1502              if(dbg >= 2) then
1503                 do lc1=1, clusterCount_all
1504                    filename = ''
1505                    write(filename,"('dbg_cluster_',I3.3,'.txt')")lc1
1506                    open(convert='big_endian',unit=201, file=filename, status='replace')
1507                    write(201,"(3x,'Time:',F10.6)") Time
1508     
1509                    cThis => clusters(lc1)
1510     
1511                    write(201,"(3x,'Cluster ',I6,' reports ',I6,&
1512                       &' particles.')") lc1, cThis%size
1513                    write(201,"(3x,' Particles in this cluster include:')")
1514     
1515                    nullify(pThis)
1516                    if(associated(cThis%particle)) pThis => cThis%particle
1517                    do while(associated(pThis))
1518                       write(201,"(3x,I8)") pThis%id
1519                       if(associated(pThis%next)) then
1520                          pThis => pThis%next
1521                       else
1522                          nullify(pThis)
1523                       endif
1524                    enddo
1525                    close(201)
1526                    nullify(cThis)
1527                 enddo
1528              endif
1529     
1530     
1531     !``````````````````````````````````````````````````````````````````````!
1532     ! Root Process Only -                                                  !
1533     !                                                                      !
1534     ! Purpose:: Default message for invalid message id.                    !
1535     !                                                                      !
1536     !``````````````````````````````````````````````````````````````````````!
1537           CASE DEFAULT; if(myPE /= clusterPE) return
1538              WRITE(*,"(3x,'No message exists for msgID: ',I4)")lmsgID
1539              write(*,*)''
1540           END SELECT
1541     
1542           return
1543           END SUBROUTINE dbg_print_clusters
1544     
1545           END SUBROUTINE INIT_PRINT_CLUSTERS
1546     
1547     
1548     !......................................................................!
1549     !  Module name: FINL_PRINT_CLUSTERS                                    !
1550     !                                                                      !
1551     !  Purpose: Deallocate storage arrays used for collecting cluster data.!
1552     !                                                                      !
1553     !  Author: J.Musser                                   Date:  Dec-12    !
1554     !......................................................................!
1555           SUBROUTINE finl_print_clusters(dbg_level)
1556     
1557           implicit none
1558     
1559     ! Flags used for debugging.
1560     !   0: No messages
1561     !   1: Screen only messages
1562     !   2: Detailed Files are written by each process
1563           INTEGER, intent(in) :: dbg_level
1564     
1565           integer lc1
1566           type(cType), pointer :: cThis
1567           type(pType), pointer :: pThis
1568     
1569           if(dbg_level >= 2)  write(202,"(3x,'check B')")
1570           call des_mpi_barrier()
1571     
1572           if(allocated(recv_cnt)) deallocate(recv_cnt)
1573           if(allocated(recv_dsp)) deallocate(recv_dsp)
1574     
1575           if(mype /= clusterPE ) then
1576              if(allocated(clusters)) deallocate(clusters)
1577           else
1578              if(allocated(clusters)) then
1579                 do lc1=1,size(clusters)
1580                    cThis => clusters(lc1)
1581                    do while(associated(cThis%particle))
1582                       if(associated(cThis%particle%next)) then
1583                          pThis => cThis%particle
1584                          cThis%particle => pThis%next
1585                          nullify(pThis%next)
1586                          deallocate(pThis)
1587                          nullify(pThis)
1588                       else
1589                          deallocate(cThis%particle)
1590                          nullify(cThis%particle)
1591                       endif
1592                       cThis%size = cThis%size - 1
1593                    enddo
1594                    if(cThis%size /= 0) then
1595                       write(*,"(3x,'Error deallocating clusters: ', I6)")&
1596                       cThis%size
1597                    endif
1598                 enddo
1599                 deallocate(clusters)
1600              endif
1601           endif
1602     
1603     ! Write debug message.
1604           if(dbg_level >= 2) close(202)
1605           if(myPE == clusterPE .and. dbg_level >=1) &
1606              write(*,"(/1x,'End data dump ',52('-'),'<'//)")
1607     
1608     ! Clean up local pointers.
1609           nullify(pThis)
1610           nullify(cThis)
1611           END SUBROUTINE finl_print_clusters
1612     
1613     
1614     !......................................................................!
1615     !  Module name: addParticle                                            !
1616     !                                                                      !
1617     !  Purpose: Given a specific cluster (this) a new particle object is   !
1618     !  is created and added to the cluster link-list.                      !
1619     !                                                                      !
1620     !  Author: J.Musser                                   Date:  Dec-12    !
1621     !......................................................................!
1622           SUBROUTINE addParticle(this, lMap, lID)
1623     
1624     ! Cluster object getting a new partilced added.
1625           Type(cType), pointer, intent(inout) :: this
1626     
1627     ! The particle map is used for locating a particle in send/recv data.
1628           integer, intent(in) :: lMap
1629     ! The global ID of the particle. Used mainly for debugging.
1630           integer, intent(in) :: lID
1631     
1632     ! The new particle being added to the cluster.
1633           Type(pType), pointer :: new
1634     
1635     ! New particle constructor: Allocate the memory and set its entries.
1636           allocate(new)
1637           nullify(new%next)
1638           new%map = lMap
1639           new%id = lID
1640     
1641     ! Update the cluster size (the number of particles in the cluster)
1642           this%size = this%size + 1
1643     
1644     ! Add the new particle to the cluster.
1645           if(this%size == 1) then
1646     ! If this is the first particle, then the cluster should have a null
1647     ! entry for the particle link-list. (Sanity check)
1648              if(associated(this%particle)) then
1649                 write(*,"(3x,'Fatal Error (000) adding particle.')")
1650                 return
1651              endif
1652           else
1653     ! If this isn't the first particle, then the cluster should have a
1654     ! valid particle in the link list.
1655              if(.not.associated(this%particle)) then
1656                 write(*,"(3x,'Fatal Error (001) adding particle.')")
1657                 return
1658              endif
1659              new%next => this%particle
1660           endif
1661     
1662           this%particle => new
1663           nullify(new)
1664     
1665           return
1666           END SUBROUTINE addParticle
1667     
1668     
1669     !......................................................................!
1670     !  Module name: getClusterParticleData_1i                              !
1671     !                                                                      !
1672     !  Purpose: Collect data from local processes and send to clusterPE.   !
1673     !  The data from cluster-bound particles is placed in a temporary      !
1674     !  array for send/recv.                                                !
1675     !                                                                      !
1676     !  Incoming data is a 1D Array of Integers                             !
1677     !                                                                      !
1678     !  Author: J.Musser                                   Date:  Dec-12    !
1679     !......................................................................!
1680           SUBROUTINE getClusterParticleData_1i(lData, lrbuff, lOut)
1681     
1682           implicit none
1683     
1684     ! Data on local process being sent to clusterPE.
1685           integer, intent(in) :: lData(:)
1686     ! Data from local processes received by clusterPE.
1687           integer, allocatable, intent(inout) :: lrbuff(:)
1688     ! Data on local process in cluster-array format.
1689           integer, allocatable, intent(inout), optional :: lOut(:)
1690     
1691     ! Temporary storage array for local data. This data is in array format.
1692           integer, allocatable :: lsbuff(:)
1693     ! Clusters detected by local search
1694           type(cluster_type),  pointer :: cluster
1695     ! Particles classified as being in a cluster on the local process.
1696           type(particle_type), pointer :: particle
1697     
1698     ! Generic loop counters
1699           integer lc1, lc2, lc3
1700     ! Error flag.
1701           integer ierr
1702     
1703     ! The receive buffer should be unallocated. Deallocate it otherwise.
1704     ! This array is only important on clusterPE.
1705           if(allocated(lrbuff)) then
1706              if(myPE == clusterPE) write(*,"(3x,&
1707                 &'Error in getClusterParticleData_1i: ',&
1708                 &'Deallocating receive buffer.')")
1709              deallocate(lrbuff)
1710           endif
1711     
1712     ! Allocate the receive buffer (and returned data set)
1713           allocate(lrbuff(recv_sum))
1714     ! Allocate the temporary send buffer.
1715           allocate(lsbuff(send_cnt))
1716     ! Allocate the local return data set
1717           if(present(lOut)) allocate(lOut(send_cnt))
1718     
1719     ! Populate the send buffer. The send buffer only contains data about
1720     ! particles that belong to clusters.
1721           if(send_cnt > 0) then
1722              NULLIFY(cluster)
1723              lc3 = 0
1724              do lc2 = 1, clusterCount
1725                 CALL getNextCluster(cluster)
1726                 NULLIFY(particle)
1727                 do lc1 = 1, cluster%ParticleCount
1728                    lc3 = lc3 + 1
1729                    CALL GetNextParticle(cluster, particle)
1730                    lsbuff(lc3) = lData(particle%ID)
1731                 enddo
1732              enddo
1733           endif
1734     
1735     ! Invoke MPI routines.
1736           call des_mpi_gatherv(lsbuff, send_cnt,             &
1737                                lrbuff, recv_cnt, recv_dsp,   &
1738                                clusterPE, ierr)
1739     
1740     ! Store the local return data set.
1741           if(present(lOut)) lOut = lsbuff
1742     
1743     ! Clean up.
1744           deallocate(lsbuff)
1745           nullify(cluster)
1746           nullify(particle)
1747     
1748           return
1749           END SUBROUTINE getClusterParticleData_1i
1750     
1751     
1752     !......................................................................!
1753     !  Module name: getClusterParticleData_2i                              !
1754     !                                                                      !
1755     !  Purpose: Collect data from local processes and send to clusterPE.   !
1756     !  The data from cluster-bound particles is placed in a temporary      !
1757     !  array for send/recv.                                                !
1758     !                                                                      !
1759     !  Incoming data is a 2D Array of Integers                             !
1760     !                                                                      !
1761     !  Author: J.Musser                                   Date:  Dec-12    !
1762     !......................................................................!
1763           SUBROUTINE getClusterParticleData_2i(lData, lrbuff, lOut)
1764     
1765           implicit none
1766     
1767     ! Data on local process being sent to clusterPE.
1768           integer, intent(in) :: lData(:,:)
1769     ! Data from local processes received by clusterPE.
1770           integer, allocatable, intent(inout) :: lrbuff(:,:)
1771     ! Data on local process in cluster-array format.
1772           integer, allocatable, intent(inout), optional :: lOut(:,:)
1773     ! Temporary storage array for local data. This data is in array format.
1774           integer, allocatable :: lsbuff(:)
1775     
1776     ! Clusters detected by local search
1777           type(cluster_type),  pointer :: cluster
1778     ! Particles classified as being in a cluster on the local process.
1779           type(particle_type), pointer :: particle
1780     
1781     ! Generic loop counters
1782           integer lc1, lc2, lc3, lc4
1783     ! Upper/Lower array bounds of lData position 2
1784           integer lbnd, ubnd
1785     ! Error flag.
1786           integer ierr
1787     
1788     ! Get the upper and lower array bounds of the incoming data.
1789           lbnd = lbound(lData,1)
1790           ubnd = ubound(lData,1)
1791     
1792     ! The receive buffer should be unallocated. Deallocate it otherwise.
1793     ! This array is only important on clusterPE.
1794           if(allocated(lrbuff)) then
1795              if(myPE == 0) write(*,"(3x,&
1796                 &'Error in getClusterParticleData_2i: ', &
1797                 &'Deallocating receive buffer.')")
1798              deallocate(lrbuff)
1799           endif
1800     
1801     ! Allocate the receive buffer (and returned data set)
1802           allocate(lrbuff(lbnd:ubnd,recv_sum))
1803     ! Allocate the temporary send buffer.
1804           allocate(lsbuff(send_cnt))
1805     ! Allocate the local return data set
1806           if(present(lOut)) allocate(lOut(lbnd:ubnd,send_cnt))
1807     
1808     ! Populate the send buffer. The send buffer only contains data about
1809     ! particles that belong to clusters.
1810           do lc4=lbnd,ubnd
1811              if(send_cnt > 0) then
1812                 lsbuff = -987654321
1813                 NULLIFY(cluster)
1814                 lc3 = 0
1815                 do lc2 = 1, clusterCount
1816                    CALL getNextCluster(cluster)
1817                    NULLIFY(particle)
1818                    do lc1 = 1, cluster%ParticleCount
1819                       lc3 = lc3 + 1
1820                       CALL GetNextParticle(cluster, particle)
1821                       lsbuff(lc3) = lData(lc4,particle%ID)
1822                    enddo
1823                 enddo
1824              endif
1825     
1826     ! Invoke MPI routines.
1827              call des_mpi_gatherv(lsbuff,        send_cnt,             &
1828                                   lrbuff(lc4, :), recv_cnt, recv_dsp,   &
1829                                   clusterPE, ierr)
1830     ! Store the local return data set.
1831              if(present(lOut)) lOut(lc4,:) = lsbuff
1832     
1833           enddo
1834     
1835     ! Clean up.
1836           deallocate(lsbuff)
1837           nullify(cluster)
1838           nullify(particle)
1839     
1840           return
1841           END SUBROUTINE getClusterParticleData_2i
1842     
1843     
1844     !......................................................................!
1845     !  Module name: getClusterParticleData_1d                              !
1846     !                                                                      !
1847     !  Purpose: Collect data from local processes and send to clusterPE.   !
1848     !  The data from cluster-bound particles is placed in a temporary      !
1849     !  array for send/recv.                                                !
1850     !                                                                      !
1851     !  Incoming data is a 1D Array of double precision values.             !
1852     !                                                                      !
1853     !  Author: J.Musser                                   Date:  Dec-12    !
1854     !......................................................................!
1855           SUBROUTINE getClusterParticleData_1d(lData, lrbuff, lOut)
1856     
1857           implicit none
1858     
1859     ! Data on local process being sent to clusterPE.
1860           double precision, intent(in) :: lData(:)
1861     ! Data from local processes received by clusterPE.
1862           double precision, allocatable, intent(inout) :: lrbuff(:)
1863     ! Data on local process in cluster-array format.
1864           double precision, allocatable, intent(inout), optional :: lOut(:)
1865     
1866     ! Temporary storage array for local data. This data is in array format.
1867           double precision, allocatable :: lsbuff(:)
1868     
1869     ! Clusters detected by local search
1870           type(cluster_type),  pointer :: cluster
1871     ! Particles classified as being in a cluster on the local process.
1872           type(particle_type), pointer :: particle
1873     
1874     ! Generic loop counters
1875           INTEGER lc1, lc2, lc3
1876     ! Error flag.
1877           integer ierr
1878     
1879     ! The receive buffer should be unallocated. Deallocate it otherwise.
1880     ! This array is only important on clusterPE.
1881           if(allocated(lrbuff)) then
1882              if(myPE == 0) write(*,"(3x,&
1883                 &'Error in getClusterParticleData_1d: ', &
1884                 &'Deallocating receive buffer.')")
1885              deallocate(lrbuff)
1886           endif
1887     
1888     ! Allocate the receive buffer (and returned data set)
1889           allocate(lrbuff(recv_sum))
1890     ! Allocate the temporary send buffer.
1891           allocate(lsbuff(send_cnt))
1892     ! Allocate the local return data set
1893           if(present(lOut)) allocate(lOut(send_cnt))
1894     
1895     ! Populate the send buffer. The send buffer only contains data about
1896     ! particles that belong to clusters.
1897           if(send_cnt > 0) then
1898              NULLIFY(cluster)
1899              lc3 = 0
1900              do lc2 = 1, clusterCount
1901                 CALL getNextCluster(cluster)
1902                 NULLIFY(particle)
1903                 do lc1 = 1, cluster%ParticleCount
1904                    lc3 = lc3 + 1
1905                    CALL GetNextParticle(cluster, particle)
1906                    lsbuff(lc3) = lData(particle%ID)
1907                 enddo
1908              enddo
1909           endif
1910     
1911     ! Invoke MPI routines.
1912           call des_mpi_gatherv(lsbuff, send_cnt,             &
1913                                lrbuff, recv_cnt, recv_dsp,   &
1914                                clusterPE, ierr)
1915     
1916     ! Store the local return data set.
1917           if(present(lOut)) lOut = lsbuff
1918     
1919     ! Clean up.
1920           deallocate(lsbuff)
1921           nullify(cluster)
1922           nullify(particle)
1923     
1924           return
1925           END SUBROUTINE getClusterParticleData_1d
1926     
1927     
1928     !......................................................................!
1929     !  Module name: getClusterParticleData_2d                              !
1930     !                                                                      !
1931     !  Purpose: Collect data from local processes and send to clusterPE.   !
1932     !  The data from cluster-bound particles is placed in a temporary      !
1933     !  array for send/recv.                                                !
1934     !                                                                      !
1935     !  Incoming data is a 2D Array of double precision values.             !
1936     !                                                                      !
1937     !  Author: J.Musser                                   Date:  Dec-12    !
1938     !......................................................................!
1939           SUBROUTINE getClusterParticleData_2d(lData, lrbuff, lOut)
1940     
1941           implicit none
1942     
1943     ! Data on local process being sent to clusterPE.
1944           double precision, intent(in) :: lData(:,:)
1945     ! Data from local processes received by clusterPE.
1946           double precision, allocatable, intent(inout) :: lrbuff(:,:)
1947     ! Data on local process in cluster-array format.
1948           double precision, allocatable, intent(inout), optional :: lOut(:,:)
1949     
1950     ! Temporary storage array for local data. This data is in array format.
1951           double precision, allocatable :: lsbuff(:) ! send buffer
1952     
1953     ! Clusters detected by local search
1954           type(cluster_type),  pointer :: cluster
1955     ! Particles classified as being in a cluster on the local process.
1956           type(particle_type), pointer :: particle
1957     
1958     ! Generic loop counters
1959           integer lc1, lc2, lc3, lc4
1960     ! Upper/Lower array bounds of lData position 2
1961           integer lbnd, ubnd
1962     ! Error flag.
1963           integer ierr
1964     
1965     ! Get the upper and lower array bounds of the incoming data.
1966           lbnd = lbound(lData,1)
1967           ubnd = ubound(lData,1)
1968     
1969     ! The receive buffer should be unallocated. Deallocate it otherwise.
1970     ! This array is only important on clusterPE.
1971           if(allocated(lrbuff)) then
1972              if(myPE == 0) write(*,"(3x,&
1973                 &'Error in getClusterParticleData_2d: ', &
1974                 &'Deallocating receive buffer.')")
1975              deallocate(lrbuff)
1976           endif
1977     
1978     ! Allocate the receive buffer (and returned data set)
1979           allocate(lrbuff(lbnd:ubnd,recv_sum))
1980     ! Allocate the temporary send buffer.
1981           allocate(lsbuff(send_cnt))
1982     ! Allocate the local return data set
1983           if(present(lOut)) allocate(lOut(lbnd:ubnd,send_cnt))
1984     
1985     ! Populate the send buffer. The send buffer only contains data about
1986     ! particles that belong to clusters.
1987           do lc4=lbnd,ubnd
1988              if(send_cnt > 0) then
1989                 lsbuff = -9.87654321
1990                 NULLIFY(cluster)
1991                 lc3 = 0
1992                 do lc2 = 1, clusterCount
1993                    CALL getNextCluster(cluster)
1994                    NULLIFY(particle)
1995                    do lc1 = 1, cluster%ParticleCount
1996                       lc3 = lc3 + 1
1997                       CALL GetNextParticle(cluster, particle)
1998                       lsbuff(lc3) = lData(lc4,particle%ID)
1999                    enddo
2000                 enddo
2001              endif
2002     
2003     ! Invoke MPI routines.
2004              call des_mpi_gatherv(lsbuff,        send_cnt,             &
2005                                   lrbuff(lc4,:), recv_cnt, recv_dsp,   &
2006                                   clusterPE, ierr)
2007     
2008     ! Store the local return data set.
2009              if(present(lOut)) lOut(lc4,:) = lsbuff
2010     
2011           enddo
2012     
2013     ! Clean up.
2014           deallocate(lsbuff)
2015           nullify(cluster)
2016           nullify(particle)
2017     
2018           return
2019           END SUBROUTINE getClusterParticleData_2d
2020     
2021     
2022     !......................................................................!
2023     !  Module name: getClusterParticleData_1l                              !
2024     !                                                                      !
2025     !  Purpose: Collect data from local processes and send to clusterPE.   !
2026     !  The data from cluster-bound particles is placed in a temporary      !
2027     !  array for send/recv.                                                !
2028     !                                                                      !
2029     !  Incoming data is a 1D Array of logicals.                            !
2030     !                                                                      !
2031     !  Note: Logicals are converted to integers for send/recv, then        !
2032     !  converted back to logical values.                                   !
2033     !                                                                      !
2034     !  Author: J.Musser                                   Date:  Dec-12    !
2035     !......................................................................!
2036           SUBROUTINE getClusterGhostData(lrbuff, lOut)
2037     
2038           implicit none
2039     
2040     ! Data from local processes received by clusterPE.
2041           logical, allocatable, intent(inout) :: lrbuff(:)
2042     ! Data on local process in cluster-array format.
2043           logical, allocatable, intent(inout), optional :: lOut(:)
2044     
2045     ! Data on local process in cluster-array format (converted to integer)
2046           integer, allocatable :: lsbuff_i(:)
2047     ! Data from local processes received by clusterPE (converted to integer)
2048           integer, allocatable :: lrbuff_i(:)
2049     
2050     ! Clusters detected by local search
2051           type(cluster_type),  pointer :: cluster
2052     ! Particles classified as being in a cluster on the local process.
2053           type(particle_type), pointer :: particle
2054     
2055     ! Generic loop counters
2056           INTEGER lc1, lc2, lc3
2057     ! Error flag.
2058           integer ierr
2059     
2060     ! The receive buffer should be unallocated. Deallocate it otherwise.
2061     ! This array is only important on clusterPE.
2062           if(allocated(lrbuff)) then
2063              if(myPE == 0) write(*,"(3x,&
2064                 &'Error in getClusterParticleData_1l: ',&
2065                 &'Deallocating receive buffer.')")
2066              deallocate(lrbuff)
2067           endif
2068     
2069     ! Allocate the receive buffer (and returned data set)
2070           allocate(lrbuff(recv_sum)); lrbuff = .false.
2071     ! Allocate the local return data set.
2072           if(present(lOut)) allocate(lOut(send_cnt))
2073     
2074     ! Allocate the temporary receiver buffer.
2075           allocate(lrbuff_i(recv_sum)); lrbuff_i = 0
2076     ! Allocate the temporary send buffer.
2077           allocate(lsbuff_i(send_cnt)); lsbuff_i = 0
2078     
2079     ! Populate the send buffer. The send buffer only contains data about
2080     ! particles that belong to clusters.
2081           if(send_cnt > 0) then
2082              NULLIFY(cluster)
2083              lc3 = 0
2084              do lc2 = 1, clusterCount
2085                 CALL getNextCluster(cluster)
2086                 NULLIFY(particle)
2087                 do lc1 = 1, cluster%ParticleCount
2088                    lc3 = lc3 + 1
2089                    CALL GetNextParticle(cluster, particle)
2090     ! Convert logical to integer.
2091                    if(IS_GHOST(particle%ID)) lsbuff_i(lc3) = 1
2092                 enddo
2093              enddo
2094           endif
2095     ! Invoke MPI routines.
2096           call des_mpi_gatherv(lsbuff_i, send_cnt,             &
2097                                lrbuff_i, recv_cnt, recv_dsp,   &
2098                                clusterPE, ierr)
2099     
2100     ! Convert integer back to logical.
2101           do lc1=1,recv_sum
2102              if(lrbuff_i(lc1) == 1) lrbuff(lc1) = .true.
2103           enddo
2104     
2105     ! Convert and store the local return data set.
2106           if(present(lOut)) then
2107              lOut = .false.
2108              do lc1 = 1,  send_cnt
2109                 if(lsbuff_i(lc1) == 1) lOut(lc1) = .true.
2110              enddo
2111           endif
2112     
2113     ! Clean up.
2114           deallocate(lrbuff_i)
2115           deallocate(lsbuff_i)
2116           nullify(cluster)
2117           nullify(particle)
2118     
2119           return
2120         END SUBROUTINE getClusterGhostData
2121     
2122     !......................................................................!
2123     !  Module name: getClusterFieldData_1d                                 !
2124     !                                                                      !
2125     !  Purpose: Collect field (continuum) data from local processes and    !
2126     !  send to clusterPE.  The data from the field associated with         !
2127     !  cluster-bound particles is accessed by the particke map and not     !
2128     !  the global IJK value.                                               !
2129     !                                                                      !
2130     !  Incoming data is a 1D Array of double precision values.             !
2131     !                                                                      !
2132     !                                                                      !
2133     !  Author: J.Musser                                   Date:  Dec-12    !
2134     !......................................................................!
2135           SUBROUTINE getClusterFieldData_1d(lData, lrbuff, lOut)
2136     
2137           implicit none
2138     
2139     ! Field data on local process being sent to clusterPE.
2140           double precision, intent(in) :: lData(:)
2141     ! Field data from local processes received by clusterPE.
2142           double precision, allocatable, intent(inout) :: lrbuff(:)
2143     ! Field ata on local process in cluster-array format.
2144           double precision, allocatable, intent(inout), optional :: lOut(:)
2145     
2146           double precision, allocatable :: lsbuff(:) ! send buffer
2147     
2148     ! Clusters detected by local search
2149           type(cluster_type),  pointer :: cluster
2150     ! Particles classified as being in a cluster on the local process.
2151           type(particle_type), pointer :: particle
2152     
2153     ! Generic loop counters
2154           integer lc1, lc2, lc3
2155     ! IJK index of fluid cell containing particle
2156           integer ijk
2157     ! Error flag.
2158           integer ierr
2159     
2160     ! The receive buffer should be unallocated. Deallocate it otherwise.
2161     ! This array is only important on clusterPE.
2162           if(allocated(lrbuff)) then
2163              if(myPE == 0) write(*,"(3x,&
2164                 &'Error in getClusterFieldData_1d: ', &
2165                 &'Deallocating receive buffer.')")
2166              deallocate(lrbuff)
2167           endif
2168     
2169     ! Allocate the receive buffer (and returned data set)
2170           allocate(lrbuff(recv_sum))
2171     ! Allocate the temporary send buffer.
2172           allocate(lsbuff(send_cnt))
2173     ! Allocate the local return data set.
2174           if(present(lOut)) allocate(lOut(send_cnt))
2175     
2176     ! Populate the send buffer. The send buffer only contains data about
2177     ! particles that belong to clusters.
2178           if(send_cnt > 0) then
2179              NULLIFY(cluster)
2180              lc3 = 0
2181              do lc2 = 1, clusterCount
2182                 CALL getNextCluster(cluster)
2183                 NULLIFY(particle)
2184                 do lc1 = 1, cluster%ParticleCount
2185                    lc3 = lc3 + 1
2186                    CALL GetNextParticle(cluster, particle)
2187     ! Use the particke fluid cell IJK to access the field data.
2188                    ijk = PIJK(particle%ID,4)
2189     ! Ghost particles have an IJK of 0. Since data from ghost particles is
2190     ! not needed, set the value as undefined.
2191                    if(ijk == 0) then
2192                       lsbuff(lc3) = undefined
2193                    else
2194                       lsbuff(lc3) = lData(ijk)
2195                    endif
2196                 enddo
2197              enddo
2198           endif
2199     
2200     ! Invoke MPI routines.
2201           call des_mpi_gatherv(lsbuff, send_cnt,             &
2202                                lrbuff, recv_cnt, recv_dsp,   &
2203                                clusterPE, ierr)
2204     
2205     ! Store the local return data set.
2206           if(present(lOut)) lOut = lsbuff
2207     
2208     ! Clean up.
2209           deallocate(lsbuff)
2210           nullify(cluster)
2211           nullify(particle)
2212     
2213           return
2214           END SUBROUTINE getClusterFieldData_1d
2215     
2216     
2217     !......................................................................!
2218     !  Module name: getClusterFieldData_3d                                 !
2219     !                                                                      !
2220     !  Purpose: Collect field (continuum) data from local processes and    !
2221     !  send to clusterPE.  The data from the field associated with         !
2222     !  cluster-bound particles is placed in a temporary and is no longer   !
2223     !  IJK bound.                                                          !
2224     !                                                                      !
2225     !  Incoming data are three 1D Array of double precision values.        !
2226     !                                                                      !
2227     !                                                                      !
2228     !  Author: J.Musser                                   Date:  Dec-12    !
2229     !......................................................................!
2230           SUBROUTINE getClusterFieldData_3d(lData_1, lData_2, lData_3, &
2231              lrbuff, lOut)
2232     
2233           implicit none
2234     
2235     ! Data on local process being sent to clusterPE.
2236           double precision, intent(in) :: lData_1(:) ! x-axis
2237           double precision, intent(in) :: lData_2(:) ! y-axis
2238           double precision, intent(in) :: lData_3(:) ! z-axis
2239     
2240     ! Data from local processes received by clusterPE.
2241           double precision, allocatable, intent(inout) :: lrbuff(:,:)
2242     ! Data on local process in cluster-array format.
2243           double precision, allocatable, intent(inout), optional :: lOut(:,:)
2244     ! Temporary storage array for local data. This data is in array format.
2245           double precision, allocatable :: lsbuff(:) ! send buffer
2246     
2247     ! Clusters detected by local search
2248           type(cluster_type),  pointer :: cluster
2249     ! Particles classified as being in a cluster on the local process.
2250           type(particle_type), pointer :: particle
2251     
2252     ! Generic loop counters
2253           integer lc1, lc2, lc3
2254     ! IJK index of fluid cell containing particle
2255           integer ijk
2256     ! Error flag.
2257           integer ierr
2258     
2259     ! The receive buffer should be unallocated. Deallocate it otherwise.
2260     ! This array is only important on clusterPE.
2261           if(allocated(lrbuff)) then
2262              if(myPE == 0) write(*,"(3x,&
2263                 &'Error in getClusterFieldData_3d: ', &
2264                 &'Deallocating receive buffer.')")
2265              deallocate(lrbuff)
2266           endif
2267     
2268     ! Allocate the receive buffer (and returned data set)
2269           allocate(lrbuff(recv_sum,1:3))
2270     ! Allocate the temporary send buffer.
2271           allocate(lsbuff(send_cnt))
2272     ! Allocate the local return data set
2273           if(present(lOut)) allocate(lOut(send_cnt,1:3))
2274     
2275     ! Populate the send buffer. The send buffer only contains data about
2276     ! particles that belong to clusters.
2277           if(send_cnt > 0) then
2278              lsbuff = -9.87654321
2279              NULLIFY(cluster)
2280              lc3 = 0
2281              do lc2 = 1, clusterCount
2282                 CALL getNextCluster(cluster)
2283                 NULLIFY(particle)
2284                 do lc1 = 1, cluster%ParticleCount
2285                    lc3 = lc3 + 1
2286                    CALL GetNextParticle(cluster, particle)
2287     ! Use the particke fluid cell IJK to access the field data.
2288                    ijk = PIJK(particle%ID,4)
2289     ! Ghost particles have an IJK of 0. Since data from ghost particles is
2290     ! not needed, set the value as undefined.
2291                    if(ijk == 0) then
2292                       lsbuff(lc3) = undefined
2293                    else
2294                       lsbuff(lc3) = lData_1(ijk)
2295                    endif
2296                 enddo
2297              enddo
2298           endif
2299     ! Invoke MPI routines.
2300           call des_mpi_gatherv(lsbuff,      send_cnt,             &
2301                                lrbuff(:,1), recv_cnt, recv_dsp,   &
2302                                clusterPE, ierr)
2303     ! Store the local return data set.
2304           if(present(lOut)) lOut(:,1) = lsbuff
2305     
2306           if(send_cnt > 0) then
2307              lsbuff = -9.87654321
2308              NULLIFY(cluster)
2309              lc3 = 0
2310              do lc2 = 1, clusterCount
2311                 CALL getNextCluster(cluster)
2312                 NULLIFY(particle)
2313                 do lc1 = 1, cluster%ParticleCount
2314                    lc3 = lc3 + 1
2315                    CALL GetNextParticle(cluster, particle)
2316     ! Use the particke fluid cell IJK to access the field data.
2317                    ijk = PIJK(particle%ID,4)
2318     ! Ghost particles have an IJK of 0. Since data from ghost particles is
2319     ! not needed, set the value as undefined.
2320                    if(ijk == 0) then
2321                       lsbuff(lc3) = undefined
2322                    else
2323                       lsbuff(lc3) = lData_2(ijk)
2324                    endif
2325                 enddo
2326              enddo
2327           endif
2328     ! Invoke MPI routines.
2329           call des_mpi_gatherv(lsbuff,      send_cnt,             &
2330                                lrbuff(:,2), recv_cnt, recv_dsp,   &
2331                                clusterPE, ierr)
2332     ! Store the local return data set.
2333           if(present(lOut)) lOut(:,2) = lsbuff
2334     
2335           if(DO_K) then
2336              if(send_cnt > 0) then
2337                 lsbuff = -9.87654321
2338                 NULLIFY(cluster)
2339                 lc3 = 0
2340                 do lc2 = 1, clusterCount
2341                    CALL getNextCluster(cluster)
2342                    NULLIFY(particle)
2343                    do lc1 = 1, cluster%ParticleCount
2344                       lc3 = lc3 + 1
2345                       CALL GetNextParticle(cluster, particle)
2346     ! Use the particke fluid cell IJK to access the field data.
2347                        ijk = PIJK(particle%ID,4)
2348     ! Ghost particles have an IJK of 0. Since data from ghost particles is
2349     ! not needed, set the value as undefined.
2350                       if(ijk == 0) then
2351                          lsbuff(lc3) = undefined
2352                       else
2353                          lsbuff(lc3) = lData_3(ijk)
2354                       endif
2355                    enddo
2356                 enddo
2357              endif
2358     ! Invoke MPI routines.
2359              call des_mpi_gatherv(lsbuff,      send_cnt,             &
2360                                   lrbuff(:,3), recv_cnt, recv_dsp,   &
2361                                   clusterPE, ierr)
2362     ! Store the local return data set.
2363              if(present(lOut)) lOut(:,3) = lsbuff
2364     
2365           endif
2366     
2367     ! Clean up.
2368           deallocate(lsbuff)
2369           nullify(cluster)
2370           nullify(particle)
2371     
2372           return
2373           END SUBROUTINE getClusterFieldData_3d
2374     
2375     
2376     !......................................................................!
2377     !  Module name: sendClusterData_1d                                     !
2378     !                                                                      !
2379     !  Purpose: Take cluster data from clusterPE and send it to each       !
2380     !  process. (Storage of output data for visualization)                 !
2381     !                                                                      !
2382     !  Note: The local storage array is initialized to zero by default.    !
2383     !  A different value may be used by specifying initValue.              !
2384     !                                                                      !
2385     !                                                                      !
2386     !  Author: J.Musser                                   Date:  Dec-12    !
2387     !......................................................................!
2388           SUBROUTINE sendClusterData_1d(lsbuff, lData)
2389     
2390           implicit none
2391     
2392     ! Data on clusterPE being distributed to the other processes.
2393           double precision, intent(in) :: lsbuff(:)
2394     ! Storage array on local process.
2395           double precision, intent(inout) :: lData(:)
2396     
2397     ! Data from clusterPE received by local process.
2398           double precision, allocatable :: lrbuff(:)
2399     
2400     ! Clusters detected by local search
2401           type(cluster_type),  pointer :: cluster
2402     ! Particles classified as being in a cluster on the local process.
2403           type(particle_type), pointer :: particle
2404     
2405     ! Generic loop counters
2406           integer lc1, lc2, lc3
2407     ! Error flag.
2408           integer ierr
2409     
2410     ! Allocate the temporary receive buffer on local process.
2411           allocate(lrbuff(send_cnt))
2412     
2413     ! Invoke MPI routines.
2414           call des_mpi_scatterv(lsbuff, recv_cnt, recv_dsp, lrbuff, &
2415              send_cnt, clusterPE, ierr)
2416     
2417     ! Map the received data back to individual particles using the local
2418     ! cluster data.
2419           if(send_cnt > 0) then
2420              NULLIFY(cluster)
2421              lc3 = 0
2422              do lc2 = 1, clusterCount
2423                 CALL getNextCluster(cluster)
2424                 NULLIFY(particle)
2425                 do lc1 = 1, cluster%ParticleCount
2426                    lc3 = lc3 + 1
2427                    CALL GetNextParticle(cluster, particle)
2428                    lData(particle%ID) = lrbuff(lc3)
2429                 enddo
2430              enddo
2431           endif
2432     
2433     ! Clean up.
2434           deallocate(lrbuff)
2435           nullify(cluster)
2436           nullify(particle)
2437     
2438           return
2439           END SUBROUTINE sendClusterData_1d
2440     
2441     
2442     
2443           END MODULE DES_CLUSTER
2444