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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: CALC_INTERP_WEIGHTS                                     !
3     !                                                                      !
4     !  Purpose: Calculate weights used for interpolation.                  !
5     !                                                                      !
6     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
7           SUBROUTINE CALC_INTERP_WEIGHTS
8     
9           use particle_filter, only: DES_INTERP_SCHEME_ENUM
10           use particle_filter, only: DES_INTERP_DPVM
11           use particle_filter, only: DES_INTERP_GAUSS
12           use particle_filter, only: DES_INTERP_LHAT
13     
14           SELECT CASE(DES_INTERP_SCHEME_ENUM)
15           CASE(DES_INTERP_DPVM);  CALL CALC_INTERP_WEIGHTS1
16           CASE(DES_INTERP_GAUSS); CALL CALC_INTERP_WEIGHTS1
17           CASE(DES_INTERP_LHAT);  CALL CALC_INTERP_WEIGHTS2
18           END SELECT
19     
20           RETURN
21           END SUBROUTINE CALC_INTERP_WEIGHTS
22     
23     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
24     !  Subroutine: CALC_INTERP_WEIGHTS1                                    !
25     !                                                                      !
26     !  Purpose: Calculate weights used for interpolation.                  !
27     !                                                                      !
28     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
29           SUBROUTINE CALC_INTERP_WEIGHTS1
30     
31           use discretelement, only: MAX_PIP
32           use discretelement, only: PIJK
33           use discretelement, only: DES_POS_NEW
34           use discretelement, only: XE, YN, ZT
35           use particle_filter, only: FILTER_CELL
36           use particle_filter, only: FILTER_WEIGHT
37           use geometry, only: DO_K
38     
39           use param1, only: ZERO, ONE
40     ! compar is needed by functions.inc
41           use compar
42     
43           IMPLICIT NONE
44     
45           INTEGER :: L, IDX
46           INTEGER :: IJK, IJKt
47     
48           INTEGER :: I, J, K
49           INTEGER :: IC, JC, KC
50           INTEGER :: Km, Kp
51     
52           DOUBLE PRECISION :: WEIGHT
53           DOUBLE PRECISION :: WEIGHT_I(-1:1)
54           DOUBLE PRECISION :: WEIGHT_J(-1:1)
55           DOUBLE PRECISION :: WEIGHT_K(-1:1)
56     
57     !$omp parallel default(none) private(L, WEIGHT_I, WEIGHT_J, WEIGHT_K,  &
58     !$omp    KM, KP, IDX, IJK, I, J, K, WEIGHT, IJKT)                      &
59     !$omp shared(MAX_PIP, PIJK, DES_POS_NEW, XE, YN, ZT, DO_K,        &
60     !$omp    FILTER_CELL, FILTER_WEIGHT)
61     !$omp do
62           DO L = 1, MAX_PIP
63     
64              IF(.NOT.IS_NORMAL(L)) CYCLE
65     
66              I = PIJK(L,1)
67              J = PIJK(L,2)
68              K = PIJK(L,3)
69     
70     ! Tentative weights for I indices to the West and East.
71              WEIGHT_I(-1) = CALC_FILTER_WEIGHTS(DES_POS_NEW(L,1), XE(I-1))
72              WEIGHT_I( 1) = CALC_FILTER_WEIGHTS(XE(I), DES_POS_NEW(L,1))
73              WEIGHT_I( 0) = ONE - WEIGHT_I(-1) - WEIGHT_I(1)
74     
75     ! Tentative weights for J indices to the South and North.
76              WEIGHT_J(-1) = CALC_FILTER_WEIGHTS(DES_POS_NEW(L,2), YN(J-1))
77              WEIGHT_J( 1) = CALC_FILTER_WEIGHTS(YN(J), DES_POS_NEW(L,2))
78              WEIGHT_J( 0) = ONE - WEIGHT_J(-1) - WEIGHT_J(1)
79     
80     ! Tentative weights for K indices to the Top and Bottom.
81              IF(DO_K) THEN
82                 Km=-1;  Kp=1
83                 WEIGHT_K(-1) = CALC_FILTER_WEIGHTS(DES_POS_NEW(L,3),ZT(K-1))
84                 WEIGHT_K( 1) = CALC_FILTER_WEIGHTS(ZT(K), DES_POS_NEW(L,3))
85                 WEIGHT_K( 0) = ONE - WEIGHT_K(-1) - WEIGHT_K(1)
86              ELSE
87                 Km= 0; Kp=0
88                 WEIGHT_K( 0) = ONE
89              ENDIF
90     
91     ! Set the default fluid cell index and loop counter.
92              IJK = PIJK(L,4)
93              IDX=0
94     
95     ! Calculate weights for ghost particles. Only store weights that the
96     ! current process owns.
97              DO KC=Km,Kp
98              DO IC=-1,+1
99              DO JC=-1,+1
100                 IDX=IDX+1
101                 WEIGHT = WEIGHT_I(IC)*WEIGHT_J(JC)*WEIGHT_K(KC)
102                 IF(IS_ON_MYPE_PLUS2LAYERS(I+IC,J+JC,K+KC)) THEN
103                    IJKt = FUNIJK(I+IC,J+JC,K+KC)
104                    IF(FLUID_AT(IJKt) .OR. CYCLIC_AT(IJKt)) THEN
105                       FILTER_CELL(IDX,L) = IJKt
106                       FILTER_WEIGHT(IDX,L) = WEIGHT
107                    ELSE
108                       FILTER_CELL(IDX,L) = IJK
109                       FILTER_WEIGHT(IDX,L) = WEIGHT
110                    ENDIF
111                 ELSE
112                    FILTER_CELL(IDX,L) = IJK
113                    FILTER_WEIGHT(IDX,L) = ZERO
114                 ENDIF
115              ENDDO
116              ENDDO
117              ENDDO
118     
119           ENDDO
120     !$omp end do
121     !$omp end parallel
122     
123           CONTAINS
124     
125     !``````````````````````````````````````````````````````````````````````!
126     !                                                                      !
127     !``````````````````````````````````````````````````````````````````````!
128           DOUBLE PRECISION FUNCTION CALC_FILTER_WEIGHTS(POS1, POS2)
129     
130           use particle_filter, only: FILTER_WIDTH_INTERP
131           use particle_filter, only: OoFILTER_VOL
132           use particle_filter, only: FILTER_WIDTH_INTERPx3
133     
134           DOUBLE PRECISION, INTENT(IN) :: POS1, POS2
135           DOUBLE PRECISION :: OVERLAP, HEIGHT
136     
137           OVERLAP = POS1 - POS2
138           IF(OVERLAP < FILTER_WIDTH_INTERP) THEN
139              HEIGHT = FILTER_WIDTH_INTERP - OVERLAP
140              CALC_FILTER_WEIGHTS = HEIGHT**2 * &
141                 (FILTER_WIDTH_INTERPx3 - HEIGHT)*OoFILTER_VOL
142           ELSE
143              CALC_FILTER_WEIGHTS = ZERO
144           ENDIF
145     
146           END FUNCTION CALC_FILTER_WEIGHTS
147     
148           INCLUDE 'functions.inc'
149     
150           END SUBROUTINE CALC_INTERP_WEIGHTS1
151     
152     
153     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
154     !  Subroutine: CALC_INTERP_WEIGHTS2                                    !
155     !                                                                      !
156     !  Purpose: Calculate weights used for interpolation.                  !
157     !                                                                      !
158     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
159           SUBROUTINE CALC_INTERP_WEIGHTS2
160     
161           use discretelement, only: MAX_PIP
162           use discretelement, only: PIJK
163           use discretelement, only: DES_POS_NEW
164           use discretelement, only: XE, YN, ZT
165           use particle_filter, only: FILTER_CELL
166           use particle_filter, only: FILTER_WEIGHT
167           use geometry, only: DX, DY, DZ
168           use geometry, only: oDX_E, oDY_N, oDZ_T, DO_K
169     
170           use param1, only: ZERO, HALF, ONE
171     
172           use compar, only: FUNIJK_MAP_C
173     
174     
175           IMPLICIT NONE
176     
177           INTEGER :: L, IDX
178           INTEGER :: IJK, IJKt
179     
180           INTEGER :: I, J, K
181           INTEGER :: IC, JC, KC
182           INTEGER :: Km, Kp
183     
184           DOUBLE PRECISION :: XC, YC, ZC
185           DOUBLE PRECISION :: WEIGHT
186           DOUBLE PRECISION :: WEIGHT_I(-1:1)
187           DOUBLE PRECISION :: WEIGHT_J(-1:1)
188           DOUBLE PRECISION :: WEIGHT_K(-1:1)
189     
190     
191           DO L = 1, MAX_PIP
192     
193              IF(.NOT.IS_NORMAL(L)) CYCLE
194     
195              I = PIJK(L,1)
196              J = PIJK(L,2)
197              K = PIJK(L,3)
198     
199     ! Tentative weights for I indices to the West and East.
200              XC = XE(I-1) + HALF*DX(I)
201              IF(DES_POS_NEW(L,1) < XE(I-1)+HALF*DX(I)) THEN
202                 WEIGHT_I(-1) = (XC - DES_POS_NEW(L,1))*oDX_E(I-1)
203                 WEIGHT_I( 0) = ONE - WEIGHT_I(-1)
204                 WEIGHT_I( 1) = ZERO
205              ELSE
206                 WEIGHT_I( 1) = (DES_POS_NEW(L,1) - XC)*oDX_E(I)
207                 WEIGHT_I( 0) = ONE - WEIGHT_I(1)
208                 WEIGHT_I(-1) = ZERO
209              ENDIF
210     
211     ! Tentative weights for J indices to the South and North.
212              YC = YN(J-1) + HALF*DY(J)
213              IF(DES_POS_NEW(L,2) < YN(J-1)+HALF*DY(J)) THEN
214                 WEIGHT_J(-1) = (YC - DES_POS_NEW(L,2))*oDY_N(J-1)
215                 WEIGHT_J( 0) = ONE - WEIGHT_J(-1)
216                 WEIGHT_J( 1) = ZERO
217              ELSE
218                 WEIGHT_J( 1) = (DES_POS_NEW(L,2) - YC)*oDY_N(J)
219                 WEIGHT_J( 0) = ONE - WEIGHT_J(1)
220                 WEIGHT_J(-1) = ZERO
221              ENDIF
222     
223     ! Tentative weights for K indices to the Top and Bottom.
224              IF(DO_K) THEN
225                 Km=-1;  Kp=1
226                 ZC = ZT(K-1) + HALF*DZ(K)
227                 IF(DES_POS_NEW(L,3) < ZT(K-1)+HALF*DZ(K)) THEN
228                    WEIGHT_K(-1) = (ZC - DES_POS_NEW(L,3))*oDZ_T(K-1)
229                    WEIGHT_K( 0) = ONE - WEIGHT_K(-1)
230                    WEIGHT_K( 1) = ZERO
231                 ELSE
232                    WEIGHT_K( 1) = (DES_POS_NEW(L,3) - ZC)*oDZ_T(K)
233                    WEIGHT_K( 0) = ONE - WEIGHT_K(1)
234                    WEIGHT_K(-1) = ZERO
235                 ENDIF
236              ELSE
237                 Km= 0; Kp=0
238                 WEIGHT_K( 0) = ONE
239              ENDIF
240     
241     ! Set the default fluid cell index and loop counter.
242              IJK = PIJK(L,4)
243              IDX=0
244     
245     ! Calculate weights for ghost particles. Only store weights that the
246     ! current process owns.
247              DO KC=Km,Kp
248              DO IC=-1,+1
249              DO JC=-1,+1
250                 IDX=IDX+1
251                 WEIGHT = WEIGHT_I(IC)*WEIGHT_J(JC)*WEIGHT_K(KC)
252                 IF(IS_ON_MYPE_PLUS2LAYERS(I+IC,J+JC,K+KC)) THEN
253                    IJKt = FUNIJK_MAP_C(I+IC,J+JC,K+KC)
254                    IF(FLUID_AT(IJKt)) THEN
255                       FILTER_CELL(IDX,L) = IJKt
256                       FILTER_WEIGHT(IDX,L) = WEIGHT
257                    ELSE
258                       FILTER_CELL(IDX,L) = IJK
259                       FILTER_WEIGHT(IDX,L) = WEIGHT
260                    ENDIF
261                 ELSE
262                    FILTER_CELL(IDX,L) = IJK
263                    FILTER_WEIGHT(IDX,L) = ZERO
264                 ENDIF
265              ENDDO
266              ENDDO
267              ENDDO
268     
269           ENDDO
270     
271           CONTAINS
272     
273           INCLUDE 'functions.inc'
274     
275           END SUBROUTINE CALC_INTERP_WEIGHTS2
276