File: /nfs/home/0/users/jenkins/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     
157           use geometry
158           use indices
159           use physprop
160           use ps
161           use run
162     
163     ! To be removed upon complete integration of point source routines.
164           use bc
165           use usr
166           use functions
167     
168           IMPLICIT NONE
169     !-----------------------------------------------
170     ! Dummy arguments
171     !-----------------------------------------------
172     
173     ! Vector b_m
174           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
175     
176     
177     ! maximum term in b_m expression
178           DOUBLE PRECISION, INTENT(IN) :: PS_PHI(DIMENSION_BC)
179     
180     ! maximum term in b_m expression
181           DOUBLE PRECISION, INTENT(IN) :: PS_FLOW(DIMENSION_BC)
182     
183           INTEGER, intent(in) :: M
184     
185     ! Septadiagonal matrix A_m
186           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
187     
188     ! Vector b_m
189           DOUBLE PRECISION, intent(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
190     
191     !-----------------------------------------------
192     ! Local Variables
193     !-----------------------------------------------
194     
195     ! Indices
196           INTEGER :: IJK, I, J, K
197           INTEGER :: PSV
198     
199     ! terms of bm expression
200           DOUBLE PRECISION pSource
201     
202     !-----------------------------------------------
203     
204     ! External function. Integrates the temperature-dependent specific
205     ! heat from zero to T.
206           DOUBLE PRECISION, EXTERNAL :: calc_ICpoR
207     
208           PSV_LP: do PSV = 1, DIMENSION_PS
209     
210              if(.NOT.PS_DEFINED(PSV)) cycle PSV_LP
211              if(abs(PS_FLOW(PSV)) < small_number) cycle PSV_LP
212     
213              do k = PS_K_B(PSV), PS_K_T(PSV)
214              do j = PS_J_S(PSV), PS_J_N(PSV)
215              do i = PS_I_W(PSV), PS_I_E(PSV)
216     
217                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
218                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
219     
220                 ijk = funijk(i,j,k)
221                 if(.NOT.fluid_at(ijk)) cycle
222     
223                 if(A_M(IJK,0,M) == -ONE .AND. B_M(IJK,M) == -PHI(IJK)) then
224                    B_M(IJK,M) = -PS_PHI(PSV)
225                 else
226     
227     ! Calculate the mass flow rate for this cell. (mass/time)
228                    pSource = PS_FLOW(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
229     
230                    A_M(IJK,0,M) = A_M(IJK,0,M) - pSource
231                    B_M(IJK,M) = B_M(IJK,M) - PS_PHI(PSV) * pSource
232     
233                 endif
234     
235              enddo
236              enddo
237              enddo
238     
239           enddo PSV_LP
240     
241     
242           RETURN
243           END SUBROUTINE POINT_SOURCE_PHI
244