File: /nfs/home/0/users/jenkins/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     
20           IMPLICIT NONE
21     
22     !-----------------------------------------------
23     ! Local variables
24     !-----------------------------------------------
25           INTEGER I, J, K
26           INTEGER L, LL
27           DOUBLE PRECISION DISTVEC(3), DIST, R_LM
28     ! Temporary variables to adjust particle position in event of periodic
29     ! boundaries
30           DOUBLE PRECISION XPOS(2), YPOS(2), ZPOS(2), TMPPOS(3)
31     ! Max loop limit in each coordinate direction
32           INTEGER II, JJ, KK
33     ! Index to track accounted for particles
34           INTEGER PC
35     ! Index to track unaccounted for particles
36           INTEGER PNPC
37     !-----------------------------------------------
38     
39           PC=1
40           DO L=1, MAX_PIP
41     
42              IF(PC .GE. PIP ) EXIT
43              IF(.NOT.PEA(L,1)) CYCLE
44     
45              PNPC = PIP - PC
46              DO LL = L+1, MAX_PIP
47     
48                 IF(PNPC .LE. 0) EXIT
49                 IF(.NOT.PEA(LL,1)) CYCLE
50     
51                 R_LM = DES_RADIUS(L) + DES_RADIUS(LL)
52                 R_LM = FACTOR_RLM*R_LM
53     
54     ! the following section adjusts the neighbor check routine when
55     ! any boundary is periodic
56     ! ------------------------------
57                 IF (DES_PERIODIC_WALLS) THEN
58                    XPOS(:) = DES_POS_NEW(1,LL)
59                    YPOS(:) = DES_POS_NEW(2,LL)
60                    II = 1
61                    JJ = 1
62                    KK = 1
63     
64                    IF(DES_PERIODIC_WALLS_X) THEN
65                       IF (DES_POS_NEW(1,L) + R_LM > XLENGTH) THEN
66                          II = 2
67                          XPOS(II) = DES_POS_NEW(1,LL) + XLENGTH
68                       ELSEIF (DES_POS_NEW(1,L) - R_LM < ZERO) THEN
69                          II = 2
70                          XPOS(II) = DES_POS_NEW(1,LL) - XLENGTH
71                       ENDIF
72                    ENDIF
73                    IF(DES_PERIODIC_WALLS_Y) THEN
74                       IF (DES_POS_NEW(2,L) + R_LM > YLENGTH) THEN
75                          JJ = 2
76                          YPOS(JJ) = DES_POS_NEW(2,LL) + YLENGTH
77                       ELSEIF (DES_POS_NEW(2,L) - R_LM < YLENGTH) THEN
78                          JJ = 2
79                          YPOS(JJ) = DES_POS_NEW(2,LL) - YLENGTH
80                       ENDIF
81                    ENDIF
82                    IF(DO_K) THEN
83                       ZPOS(:) = DES_POS_NEW(3,LL)
84                       IF(DES_PERIODIC_WALLS_Z) THEN
85                          IF (DES_POS_NEW(3,L) + R_LM > ZLENGTH) THEN
86                             KK = 2
87                             ZPOS(KK) = DES_POS_NEW(3,LL) + ZLENGTH
88                          ELSEIF (DES_POS_NEW(3,L) - R_LM < ZERO) THEN
89                             KK = 2
90                             ZPOS(KK) = DES_POS_NEW(3,LL) - ZLENGTH
91                          ENDIF
92                       ENDIF
93                    ENDIF
94     
95     ! if particle L is within R_LM of a periodic boundary then check
96     ! particles LL current position and its position shifted to the
97     ! opposite boundary for neighbor contact
98                    OUTER: DO I = 1,II
99                       DO J = 1,JJ
100                          TMPPOS(1) = XPOS(I)
101                          TMPPOS(2) = YPOS(J)
102                          IF (DO_K) THEN
103                             DO K = 1,KK
104                                TMPPOS(3) = ZPOS(K)
105                                DISTVEC(:) = TMPPOS(:) - DES_POS_NEW(:,L)
106                                DIST = dot_product(DISTVEC,DISTVEC)
107                                IF (DIST.LE.R_LM) EXIT OUTER
108                             ENDDO
109                          ELSE
110                             DISTVEC(:) = TMPPOS(:) - DES_POS_NEW(:,L)
111                             DIST = dot_product(DISTVEC,DISTVEC)
112                             IF (DIST.LE.R_LM) EXIT OUTER
113                          ENDIF
114                       ENDDO
115                    ENDDO OUTER
116     
117                 ELSE   ! if .not.des_periodic_walls
118                    DISTVEC(:) = DES_POS_NEW(:,LL) - DES_POS_NEW(:,L)
119                    DIST = dot_product(DISTVEC,DISTVEC)
120                 ENDIF    ! endif des_periodic_walls
121     ! ------------------------------
122     
123                 IF (DIST < R_LM**2) THEN
124                    call add_pair(L, LL)
125                 ENDIF
126                 PNPC = PNPC - 1
127              ENDDO   ! end loop over LL
128     
129              PC = PC + 1
130           ENDDO   ! end loop over L
131     
132           RETURN
133           END SUBROUTINE NSQUARE
134     
135     
136     
137