File: /nfs/home/0/users/jenkins/mfix.git/model/source_pp_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_Pp_g                                             C
4     !  Purpose: Determine source terms for Pressure correction equation.   C
5     !                                                                      C
6     !  Notes: The off-diagonal coefficients are positive. The center       C
7     !         coefficient and the source vector are negative. See          C
8     !         conv_Pp_g                                                    C
9     !                                                                      C
10     !  Author: M. Syamlal                                 Date: 21-JUN-96  C
11     !  Reviewer:                                          Date:            C
12     !                                                                      C
13     !  Revision Number: 1                                                  C
14     !  Purpose: To incorporate Cartesian grid modifications                C
15     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
16     !                                                                      C
17     !  Literature/Document References:                                     C
18     !  Variables referenced:                                               C
19     !  Variables modified:                                                 C
20     !  Local variables:                                                    C
21     !                                                                      C
22     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
23     
24           SUBROUTINE SOURCE_PP_G(A_M, B_M, B_MMAX)
25     
26     !-----------------------------------------------
27     ! Modules
28     !-----------------------------------------------
29           USE bc
30           USE compar
31           USE constant
32           USE cutcell
33           USE eos, ONLY: DROODP_G
34           USE fldvar
35           USE functions
36           USE geometry
37           USE indices
38           USE matrix
39           USE parallel
40           USE param
41           USE param1
42           USE pgcor
43           USE physprop
44           USE quadric
45           USE run
46           USE rxns
47           USE ur_facs
48           USE vshear
49           Use xsi_array
50     
51           IMPLICIT NONE
52     !-----------------------------------------------
53     ! Dummy arguments
54     !-----------------------------------------------
55     ! Septadiagonal matrix A_m
56           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
57     ! Vector b_m
58           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
59     ! maximum term in b_m expression
60           DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M)
61     !-----------------------------------------------
62     ! Local Variables
63     !-----------------------------------------------
64     ! solids phase index
65           INTEGER :: M
66     ! Indices
67           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
68     ! under relaxation factor for pressure
69           DOUBLE PRECISION fac
70     ! terms of bm expression
71           DOUBLE PRECISION bma, bme, bmw, bmn, bms, bmt, bmb, bmr
72     ! loezos: used for including shearing
73           integer :: incr
74     ! error message
75           CHARACTER(LEN=80) :: LINE(1)
76     ! temporary use of global arrays:
77     ! xsi_array: convection weighting factors
78     !      DOUBLE PRECISION XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),&
79     !                       XSI_t(DIMENSION_3)
80     !-----------------------------------------------
81           call lock_xsi_array
82     ! loezos
83     ! update to true velocity
84           IF (SHEAR) THEN
85     !!$omp parallel do private(IJK)
86              DO IJK = IJKSTART3, IJKEND3
87                 IF (FLUID_AT(IJK)) THEN
88                    V_G(IJK)=V_G(IJK)+VSH(IJK)
89                 ENDIF
90              ENDDO
91           ENDIF
92     
93     ! Calculate convection-diffusion fluxes through each of the faces
94     
95     !$omp parallel default(none) &
96     !$omp          private(IJK, IMJK, IJMK, IJKM, M, bma,bme,bmw,bmn,bms,bmt,bmb,bmr,line)  &
97     !$omp          shared(ijkstart3,ijkend3,cartesian_grid,rop_g,rop_go,rop_s,rop_so,vol,odt,u_g,v_g,w_g,u_s,v_s,w_s,b_m, &
98     !$omp                 b_mmax,d_e,d_n,d_t,a_m,a_upg_e,a_vpg_n,a_wpg_t,mmax,close_packed,sum_r_s,sum_r_g,ro_g0)
99     !$omp do
100           DO IJK = ijkstart3, ijkend3
101              IF (FLUID_AT(IJK)) THEN
102                 IMJK = IM_OF(IJK)
103                 IJMK = JM_OF(IJK)
104                 IJKM = KM_OF(IJK)
105     
106                 bma = (ROP_G(IJK)-ROP_GO(IJK))*VOL(IJK)*ODT
107                 bme = A_M(IJK,E,0)*U_G(IJK)
108                 bmw = A_M(IJK,W,0)*U_G(IMJK)
109                 bmn = A_M(IJK,N,0)*V_G(IJK)
110                 bms = A_M(IJK,S,0)*V_G(IJMK)
111                 bmt = A_M(IJK,T,0)*W_G(IJK)
112                 bmb = A_M(IJK,B,0)*W_G(IJKM)
113                 bmr = SUM_R_G(IJK)*VOL(IJK)
114                 B_M(IJK,0) = -((-(bma + bme - bmw + bmn - bms + bmt - bmb ))+ bmr )
115                 B_MMAX(IJK,0) = max(abs(bma), abs(bme), abs(bmw), abs(bmn), abs(bms), abs(bmt), abs(bmb), abs(bmr) )
116     
117                 A_M(IJK,E,0) = A_M(IJK,E,0)*D_E(IJK,0)
118                 A_M(IJK,W,0) = A_M(IJK,W,0)*D_E(IMJK,0)
119                 A_M(IJK,N,0) = A_M(IJK,N,0)*D_N(IJK,0)
120                 A_M(IJK,S,0) = A_M(IJK,S,0)*D_N(IJMK,0)
121                 A_M(IJK,T,0) = A_M(IJK,T,0)*D_T(IJK,0)
122                 A_M(IJK,B,0) = A_M(IJK,B,0)*D_T(IJKM,0)
123     
124                 IF(CARTESIAN_GRID) THEN
125                    A_M(IJK,E,0) = A_M(IJK,E,0) * A_UPG_E(IJK)
126                    A_M(IJK,W,0) = A_M(IJK,W,0) * A_UPG_E(IMJK)
127                    A_M(IJK,N,0) = A_M(IJK,N,0) * A_VPG_N(IJK)
128                    A_M(IJK,S,0) = A_M(IJK,S,0) * A_VPG_N(IJMK)
129                    A_M(IJK,T,0) = A_M(IJK,T,0) * A_WPG_T(IJK)
130                    A_M(IJK,B,0) = A_M(IJK,B,0) * A_WPG_T(IJKM)
131                 ENDIF
132     
133     ! solids phase terms are needed for those solids phases that do not
134     ! become close packed. these terms are used to essentially make the
135     ! matrix equation for gas pressure correction a mixture pressure
136     ! correction equation by adding solids phases that do not become
137     ! close packed
138                 DO M = 1, MMAX
139                    IF (.NOT.CLOSE_PACKED(M)) THEN
140                       B_M(IJK,0) = B_M(IJK,0) - &
141                          ((-((ROP_S(IJK,M)-ROP_SO(IJK,M))*VOL(IJK)*ODT+&
142                          A_M(IJK,E,M)*U_S(IJK,M)-A_M(IJK,W,M)*U_S(IMJK,M)+&
143                          A_M(IJK,N,M)*V_S(IJK,M)-A_M(IJK,S,M)*V_S(IJMK,M)+&
144                          A_M(IJK,T,M)*W_S(IJK,M)-A_M(IJK,B,M)*W_S(IJKM,M)))+&
145                          SUM_R_S(IJK,M)*VOL(IJK))
146     
147                       IF(.NOT.CARTESIAN_GRID) THEN
148                          A_M(IJK,E,0) = A_M(IJK,E,0) + A_M(IJK,E,M)*D_E(IJK,M)
149                          A_M(IJK,W,0) = A_M(IJK,W,0) + A_M(IJK,W,M)*D_E(IMJK,M)
150                          A_M(IJK,N,0) = A_M(IJK,N,0) + A_M(IJK,N,M)*D_N(IJK,M)
151                          A_M(IJK,S,0) = A_M(IJK,S,0) + A_M(IJK,S,M)*D_N(IJMK,M)
152                          A_M(IJK,T,0) = A_M(IJK,T,0) + A_M(IJK,T,M)*D_T(IJK,M)
153                          A_M(IJK,B,0) = A_M(IJK,B,0) + A_M(IJK,B,M)*D_T(IJKM,M)
154                       ELSE
155                          A_M(IJK,E,0) = A_M(IJK,E,0) + A_M(IJK,E,M)*D_E(IJK,M)  * A_UPG_E(IJK)
156                          A_M(IJK,W,0) = A_M(IJK,W,0) + A_M(IJK,W,M)*D_E(IMJK,M) * A_UPG_E(IMJK)
157                          A_M(IJK,N,0) = A_M(IJK,N,0) + A_M(IJK,N,M)*D_N(IJK,M)  * A_VPG_N(IJK)
158                          A_M(IJK,S,0) = A_M(IJK,S,0) + A_M(IJK,S,M)*D_N(IJMK,M) * A_VPG_N(IJMK)
159                          A_M(IJK,T,0) = A_M(IJK,T,0) + A_M(IJK,T,M)*D_T(IJK,M)  * A_WPG_T(IJK)
160                          A_M(IJK,B,0) = A_M(IJK,B,0) + A_M(IJK,B,M)*D_T(IJKM,M) * A_WPG_T(IJKM)
161                       ENDIF
162                    ENDIF
163                 ENDDO
164     
165                 A_M(IJK,0,0) = -(A_M(IJK,E,0)+A_M(IJK,W,0)+A_M(IJK,N,0)+A_M(IJK,S,0&
166                    )+A_M(IJK,T,0)+A_M(IJK,B,0))
167     
168                 IF (ABS(A_M(IJK,0,0)) < SMALL_NUMBER) THEN
169                    IF (ABS(B_M(IJK,0)) < SMALL_NUMBER) THEN
170                       A_M(IJK,0,0) = -ONE
171                       B_M(IJK,0) = ZERO
172                    ELSEIF (RO_G0 .NE. UNDEFINED) THEN !This is an error only in incompressible flow
173     !!$omp             critical
174                       WRITE (LINE, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', IJK, &
175                          ' M = ', 0, ' A = 0 and b = ', B_M(IJK,0)
176                       CALL WRITE_ERROR ('SOURCE_Pp_g', LINE, 1)
177     !!$omp             end critical
178                    ENDIF
179                 ENDIF
180     
181              ELSE   ! if/else branch .not.fluid_at(ijk)
182     ! set the value (correction) in all wall and flow boundary cells to zero
183     ! note the matrix coefficients and source vector should already be zero
184     ! from the initialization of A and B but the following ensures the zero
185     ! value.
186                 A_M(IJK,E,0) = ZERO
187                 A_M(IJK,W,0) = ZERO
188                 A_M(IJK,N,0) = ZERO
189                 A_M(IJK,S,0) = ZERO
190                 A_M(IJK,T,0) = ZERO
191                 A_M(IJK,B,0) = ZERO
192                 A_M(IJK,0,0) = -ONE
193                 B_M(IJK,0) = ZERO
194              ENDIF   ! end if/else branch fluid_at(ijk)
195           ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)
196     !$omp end parallel
197     
198     
199     ! loezos
200           IF (SHEAR) THEN
201     !!$omp parallel do private(IJK)
202              DO IJK = IJKSTART3, IJKEND3
203                 IF (FLUID_AT(IJK)) THEN
204                    V_G(IJK)=V_G(IJK)-VSH(IJK)
205                 ENDIF
206              ENDDO
207           ENDIF
208     
209     
210     ! make correction for compressible flows
211           IF (RO_G0 == UNDEFINED) THEN
212              fac = UR_FAC(1)  !since p_g = p_g* + ur_fac * pp_g
213     
214     ! loezos
215              incr=0
216     !         CALL CALC_XSI(DISCRETIZE(1),ROP_G,U_G,V_G,W_G,XSI_E,XSI_N,XSI_T,incr)
217     
218     !!$omp    parallel do                                                     &
219     !!$omp&   private(IJK,I,J,K,                                       &
220     !!$omp&            IMJK,IJMK,IJKM,IJKE,IJKW,IJKN,IJKS,IJKT,IJKB)
221              DO IJK = ijkstart3, ijkend3
222                 IF (FLUID_AT(IJK)) THEN
223     
224                    A_M(IJK,0,0) = A_M(IJK,0,0) - &
225                       fac*DROODP_G(RO_G(IJK),P_G(IJK))*&
226                       EP_G(IJK)*VOL(IJK)*ODT
227     
228     ! Although the following is a better approximation for high speed flows because
229     ! it considers density changes in the neighboring cells, the code runs faster
230     ! without it for low speed flows.  The gas phase mass balance cannot be
231     ! maintained to machine precision with the following approximation. If the
232     ! following lines are uncommented, the calc_xsi call above should also be
233     ! uncommented
234     !               IMJK = IM_OF(IJK)
235     !               IJMK = JM_OF(IJK)
236     !               IJKM = KM_OF(IJK)
237     !               IJKE = EAST_OF(IJK)
238     !               IJKW = WEST_OF(IJK)
239     !               IJKN = NORTH_OF(IJK)
240     !               IJKS = SOUTH_OF(IJK)
241     !               IJKT = TOP_OF(IJK)
242     !               IJKB = BOTTOM_OF(IJK)
243     !               A_M(IJK,0,0) = A_M(IJK,0,0) - fac*DROODP_G(RO_G(IJK),P_G(IJK))*EP_G(&
244     !                  IJK)*((ONE - XSI_E(IJK))*U_G(IJK)*AYZ(IJK)-XSI_E(IMJK)*U_G(&
245     !                  IMJK)*AYZ(IMJK)+(ONE-XSI_N(IJK))*V_G(IJK)*AXZ(IJK)-XSI_N(IJMK&
246     !                  )*V_G(IJMK)*AXZ(IJMK))
247     
248     !               A_M(IJK,E,0) = A_M(IJK,E,0) - EP_G(IJKE)*fac*DROODP_G(RO_G(IJKE),P_G&
249     !                  (IJKE))*XSI_E(IJK)*U_G(IJK)*AYZ(IJK)
250     !               A_M(IJK,W,0) = A_M(IJK,W,0) + EP_G(IJKW)*fac*DROODP_G(RO_G(IJKW),P_G&
251     !                  (IJKW))*(ONE - XSI_E(IMJK))*U_G(IMJK)*AYZ(IMJK)
252     !               A_M(IJK,N,0) = A_M(IJK,N,0) - EP_G(IJKN)*fac*DROODP_G(RO_G(IJKN),P_G&
253     !                  (IJKN))*XSI_N(IJK)*V_G(IJK)*AXZ(IJK)
254     !               A_M(IJK,S,0) = A_M(IJK,S,0) + EP_G(IJKS)*fac*DROODP_G(RO_G(IJKS),P_G&
255     !                  (IJKS))*(ONE - XSI_N(IJMK))*V_G(IJMK)*AXZ(IJMK)
256     !               IF (DO_K) THEN
257     !                  A_M(IJK,0,0) = A_M(IJK,0,0) - fac*DROODP_G(RO_G(IJK),P_G(IJK))*&
258     !                     EP_G(IJK)*((ONE - XSI_T(IJK))*W_G(IJK)*AXY(IJK)-XSI_T(IJKM&
259     !                     )*W_G(IJKM)*AXY(IJKM))
260     !                  A_M(IJK,T,0) = A_M(IJK,T,0) - EP_G(IJKT)*fac*DROODP_G(RO_G(IJKT),&
261     !                     P_G(IJKT))*XSI_T(IJK)*W_G(IJK)*AXY(IJK)
262     !                  A_M(IJK,B,0) = A_M(IJK,B,0) + EP_G(IJKB)*fac*DROODP_G(RO_G(IJKB),&
263     !                     P_G(IJKB))*(ONE - XSI_T(IJKM))*W_G(IJKM)*AXY(IJKM)
264     !               ENDIF
265     !
266                 ENDIF   !end if (fluid_at(ijk))
267              ENDDO    ! end do (ijk=ijkstart3,ijkend3)
268           ENDIF   ! end if (ro_g0 == undefined); i.e., compressible flow
269     
270     
271     ! Remove the asymmetry in matrix caused by the pressure outlet or inlet
272     ! boundaries.  Because the P' at such boundaries is zero we may set the
273     ! coefficient in the neighboring fluid cell to zero without affecting
274     ! the linear equation set.
275     !!$omp    parallel do                                                     &
276     !!$omp&   private(IJK,IMJK, IPJK, IJMK, IJPK, IJKM, IJKP)
277           DO IJK = ijkstart3, ijkend3
278              IF (FLUID_AT(IJK)) THEN
279                 IMJK = IM_OF(IJK)
280                 IPJK = IP_OF(IJK)
281                 IJMK = JM_OF(IJK)
282                 IJPK = JP_OF(IJK)
283                 IJKM = KM_OF(IJK)
284                 IJKP = KP_OF(IJK)
285     ! Cutting the neighbor link between fluid cell and adjacent p_flow_at cell
286                 if(p_flow_at(imjk)) A_m(IJK, w, 0) = ZERO
287                 if(p_flow_at(ipjk)) A_m(IJK, e, 0) = ZERO
288                 if(p_flow_at(ijmk)) A_m(IJK, s, 0) = ZERO
289                 if(p_flow_at(ijpk)) A_m(IJK, n, 0) = ZERO
290                 if(p_flow_at(ijkm)) A_m(IJK, b, 0) = ZERO
291                 if(p_flow_at(ijkp)) A_m(IJK, t, 0) = ZERO
292              ENDIF
293           ENDDO
294     
295     ! Specify P' to zero for incompressible flows. Check set_bc0
296     ! for details on selection of IJK_P_g.
297           IF (IJK_P_G /= UNDEFINED_I) THEN
298              B_M(IJK_P_G,0) = ZERO
299              A_M(IJK_P_G,:,0) = ZERO
300              A_M(IJK_P_G,0,0) = -ONE
301           ENDIF
302     
303           call unlock_xsi_array
304     
305           RETURN
306           END SUBROUTINE SOURCE_PP_G
307