File: RELATIVE:/../../../mfix.git/model/source_phi.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_phi                                              C
4     !  Purpose: Determine source terms for phi eq. The terms               C
5     !     appear in the center coefficient and RHS vector. The center      C
6     !     coefficient and source vector are negative. The off-diagonal     C
7     !     coefficients are positive. S_p must be positive.                 C
8     !                                                                      C
9     !  Note: This routine is now restricted to Non-Negative scalers when   C
10     !     using deferred correction.                                       C
11     !                                                                      C
12     !                                                                      C
13     !  Author: M. Syamlal                                 Date: 30-APR-97  C
14     !  Reviewer:                                          Date:            C
15     !                                                                      C
16     !                                                                      C
17     !  Literature/Document References:                                     C
18     !                                                                      C
19     !  Variables referenced:                                               C
20     !  Variables modified:                                                 C
21     !  Local variables:                                                    C
22     !                                                                      C
23     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
24     
25           SUBROUTINE SOURCE_PHI(S_P, S_C, EP, PHI, M, A_M, B_M)
26     
27     !-----------------------------------------------
28     ! Modules
29     !-----------------------------------------------
30           USE param
31           USE param1
32           USE parallel
33           USE matrix
34           USE scales
35           USE physprop
36           USE fldvar
37           USE visc_s
38           USE rxns
39           USE run
40           USE toleranc
41           USE geometry
42           USE indices
43           USE is
44           USE tau_s
45           USE compar
46           USE fun_avg
47           USE functions
48           IMPLICIT NONE
49     !-----------------------------------------------
50     ! Dummy Arguments
51     !-----------------------------------------------
52     ! Source term on LHS.  Must be positive.
53           DOUBLE PRECISION, INTENT(IN) :: S_p(DIMENSION_3)
54     ! Source term on RHS
55           DOUBLE PRECISION, INTENT(IN) :: S_C(DIMENSION_3)
56     ! Phase volume fraction
57           DOUBLE PRECISION, INTENT(IN) :: EP(DIMENSION_3)
58     ! Phi
59           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
60     ! phase index
61           INTEGER, INTENT(IN) :: M
62     ! Septadiagonal matrix A_m
63           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
64     ! Vector b_m
65           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
66     !-----------------------------------------------
67     ! Local variables
68     !-----------------------------------------------
69     ! Indices
70           INTEGER :: IJK, IMJK, IJMK, IJKM
71     !-----------------------------------------------
72     
73           DO IJK = ijkstart3, ijkend3
74     
75              IF (FLUID_AT(IJK)) THEN
76     
77     ! dilute flow
78                 IF (EP(IJK) <= DIL_EP_S) THEN
79                    A_M(IJK,E,M) = ZERO
80                    A_M(IJK,W,M) = ZERO
81                    A_M(IJK,N,M) = ZERO
82                    A_M(IJK,S,M) = ZERO
83                    A_M(IJK,T,M) = ZERO
84                    A_M(IJK,B,M) = ZERO
85                    A_M(IJK,0,M) = -ONE
86                    B_M(IJK,M) = ZERO
87     
88     ! use the average boundary cell values to compute phi (sof, Aug 23 2005)
89                    IMJK = IM_OF(IJK)
90                    IJMK = JM_OF(IJK)
91                    IJKM = KM_OF(IJK)
92                    IF (EP(WEST_OF(IJK)) > DIL_EP_S .AND. &
93                        .NOT.IS_AT_E(IMJK)) A_M(IJK,W,M) = ONE
94                    IF (EP(EAST_OF(IJK)) > DIL_EP_S .AND. &
95                        .NOT.IS_AT_E(IJK)) A_M(IJK,E,M) = ONE
96                    IF (EP(SOUTH_OF(IJK)) > DIL_EP_S .AND. &
97                        .NOT.IS_AT_N(IJMK)) A_M(IJK,S,M) = ONE
98                    IF (EP(NORTH_OF(IJK)) > DIL_EP_S .AND. &
99                        .NOT.IS_AT_N(IJK)) A_M(IJK,N,M) = ONE
100                    IF(.NOT. NO_K) THEN
101                       IF (EP(BOTTOM_OF(IJK)) > DIL_EP_S .AND. &
102                          .NOT.IS_AT_T(IJKM)) A_M(IJK,B,M) = ONE
103                       IF (EP(TOP_OF(IJK)) > DIL_EP_S .AND. &
104                          .NOT.IS_AT_T(IJK)) A_M(IJK,T,M) = ONE
105                    ENDIF
106     
107                    IF((A_M(IJK,W,M)+A_M(IJK,E,M)+&
108                       A_M(IJK,S,M)+A_M(IJK,N,M)+ &
109                       A_M(IJK,B,M)+A_M(IJK,T,M)) == ZERO) THEN
110                       B_M(IJK,M) = -PHI(IJK)
111                    ELSE
112                       A_M(IJK,0,M) = -(A_M(IJK,E,M) + A_M(IJK,W,M) +&
113                          A_M(IJK,N,M) + A_M(IJK,S,M) + &
114                          A_M(IJK,T,M)+A_M(IJK,B,M))
115                    ENDIF
116     
117     ! Normal case
118                 ELSE
119     
120     ! Collect the terms
121                    A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+&
122                                     A_M(IJK,N,M)+A_M(IJK,S,M)+&
123                                     A_M(IJK,T,M)+A_M(IJK,B,M)+S_P(IJK))
124     
125     ! B_m and A_m are corrected in case deferred corrections computes B_m > S_c
126     ! see CONV_DIF_PHI_DC.
127                    IF(B_M(IJK,M) < S_C(IJK) .OR. PHI(IJK) == ZERO) THEN
128                       B_M(IJK,M) = -S_C(IJK)+B_M(IJK,M)
129                    ELSE ! disable ELSE statememt if PHI can be negative
130                       A_M(IJK,0,M) = A_M(IJK,0,M) - B_M(IJK,M)/PHI(IJK)
131                       B_M(IJK,M) = -S_C(IJK)
132                    ENDIF
133     
134                 ENDIF   ! end if/else (ep_g(ijk)<=dil_ep_s)
135              ENDIF    ! end if (fluid_at(ijk))
136           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
137     
138     
139           RETURN
140           END SUBROUTINE SOURCE_PHI
141     
142     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
143     !                                                                      C
144     !  Subroutine: POINT_SOURCE_PHI                                        C
145     !                                                                      C
146     !  Purpose: Adds point sources to the scalar equations.                C
147     !                                                                      C
148     !  Author: J. Musser                                  Date: 10-JUN-13  C
149     !  Reviewer:                                          Date:            C
150     !                                                                      C
151     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
152           SUBROUTINE POINT_SOURCE_PHI(PHI, PS_PHI, PS_FLOW,  &
153              M, A_M, B_M)
154     
155           use compar
156           use geometry
157           use indices
158           use physprop
159           use ps
160           use run
161     
162     ! To be removed upon complete integration of point source routines.
163           use bc
164           use usr
165           use functions
166     
167           IMPLICIT NONE
168     !-----------------------------------------------
169     ! Dummy arguments
170     !-----------------------------------------------
171     
172     ! Vector b_m
173           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
174     
175     
176     ! maximum term in b_m expression
177           DOUBLE PRECISION, INTENT(IN) :: PS_PHI(DIMENSION_BC)
178     
179     ! maximum term in b_m expression
180           DOUBLE PRECISION, INTENT(IN) :: PS_FLOW(DIMENSION_BC)
181     
182           INTEGER, intent(in) :: M
183     
184     ! Septadiagonal matrix A_m
185           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
186     
187     ! Vector b_m
188           DOUBLE PRECISION, intent(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
189     
190     !-----------------------------------------------
191     ! Local Variables
192     !-----------------------------------------------
193     
194     ! Indices
195           INTEGER :: IJK, I, J, K
196           INTEGER :: PSV
197     
198     ! terms of bm expression
199           DOUBLE PRECISION pSource
200     
201     !-----------------------------------------------
202     
203           PSV_LP: do PSV = 1, DIMENSION_PS
204     
205              if(.NOT.PS_DEFINED(PSV)) cycle PSV_LP
206              if(abs(PS_FLOW(PSV)) < small_number) cycle PSV_LP
207     
208              do k = PS_K_B(PSV), PS_K_T(PSV)
209              do j = PS_J_S(PSV), PS_J_N(PSV)
210              do i = PS_I_W(PSV), PS_I_E(PSV)
211     
212                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
213                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
214     
215                 ijk = funijk(i,j,k)
216                 if(.NOT.fluid_at(ijk)) cycle
217     
218                 if(A_M(IJK,0,M) == -ONE .AND. B_M(IJK,M) == -PHI(IJK)) then
219                    B_M(IJK,M) = -PS_PHI(PSV)
220                 else
221     
222     ! Calculate the mass flow rate for this cell. (mass/time)
223                    pSource = PS_FLOW(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
224     
225                    A_M(IJK,0,M) = A_M(IJK,0,M) - pSource
226                    B_M(IJK,M) = B_M(IJK,M) - PS_PHI(PSV) * pSource
227     
228                 endif
229     
230              enddo
231              enddo
232              enddo
233     
234           enddo PSV_LP
235     
236           RETURN
237           END SUBROUTINE POINT_SOURCE_PHI
238