File: RELATIVE:/../../../mfix.git/model/des/nsquare.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: NSQUARE                                                C
4     !>  Purpose: DES - N-Square neighbor search
5     !                                                                      C
6     !                                                                      C
7     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
8     !  Reviewer: Sreekanth Pannala                        Date: 09-Nov-06  C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11     
12           SUBROUTINE NSQUARE
13     
14           use discretelement
15           use des_bc
16           use des_thermo
17           use des_allocate
18           use geometry, only: DO_K, xlength, ylength, zlength
19           use functions
20     
21           IMPLICIT NONE
22     
23     !-----------------------------------------------
24     ! Local variables
25     !-----------------------------------------------
26           INTEGER I, J, K
27           INTEGER L, LL, CC
28           DOUBLE PRECISION DISTVEC(3), DIST, R_LM
29     ! Temporary variables to adjust particle position in event of periodic
30     ! boundaries
31           DOUBLE PRECISION XPOS(2), YPOS(2), ZPOS(2), TMPPOS(3)
32     ! Max loop limit in each coordinate direction
33           INTEGER II, JJ, KK
34     ! Index to track accounted for particles
35           INTEGER PC
36     ! Index to track unaccounted for particles
37           INTEGER PNPC
38     !-----------------------------------------------
39     
40           PC=1
41           DO L=1, MAX_PIP
42     
43              NEIGHBOR_INDEX(L) = 1
44              if (L .gt. 1) NEIGHBOR_INDEX(L) = NEIGHBOR_INDEX(L-1)
45     
46              IF(PC .GE. PIP ) EXIT
47              IF(IS_NONEXISTENT(L)) CYCLE
48     
49              PNPC = PIP - PC
50              DO LL = L+1, MAX_PIP
51     
52                 IF(PNPC .LE. 0) EXIT
53                 IF(IS_NONEXISTENT(LL)) CYCLE
54     
55                 R_LM = DES_RADIUS(L) + DES_RADIUS(LL)
56                 R_LM = FACTOR_RLM*R_LM
57     
58     ! the following section adjusts the neighbor check routine when
59     ! any boundary is periodic
60     ! ------------------------------
61                 IF (DES_PERIODIC_WALLS) THEN
62                    XPOS(:) = DES_POS_NEW(1,LL)
63                    YPOS(:) = DES_POS_NEW(2,LL)
64                    II = 1
65                    JJ = 1
66                    KK = 1
67     
68                    IF(DES_PERIODIC_WALLS_X) THEN
69                       IF (DES_POS_NEW(1,L) + R_LM > XLENGTH) THEN
70                          II = 2
71                          XPOS(II) = DES_POS_NEW(1,LL) + XLENGTH
72                       ELSEIF (DES_POS_NEW(1,L) - R_LM < ZERO) THEN
73                          II = 2
74                          XPOS(II) = DES_POS_NEW(1,LL) - XLENGTH
75                       ENDIF
76                    ENDIF
77                    IF(DES_PERIODIC_WALLS_Y) THEN
78                       IF (DES_POS_NEW(2,L) + R_LM > YLENGTH) THEN
79                          JJ = 2
80                          YPOS(JJ) = DES_POS_NEW(2,LL) + YLENGTH
81                       ELSEIF (DES_POS_NEW(2,L) - R_LM < YLENGTH) THEN
82                          JJ = 2
83                          YPOS(JJ) = DES_POS_NEW(2,LL) - YLENGTH
84                       ENDIF
85                    ENDIF
86                    IF(DO_K) THEN
87                       ZPOS(:) = DES_POS_NEW(3,LL)
88                       IF(DES_PERIODIC_WALLS_Z) THEN
89                          IF (DES_POS_NEW(3,L) + R_LM > ZLENGTH) THEN
90                             KK = 2
91                             ZPOS(KK) = DES_POS_NEW(3,LL) + ZLENGTH
92                          ELSEIF (DES_POS_NEW(3,L) - R_LM < ZERO) THEN
93                             KK = 2
94                             ZPOS(KK) = DES_POS_NEW(3,LL) - ZLENGTH
95                          ENDIF
96                       ENDIF
97                    ENDIF
98     
99     ! if particle L is within R_LM of a periodic boundary then check
100     ! particles LL current position and its position shifted to the
101     ! opposite boundary for neighbor contact
102                    OUTER: DO I = 1,II
103                       DO J = 1,JJ
104                          TMPPOS(1) = XPOS(I)
105                          TMPPOS(2) = YPOS(J)
106                          IF (DO_K) THEN
107                             DO K = 1,KK
108                                TMPPOS(3) = ZPOS(K)
109                                DISTVEC(:) = TMPPOS(:) - DES_POS_NEW(:,L)
110                                DIST = dot_product(DISTVEC,DISTVEC)
111                                IF (DIST.LE.R_LM) EXIT OUTER
112                             ENDDO
113                          ELSE
114                             DISTVEC(:) = TMPPOS(:) - DES_POS_NEW(:,L)
115                             DIST = dot_product(DISTVEC,DISTVEC)
116                             IF (DIST.LE.R_LM) EXIT OUTER
117                          ENDIF
118                       ENDDO
119                    ENDDO OUTER
120     
121                 ELSE   ! if .not.des_periodic_walls
122                    DISTVEC(:) = DES_POS_NEW(:,LL) - DES_POS_NEW(:,L)
123                    DIST = dot_product(DISTVEC,DISTVEC)
124                 ENDIF    ! endif des_periodic_walls
125     ! ------------------------------
126     
127                 IF (DIST < R_LM**2) THEN
128                    cc = add_pair(L, LL)
129                 ENDIF
130                 PNPC = PNPC - 1
131              ENDDO   ! end loop over LL
132     
133              PC = PC + 1
134           ENDDO   ! end loop over L
135     
136           RETURN
137           END SUBROUTINE NSQUARE
138