File: N:\mfix\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 scales
34           USE physprop
35           USE fldvar
36           USE visc_s
37           USE rxns
38           USE run
39           USE toleranc
40           USE geometry
41           USE indices
42           USE is
43           USE tau_s
44           USE compar
45           USE fun_avg
46           USE functions
47           IMPLICIT NONE
48     !-----------------------------------------------
49     ! Dummy Arguments
50     !-----------------------------------------------
51     ! Source term on LHS.  Must be positive.
52           DOUBLE PRECISION, INTENT(IN) :: S_p(DIMENSION_3)
53     ! Source term on RHS
54           DOUBLE PRECISION, INTENT(IN) :: S_C(DIMENSION_3)
55     ! Phase volume fraction
56           DOUBLE PRECISION, INTENT(IN) :: EP(DIMENSION_3)
57     ! Phi
58           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
59     ! phase index
60           INTEGER, INTENT(IN) :: M
61     ! Septadiagonal matrix A_m
62           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
63     ! Vector b_m
64           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
65     !-----------------------------------------------
66     ! Local variables
67     !-----------------------------------------------
68     ! Indices
69           INTEGER :: IJK, IMJK, IJMK, IJKM
70     !-----------------------------------------------
71     
72           DO IJK = ijkstart3, ijkend3
73     
74              IF (FLUID_AT(IJK)) THEN
75     
76     ! dilute flow
77                 IF (EP(IJK) <= DIL_EP_S) THEN
78                    A_M(IJK,east,M) = ZERO
79                    A_M(IJK,west,M) = ZERO
80                    A_M(IJK,north,M) = ZERO
81                    A_M(IJK,south,M) = ZERO
82                    A_M(IJK,top,M) = ZERO
83                    A_M(IJK,bottom,M) = ZERO
84                    A_M(IJK,0,M) = -ONE
85                    B_M(IJK,M) = ZERO
86     
87     ! use the average boundary cell values to compute phi (sof, Aug 23 2005)
88                    IMJK = IM_OF(IJK)
89                    IJMK = JM_OF(IJK)
90                    IJKM = KM_OF(IJK)
91                    IF (EP(WEST_OF(IJK)) > DIL_EP_S .AND. &
92                        .NOT.IS_AT_E(IMJK)) A_M(IJK,west,M) = ONE
93                    IF (EP(EAST_OF(IJK)) > DIL_EP_S .AND. &
94                        .NOT.IS_AT_E(IJK)) A_M(IJK,east,M) = ONE
95                    IF (EP(SOUTH_OF(IJK)) > DIL_EP_S .AND. &
96                        .NOT.IS_AT_N(IJMK)) A_M(IJK,south,M) = ONE
97                    IF (EP(NORTH_OF(IJK)) > DIL_EP_S .AND. &
98                        .NOT.IS_AT_N(IJK)) A_M(IJK,north,M) = ONE
99                    IF(.NOT. NO_K) THEN
100                       IF (EP(BOTTOM_OF(IJK)) > DIL_EP_S .AND. &
101                          .NOT.IS_AT_T(IJKM)) A_M(IJK,bottom,M) = ONE
102                       IF (EP(TOP_OF(IJK)) > DIL_EP_S .AND. &
103                          .NOT.IS_AT_T(IJK)) A_M(IJK,top,M) = ONE
104                    ENDIF
105     
106                    IF((A_M(IJK,west,M)+A_M(IJK,east,M)+&
107                       A_M(IJK,south,M)+A_M(IJK,north,M)+ &
108                       A_M(IJK,bottom,M)+A_M(IJK,top,M)) == ZERO) THEN
109                       B_M(IJK,M) = -PHI(IJK)
110                    ELSE
111                       A_M(IJK,0,M) = -(A_M(IJK,east,M) + A_M(IJK,west,M) +&
112                          A_M(IJK,north,M) + A_M(IJK,south,M) + &
113                          A_M(IJK,top,M)+A_M(IJK,bottom,M))
114                    ENDIF
115     
116     ! Normal case
117                 ELSE
118     
119     ! Collect the terms
120                    A_M(IJK,0,M) = -(A_M(IJK,east,M)+A_M(IJK,west,M)+&
121                                     A_M(IJK,north,M)+A_M(IJK,south,M)+&
122                                     A_M(IJK,top,M)+A_M(IJK,bottom,M)+S_P(IJK))
123     
124     ! B_m and A_m are corrected in case deferred corrections computes B_m > S_c
125     ! see CONV_DIF_PHI_DC.
126                    IF(B_M(IJK,M) < S_C(IJK) .OR. PHI(IJK) == ZERO) THEN
127                       B_M(IJK,M) = -S_C(IJK)+B_M(IJK,M)
128                    ELSE ! disable ELSE statememt if PHI can be negative
129                       A_M(IJK,0,M) = A_M(IJK,0,M) - B_M(IJK,M)/PHI(IJK)
130                       B_M(IJK,M) = -S_C(IJK)
131                    ENDIF
132     
133                 ENDIF   ! end if/else (ep_g(ijk)<=dil_ep_s)
134              ENDIF    ! end if (fluid_at(ijk))
135           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
136     
137     
138           RETURN
139           END SUBROUTINE SOURCE_PHI
140     
141     
142     
143     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
144     !                                                                      C
145     !  Subroutine: POINT_SOURCE_PHI                                        C
146     !                                                                      C
147     !  Purpose: Adds point sources to the scalar equations.                C
148     !                                                                      C
149     !  Author: J. Musser                                  Date: 10-JUN-13  C
150     !  Reviewer:                                          Date:            C
151     !                                                                      C
152     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
153           SUBROUTINE POINT_SOURCE_PHI(PHI, PS_PHI, PS_FLOW,  &
154              M, A_M, B_M)
155     
156           use compar
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           use param, only: dimension_3, dimension_m
168           use param1, only: one, small_number
169     
170           IMPLICIT NONE
171     !-----------------------------------------------
172     ! Dummy arguments
173     !-----------------------------------------------
174     
175     ! Vector b_m
176           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
177     
178     
179     ! maximum term in b_m expression
180           DOUBLE PRECISION, INTENT(IN) :: PS_PHI(DIMENSION_BC)
181     
182     ! maximum term in b_m expression
183           DOUBLE PRECISION, INTENT(IN) :: PS_FLOW(DIMENSION_BC)
184     
185           INTEGER, intent(in) :: M
186     
187     ! Septadiagonal matrix A_m
188           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
189     
190     ! Vector b_m
191           DOUBLE PRECISION, intent(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
192     
193     !-----------------------------------------------
194     ! Local Variables
195     !-----------------------------------------------
196     
197     ! Indices
198           INTEGER :: IJK, I, J, K
199           INTEGER :: PSV
200     
201     ! terms of bm expression
202           DOUBLE PRECISION pSource
203     
204     !-----------------------------------------------
205     
206           PSV_LP: do PSV = 1, DIMENSION_PS
207     
208              if(.NOT.PS_DEFINED(PSV)) cycle PSV_LP
209              if(abs(PS_FLOW(PSV)) < small_number) cycle PSV_LP
210     
211              do k = PS_K_B(PSV), PS_K_T(PSV)
212              do j = PS_J_S(PSV), PS_J_N(PSV)
213              do i = PS_I_W(PSV), PS_I_E(PSV)
214     
215                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
216                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
217     
218                 ijk = funijk(i,j,k)
219                 if(.NOT.fluid_at(ijk)) cycle
220     
221                 if(A_M(IJK,0,M) == -ONE .AND. B_M(IJK,M) == -PHI(IJK)) then
222                    B_M(IJK,M) = -PS_PHI(PSV)
223                 else
224     
225     ! Calculate the mass flow rate for this cell. (mass/time)
226                    pSource = PS_FLOW(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
227     
228                    A_M(IJK,0,M) = A_M(IJK,0,M) - pSource
229                    B_M(IJK,M) = B_M(IJK,M) - PS_PHI(PSV) * pSource
230     
231                 endif
232     
233              enddo
234              enddo
235              enddo
236     
237           enddo PSV_LP
238     
239           RETURN
240           END SUBROUTINE POINT_SOURCE_PHI
241