File: /nfs/home/0/users/jenkins/mfix.git/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     
13           SELECT CASE(DES_INTERP_SCHEME_ENUM)
14           CASE(DES_INTERP_DPVM);  CALL CALC_INTERP_WEIGHTS1
15           CASE(DES_INTERP_GAUSS); CALL CALC_INTERP_WEIGHTS1
16           END SELECT
17     
18           RETURN
19           END SUBROUTINE CALC_INTERP_WEIGHTS
20     
21     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
22     !  Subroutine: CALC_INTERP_WEIGHTS0                                    !
23     !                                                                      !
24     !  Purpose: Calculate weights used for interpolation.                  !
25     !                                                                      !
26     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
27           SUBROUTINE CALC_INTERP_WEIGHTS1
28     
29           use discretelement, only: MAX_RADIUS
30           use discretelement, only: PEA, MAX_PIP
31           use discretelement, only: PIJK
32           use discretelement, only: DES_POS_NEW
33           use discretelement, only: XE, YN, ZT
34           use particle_filter, only: FILTER_CELL
35           use particle_filter, only: FILTER_WEIGHT
36           use geometry, only: DO_K
37           use functions, only: FUNIJK
38           use functions, only: FLUID_AT
39           use functions, only: CYCLIC_AT
40           use functions, only: IS_ON_MYPE_PLUS2LAYERS
41     
42           use param1, only: ZERO, ONE
43     
44           IMPLICIT NONE
45     
46           INTEGER :: L, IDX
47           INTEGER :: IJK, IJKt
48     
49           INTEGER :: I, J, K
50           INTEGER :: IC, JC, KC
51           INTEGER :: Km, Kp
52     
53           DOUBLE PRECISION :: WEIGHT
54           DOUBLE PRECISION :: WEIGHT_I(-1:1)
55           DOUBLE PRECISION :: WEIGHT_J(-1:1)
56           DOUBLE PRECISION :: WEIGHT_K(-1:1)
57     
58     !$omp parallel default(none) private(L, WEIGHT_I, WEIGHT_J, WEIGHT_K,  &
59     !$omp    KM, KP, IDX, IJK, I, J, K, WEIGHT, IJKT)                      &
60     !$omp shared(MAX_PIP, PEA, PIJK, DES_POS_NEW, XE, YN, ZT, DO_K,        &
61     !$omp    FILTER_CELL, FILTER_WEIGHT)
62     !$omp do
63           DO L = 1, MAX_PIP
64     
65              IF(.NOT.PEA(L,1)) CYCLE
66              IF(any(PEA(L,2:3))) CYCLE
67     
68              I = PIJK(L,1)
69              J = PIJK(L,2)
70              K = PIJK(L,3)
71     
72     ! Tentative weights for I indices to the West and East.
73              WEIGHT_I(-1) = CALC_FILTER_WEIGHTS(DES_POS_NEW(1,L), XE(I-1))
74              WEIGHT_I( 1) = CALC_FILTER_WEIGHTS(XE(I), DES_POS_NEW(1,L))
75              WEIGHT_I( 0) = ONE - WEIGHT_I(-1) - WEIGHT_I(1)
76     
77     ! Tentative weights for J indices to the South and North.
78              WEIGHT_J(-1) = CALC_FILTER_WEIGHTS(DES_POS_NEW(2,L), YN(J-1))
79              WEIGHT_J( 1) = CALC_FILTER_WEIGHTS(YN(J), DES_POS_NEW(2,L))
80              WEIGHT_J( 0) = ONE - WEIGHT_J(-1) - WEIGHT_J(1)
81     
82     ! Tentative weights for K indices to the Top and Bottom.
83              IF(DO_K) THEN
84                 Km=-1;  Kp=1
85                 WEIGHT_K(-1) = CALC_FILTER_WEIGHTS(DES_POS_NEW(3,L),ZT(K-1))
86                 WEIGHT_K( 1) = CALC_FILTER_WEIGHTS(ZT(K), DES_POS_NEW(3,L))
87                 WEIGHT_K( 0) = ONE - WEIGHT_K(-1) - WEIGHT_K(1)
88              ELSE
89                 Km= 0; Kp=0
90                 WEIGHT_K( 0) = ONE
91              ENDIF
92     
93     ! Set the default fluid cell index and loop counter.
94              IJK = PIJK(L,4)
95              IDX=0
96     
97     ! Calculate weights for ghost particles. Only store weights that the
98     ! current process owns.
99              DO KC=Km,Kp
100              DO IC=-1,+1
101              DO JC=-1,+1
102                 IDX=IDX+1
103                 WEIGHT = WEIGHT_I(IC)*WEIGHT_J(JC)*WEIGHT_K(KC)
104                 IF(IS_ON_MYPE_PLUS2LAYERS(I+IC,J+JC,K+KC)) THEN
105                    IJKt = FUNIJK(I+IC,J+JC,K+KC)
106                    IF(FLUID_AT(IJKt) .OR. CYCLIC_AT(IJKt)) THEN
107                       FILTER_CELL(IDX,L) = IJKt
108                       FILTER_WEIGHT(IDX,L) = WEIGHT
109                    ELSE
110                       FILTER_CELL(IDX,L) = IJK
111                       FILTER_WEIGHT(IDX,L) = WEIGHT
112                    ENDIF
113                 ELSE
114                    FILTER_CELL(IDX,L) = IJK
115                    FILTER_WEIGHT(IDX,L) = ZERO
116                 ENDIF
117              ENDDO
118              ENDDO
119              ENDDO
120     
121           ENDDO
122     !$omp end do
123     !$omp end parallel
124     
125           CONTAINS
126     
127     !``````````````````````````````````````````````````````````````````````!
128     !                                                                      !
129     !``````````````````````````````````````````````````````````````````````!
130           DOUBLE PRECISION FUNCTION CALC_FILTER_WEIGHTS(POS1, POS2)
131     
132           use particle_filter, only: FILTER_WIDTH_INTERP
133           use particle_filter, only: OoFILTER_VOL
134           use particle_filter, only: FILTER_WIDTH_INTERPx3
135     
136           DOUBLE PRECISION, INTENT(IN) :: POS1, POS2
137           DOUBLE PRECISION :: OVERLAP, HEIGHT
138     
139           OVERLAP = POS1 - POS2
140           IF(OVERLAP < FILTER_WIDTH_INTERP) THEN
141              HEIGHT = FILTER_WIDTH_INTERP - OVERLAP
142              CALC_FILTER_WEIGHTS = HEIGHT**2 * &
143                 (FILTER_WIDTH_INTERPx3 - HEIGHT)*OoFILTER_VOL
144           ELSE
145              CALC_FILTER_WEIGHTS = ZERO
146           ENDIF
147     
148           END FUNCTION CALC_FILTER_WEIGHTS
149     
150           END SUBROUTINE CALC_INTERP_WEIGHTS1
151