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