File: N:\mfix\model\des\pic\calc_ps_grad_pic.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: CALC_PS_GRAD_PIC                                        !
3     !  Author: R. Garg                                                     !
4     !                                                                      !
5     !  Purpose: Calculate the solids pressure graident.                    !
6     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
7           SUBROUTINE CALC_PS_GRAD_PIC
8     
9     ! Modules
10     !---------------------------------------------------------------------//
11           use particle_filter, only: DES_INTERP_SCHEME_ENUM
12           use particle_filter, only: DES_INTERP_GARG
13           use mfix_pic, only: PIC_P_S
14           use mfix_pic, only: PS_FORCE_PIC
15           IMPLICIT NONE
16     
17     !......................................................................!
18     
19           SELECT CASE(DES_INTERP_SCHEME_ENUM)
20           CASE(DES_INTERP_GARG); CALL CALC_PS_GRAD_PIC0
21           CASE DEFAULT; CALL CALC_GRAD_DES(PIC_P_S(:,1), PS_FORCE_PIC)
22           END SELECT
23     
24           RETURN
25           END SUBROUTINE CALC_PS_GRAD_PIC
26     
27     
28     
29     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
30     !  Subroutine: CALC_PS_GRAD_PIC0                                       !
31     !  Author: R. Garg                                                     !
32     !                                                                      !
33     !  Purpose: Calculate the solids pressure graident. This routine       !
34     !  stores the solid pressure graident at cell faces.                   !
35     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
36           SUBROUTINE CALC_PS_GRAD_PIC0
37     
38     ! Modules
39     !---------------------------------------------------------------------//
40           use compar, only: IJKSTART3, IJKEND3
41           use compar, only: istart2, jstart2, kstart2
42           use compar, only: istart3, jstart3, kstart3
43           use compar, only: iend3, jend3, kend3
44           use functions, only: FLUID_AT
45           use indices, only: I_OF, J_OF, K_OF
46           use functions, only: IP_OF, JP_OF, KP_OF
47           use functions, only: funijk
48           use functions, only: is_on_mype_owns
49           use geometry, only: DO_K, NO_K
50           use geometry, only: DX, DY, DZ
51           use geometry, only: imin2, jmin2, kmin2
52           USE mfix_pic, only: pic_p_s, ps_force_pic
53           use param1, only: ZERO
54           use sendrecv, only: send_recv
55           implicit none
56     
57     ! Local variables
58     !---------------------------------------------------------------------//
59     ! general i, j, k indices
60           INTEGER I, J, K, IJK, IPJK, IJPK, IJKP, IDIM
61           integer :: I1, J1, K1
62     !......................................................................!
63     
64     ! Since EP_G is already shared across the processors, the pressure
65     ! gradient calculation can be made a function call so that the extra
66     ! communication of P_S can be avoided.
67     
68           DO IJK = IJKSTART3, IJKEND3
69              PS_FORCE_PIC(:,IJK) = ZERO
70              IF(.NOT.FLUID_AT(IJK)) CYCLE
71     
72              I = I_OF(IJK)
73              J = J_OF(IJK)
74              K = K_OF(IJK)
75     
76              IPJK = IP_OF(IJK)
77              IJPK = JP_OF(IJK)
78              IJKP = KP_OF(IJK)
79     
80     ! Calculate the solids pressure gradient at east face.
81              IF(FLUID_AT(IPJK)) THEN
82                 PS_FORCE_PIC(1,IJK) = 2.0d0 *                              &
83                    (PIC_P_S(IPJK,1) - PIC_P_S(IJK,1)) /                    &
84                    (DX(I) + DX(I_OF(IPJK)))
85              ELSE
86                 IF(PIC_P_S(IJK,1) > ZERO) THEN
87                    PS_FORCE_PIC(1,IJK) = 2.0d0*                            &
88                       (PIC_P_S(IPJK,1) - PIC_P_S(IJK,1)) /                 &
89                       (DX(I) + DX(I_OF(IPJK)))
90                 ELSE
91                    PS_FORCE_PIC(1,IJK) = ZERO
92                 ENDIF
93              ENDIF
94     
95     ! Calculate the solids pressure graident at the north face.
96              IF(FLUID_AT(IJPK)) THEN
97                 PS_FORCE_PIC(2,IJK) = 2.0d0*                               &
98                    (PIC_P_S(IJPK,1) - PIC_P_S(IJK,1)) /                    &
99                    (DY(J)+DY(J_OF(IJPK)))
100              ELSE
101                 IF(PIC_P_S(IJK,1) > ZERO) THEN
102                    PS_FORCE_PIC(2,IJK) = 2.0d0*                            &
103                       (PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/                  &
104                       (DY(j)+Dy(j_of(ijpk)))
105                 ELSE
106                    PS_FORCE_PIC(2,IJK) = ZERO
107                 ENDIF
108              ENDIF
109     
110     ! Calculate the solids pressure graident at the top face.
111              IF(DO_K) THEN
112                 IF(FLUID_AT(IJKP)) THEN
113                    PS_FORCE_PIC(3,IJK) = 2.0d0*                            &
114                       (PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/                  &
115                       (DZ(K)+DZ(K_OF(IJKP)))
116                 ELSE
117                    IF(PIC_P_S(IJK,1).GT.ZERO) then
118                       PS_FORCE_PIC(3,IJK) = 2.0d0*&
119                          (PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/               &
120                          (DZ(K)+DZ(K_OF(IJKP)))
121                    ELSE
122                       PS_FORCE_PIC(3,IJK) = ZERO
123                    ENDIF
124                 ENDIF
125              ENDIF
126           ENDDO
127     
128     ! Compute the pressure gradients along the west domain bondary which
129     ! is previously skipped.
130           I1 = IMIN2
131           IF(I1 == ISTART2) THEN
132              DO K1 = KSTART3, KEND3
133              DO J1 = JSTART3, JEND3
134                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
135                 IJK = FUNIJK(I1,J1,K1)
136                 IPJK = IP_OF(IJK)
137                 IF(PIC_P_S(IPJK,1) > ZERO) THEN
138                    PS_FORCE_PIC(1,IJK) = 2.0d0*                            &
139                       (PIC_P_S(IPJK,1) - PIC_P_S(IJK,1)) /                 &
140                       (DX(I1) + DX(I_OF(IPJK)))
141                 ELSE
142                    PS_FORCE_PIC(1,IJK) = ZERO
143                 ENDIF
144              ENDDO
145              ENDDO
146           ENDIF
147     
148     ! Compute the pressure gradients along the south domain bondary which
149     ! is previously skipped.
150           J1 = JMIN2
151           IF(J1 == JSTART2) THEN
152              DO K1 = KSTART3, KEND3
153              DO I1 = ISTART3, IEND3
154                 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
155                 IJK = FUNIJK(I1,J1,K1)
156                 IJPK = JP_OF(IJK)
157                 IF(PIC_P_S(IJPK,1).GT.ZERO) then
158                    PS_FORCE_PIC(2,IJK) = 2.0d0*                            &
159                       (PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/                  &
160                       (DY(J) + DY(J_OF(IJPK)))
161                 ELSE
162                    PS_FORCE_PIC(2,IJK) = ZERO
163                 ENDIF
164              ENDDO
165              ENDDO
166           ENDIF
167     
168     ! Compute the pressure gradients along the south domain bondary which
169     ! is previously skipped.
170           IF(DO_K) then
171              K1 = KMIN2
172              IF(K1 == KSTART2) THEN
173                 DO J1 = JSTART3, JEND3
174                 DO I1 = ISTART3, IEND3
175                    IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
176                    IJK = FUNIJK(I1,J1,K1)
177                    IJKP = KP_OF(IJK)
178                    IF(PIC_P_S(IJKP,1).GT.ZERO) THEN
179                       PS_FORCE_PIC(3,IJK) = 2.0d0*                         &
180                          (PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/               &
181                          (DZ(K)+DZ(K_OF(IJKP)))
182                    ELSE
183                       PS_FORCE_PIC(3,IJK) = ZERO
184                    ENDIF
185                 ENDDO
186                 ENDDO
187              ENDIF
188           ENDIF
189     
190           DO IDIM = 1, merge(2,3,NO_K)
191              CALL SEND_RECV(PS_FORCE_PIC(IDIM,:),1)
192           ENDDO
193     
194           END SUBROUTINE CALC_PS_GRAD_PIC0
195