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