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