File: /nfs/home/0/users/jenkins/mfix.git/model/source_phi.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 SUBROUTINE SOURCE_PHI(S_P, S_C, EP, PHI, M, A_M, B_M)
26
27
28
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
51
52
53 DOUBLE PRECISION, INTENT(IN) :: S_p(DIMENSION_3)
54
55 DOUBLE PRECISION, INTENT(IN) :: S_C(DIMENSION_3)
56
57 DOUBLE PRECISION, INTENT(IN) :: EP(DIMENSION_3)
58
59 DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
60
61 INTEGER, INTENT(IN) :: M
62
63 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
64
65 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
66
67
68
69
70 INTEGER :: IJK, IMJK, IJMK, IJKM
71
72
73 DO IJK = ijkstart3, ijkend3
74
75 IF (FLUID_AT(IJK)) THEN
76
77
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
89 = 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
118 ELSE
119
120
121 (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
126
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
130 (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
135 ENDIF
136 ENDDO
137
138
139 RETURN
140 END SUBROUTINE SOURCE_PHI
141
142
143
144
145
146
147
148
149
150
151
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
164 use bc
165 use usr
166 use functions
167
168 IMPLICIT NONE
169
170
171
172
173
174 DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
175
176
177
178 DOUBLE PRECISION, INTENT(IN) :: PS_PHI(DIMENSION_BC)
179
180
181 DOUBLE PRECISION, INTENT(IN) :: PS_FLOW(DIMENSION_BC)
182
183 INTEGER, intent(in) :: M
184
185
186 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
187
188
189 DOUBLE PRECISION, intent(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
190
191
192
193
194
195
196 INTEGER :: IJK, I, J, K
197 INTEGER :: PSV
198
199
200 DOUBLE PRECISION pSource
201
202
203
204
205
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
219
220 = 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
228 = 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