File: /nfs/home/0/users/jenkins/mfix.git/model/set_ps.f
1
2
3
4
5
6
7
8
9
10
11 SUBROUTINE SET_PS
12
13 use param
14 use run
15 use physprop
16 use ps
17 use compar
18 use geometry
19 use mpi_utility
20 use functions
21
22 implicit none
23
24 INTEGER :: IJK, I, J, K, M, N
25
26 INTEGER PSV
27
28 CHARACTER(LEN=64) :: eMsg
29
30 INTEGER :: PS_SIZE
31
32 DOUBLE PRECISION, allocatable :: lData_dp(:)
33 DOUBLE PRECISION, allocatable :: gData_dp(:)
34
35 LOGICAL, parameter :: dbg_PS = .FALSE.
36
37 if(.NOT.POINT_SOURCE) return
38
39
40 L50: do PSV = 1, DIMENSION_PS
41
42 IF(.NOT.PS_DEFINED(PSV)) cycle L50
43
44
45
46 CALL CALC_PS_VEL_MAG(PS_VEL_MAG_g(PSV), PS_U_g(PSV), &
47 PS_V_g(PSV), PS_W_g(PSV))
48
49 CALL CALC_PS_CpxMFLOW(PS_CpxMFLOW_g(PSV), PS_MASSFLOW_g(PSV), &
50 PS_T_g(PSV), PS_X_g(PSV,:), 0, C_PG0, DIM_N_g, MW_g)
51
52 do M=1, MMAX
53 CALL CALC_PS_VEL_MAG(PS_VEL_MAG_s(PSV,M), PS_U_s(PSV,M), &
54 PS_V_s(PSV,M), PS_W_s(PSV,M))
55
56 CALL CALC_PS_CpxMFLOW(PS_CpxMFLOW_s(PSV,M), &
57 PS_MASSFLOW_s(PSV,M), PS_T_s(PSV,M), PS_X_s(PSV,M,:), M,&
58 C_PS0(M), DIM_N_s, MW_s(M,:))
59 enddo
60
61
62
63
64
65 = (PS_I_E(PSV) - PS_I_W(PSV) + 1) * &
66 (PS_J_N(PSV) - PS_J_S(PSV) + 1)
67 if(DO_K) PS_SIZE = PS_SIZE * (PS_K_T(PSV) - PS_K_B(PSV) + 1)
68
69 if(PS_SIZE < 1) then
70 eMsg = ''; write(eMsg,"('Invalid PS size: ', I4)")PS_SIZE
71 goto 500
72 endif
73
74
75 allocate(lData_dp( 0:numPEs-1)); lData_dp = ZERO
76 allocate(gData_dp( 0:numPEs-1)); gData_dp = ZERO
77
78
79
80
81
82 (PSV) = ZERO
83
84 do k = PS_K_B(PSV), PS_K_T(PSV)
85 do j = PS_J_S(PSV), PS_J_N(PSV)
86 do i = PS_I_W(PSV), PS_I_E(PSV)
87
88 if(IS_ON_myPE_owns(I, J, K)) then
89 IF (DEAD_CELL_AT(I,J,K)) CYCLE
90 = funijk(i,j,k)
91
92 if(fluid_at(ijk)) &
93 lData_dp(myPE) = lData_dp(myPE) + VOL(IJK)
94 endif
95
96 enddo
97 enddo
98 enddo
99
100
101
102
103 CALL global_all_sum(lData_dp, gData_dp)
104
105 PS_VOLUME(PSV) = sum(gData_dp)
106 if(PS_VOLUME(PSV) == ZERO) then
107 eMsg = 'No PS_VOLUME == ZERO'
108 CALL DEBUG_PS(PSV, PS_SIZE)
109 goto 501
110 endif
111
112 if(allocated(lData_dp)) deallocate(lData_dp)
113 if(allocated(gData_dp)) deallocate(gData_dp)
114
115
116 IF(dbg_PS) CALL DEBUG_PS(PSV, PS_SIZE)
117
118
119 enddo L50
120
121 return
122
123 500 continue
124 if(myPE == PE_IO) &
125 write(*,"('PointSource setup Error: ',A)") trim(eMsg)
126
127 call mfix_exit(myPE)
128
129
130 501 continue
131 write(*,"('PointSource setup Error: ',I3,2x,A)") myPE, trim(eMsg)
132
133 call mfix_exit(myPE)
134
135
136 RETURN
137
138 contains
139
140
141
142
143
144
145
146
147
148
149
150
151 SUBROUTINE CALC_PS_CpxMFLOW(CpxMFLOW, PS_MFLOW, PS_T, PS_X, lM, &
152 Cp0, lDIM_N, lMW)
153
154 use constant, only: GAS_CONST_cal
155 use read_thermochemical, only: calc_CpoR
156
157 DOUBLE PRECISION, intent(out) :: CpxMFLOW
158
159 INTEGER, intent(in) :: lM
160 INTEGER, intent(in) :: lDIM_N
161
162 DOUBLE PRECISION, intent(in) :: Cp0
163 DOUBLE PRECISION, intent(in) :: PS_MFLOW
164 DOUBLE PRECISION, intent(in) :: PS_T
165 DOUBLE PRECISION, intent(in) :: PS_X(lDIM_N)
166
167 DOUBLE PRECISION, intent(in) :: lMW(lDIM_N)
168
169 INTEGER :: IER
170
171
172
173 if(.NOT.ENERGY_EQ .OR. PS_MFLOW == ZERO) then
174 CpxMFLOW = ZERO
175 return
176 endif
177
178
179 if(Cp0 == UNDEFINED) then
180 IF(.NOT.database_read) call read_database0(IER)
181 CpxMFLOW = ZERO
182 do N = 1, NMAX(lM)
183 CpxMFLOW = CpxMFLOW + PS_X(N) * (GAS_CONST_cal / lMW(N)) * &
184 calc_CpoR(PS_T, lM, N, IER)
185 enddo
186 else
187 CpxMFLOW = Cp0
188 endif
189
190
191 if(UNITS == 'SI') CpxMFLOW = 4.183925d3*CpxMFLOW
192 CpxMFLOW = CpxMFLOW*PS_MFLOW
193
194 END SUBROUTINE CALC_PS_CpxMFLOW
195
196
197
198
199
200
201
202
203
204
205
206 SUBROUTINE CALC_PS_VEL_MAG(VEL_MAG, lU, lV, lW)
207
208 DOUBLE PRECISION, intent(inout) :: VEL_MAG
209 DOUBLE PRECISION, intent(inout) :: lU, lV, lW
210
211
212 = lU**2 + lV**2
213 if(DO_K) VEL_MAG = VEL_MAG + lW**2
214
215 VEL_MAG = sqrt(VEL_MAG)
216
217 if(VEL_MAG > small_number) then
218 lU = lU/VEL_MAG
219 lV = lV/VEL_MAG
220 lW = lW/VEL_MAG
221 else
222 VEL_MAG = ZERO
223 lU = ZERO
224 lV = ZERO
225 lW = ZERO
226 endif
227
228
229 END SUBROUTINE CALC_PS_VEL_MAG
230
231
232 END SUBROUTINE SET_PS
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248 SUBROUTINE DEBUG_PS(lPSV, lPS_SIZE)
249
250 use bc
251 use compar
252 use constant
253 use cutcell
254 use fldvar
255 use geometry
256 use ic
257 use indices
258 use mflux
259 use parallel
260 use param
261 use param1
262 use physprop
263 use run
264 use sendrecv
265 use toleranc
266 use usr
267 use rxns
268 use ps
269 use mpi_utility
270 use functions
271
272 IMPLICIT NONE
273
274
275 INTEGER, intent(in) :: lPSV
276
277 INTEGER, intent(in) :: lPS_SIZE
278
279 INTEGER :: IJK, I, J, K, M, N
280
281 INTEGER :: lc1
282
283 INTEGER, allocatable :: lFlags_i(:,:)
284 INTEGER, allocatable :: gFlags_i(:,:)
285
286 if(myPE == PE_IO) then
287 write(*,"(3/,3x,'Debug Point Source Index: ',I3)") lPSV
288 write(*,"(/3x,'Size: ',I4)") lPS_SIZE
289 endif
290
291 allocate(lFlags_i(lPS_SIZE,1:2) ); lFlags_i = 0
292 allocate(gFlags_i(lPS_SIZE,1:2) ); gFlags_i = 0
293
294 lc1 = 0
295
296 do k = PS_K_B(lPSV), PS_K_T(lPSV)
297 do j = PS_J_S(lPSV), PS_J_N(lPSV)
298 do i = PS_I_W(lPSV), PS_I_E(lPSV)
299
300 lc1 = lc1 + 1
301 if(IS_ON_myPE_owns(I, J, K)) then
302 IF (DEAD_CELL_AT(I,J,K)) CYCLE
303 = funijk(i,j,k)
304 if(fluid_at(ijk)) then
305 lFlags_i(lc1,1) = myPE
306 lFlags_i(lc1,2) = FLAG(IJK)
307 endif
308 endif
309 enddo
310 enddo
311 enddo
312
313
314 CALL global_sum(lFlags_i, gFlags_i)
315
316
317 if(myPE == PE_IO) then
318 write(*,"(/5x,'Location:')")
319 write(*,"( 5x,'X:',2(2x,g12.5),' :: ',2(2x,I4))")&
320 PS_X_w(lPSV), PS_X_e(lPSV), PS_I_w(lPSV), PS_I_e(lPSV)
321 write(*,"( 5x,'Y:',2(2x,g12.5),' :: ',2(2x,I4))")&
322 PS_Y_s(lPSV), PS_Y_n(lPSV), PS_J_s(lPSV), PS_J_n(lPSV)
323 if(DO_K)write(*,"( 5x,'Z:',2(2x,g12.5),' :: ',2(2x,I4))")&
324 PS_Z_b(lPSV), PS_Z_t(lPSV), PS_K_b(lPSV), PS_K_t(lPSV)
325
326 write(*,"(/5x,'Volume: ',g12.5)") PS_VOLUME(lPSV)
327
328
329 if(PS_MASSFLOW_G(lPSV) > small_number) then
330 write(*,"(//5x,'Point Source Gas Phase:')")
331 write(*,"(7x,'Mass Flow Rate: ',g12.5)")PS_MASSFLOW_G(lPSV)
332 write(*,"(7x,'Velocity Magnitude: ',g12.5)") PS_VEL_MAG_g(lPSV)
333 write(*,"(7x,'Normal:')")
334 write(*,"(9x,'x-Axis: ',g12.5)")PS_U_g(lPSV)
335 write(*,"(9x,'y-Axis: ',g12.5)")PS_V_g(lPSV)
336 if(DO_K) write(*,"(9x,'z-Axis: ',g12.5)")PS_W_g(lPSV)
337 if(energy_eq) &
338 write(*,"(7x,'Temperature: ',g12.5)")PS_T_g(lPSV)
339 if(species_eq(0)) then
340 write(*,"(7x,'Species Composition:')")
341 do n=1,nmax(0)
342 write(*,"(9x,A,': ',g12.5)") &
343 trim(SPECIES_ALIAS_g(n)), PS_X_g(lPSV,N)
344 enddo
345 endif
346 else
347 write(*,"(//5x,'No gas phase point source.')")
348 endif
349
350
351 do m=1,mmax
352 if(PS_MASSFLOW_S(lPSV,M) > small_number) then
353 write(*,"(//5x,'Point Source Solids Phase ',I1,':')")M
354 write(*,"(7x,'Mass Flow Rate: ',g12.5)")PS_MASSFLOW_S(lPSV,M)
355 write(*,"(7x,'Velocity Magnitude: ',g12.5)") PS_VEL_MAG_s(lPSV,M)
356 write(*,"(7x,'Normal:')")
357 write(*,"(9x,'x-Axis: ',g12.5)")PS_U_s(lPSV,M)
358 write(*,"(9x,'y-Axis: ',g12.5)")PS_V_s(lPSV,M)
359 if(DO_K) write(*,"(9x,'z-Axis: ',g12.5)")PS_W_s(lPSV,M)
360 if(energy_eq) &
361 write(*,"(7x,'Temperature: ',g12.5)")PS_T_s(lPSV,M)
362 if(species_eq(m)) then
363 write(*,"(7x,'Species Composition:')")
364 do n=1,nmax(m)
365 write(*,"(9x,A,': ',g12.5)") &
366 trim(SPECIES_ALIAS_s(m,n)), PS_X_s(lPSV,m,N)
367 enddo
368 endif
369 else
370 write(*,"(//5x,'No solids phase ',I1,' point source.')") m
371 endif
372 enddo
373
374 write(*,"(//5x,'Point Source Cells:')")
375 write(*,"(9x,'IJK',3(6x,A1),3x,'OWNS',3x,'FLAG')") 'I','J','K'
376
377 lc1 = 0
378 do k = PS_K_B(lPSV), PS_K_T(lPSV)
379 do j = PS_J_S(lPSV), PS_J_N(lPSV)
380 do i = PS_I_W(lPSV), PS_I_E(lPSV)
381 IF (DEAD_CELL_AT(I,J,K)) CYCLE
382 = lc1 + 1
383 write(*,"(4x,I8,5(3x,I4))") IJK, I, J, K, gFlags_i(lc1,:)
384 enddo
385 enddo
386 enddo
387
388 endif
389
390 if(allocated(lFlags_i)) deallocate(lFlags_i)
391 if(allocated(gFlags_i)) deallocate(gFlags_i)
392
393 END SUBROUTINE DEBUG_PS
394
395