File: /nfs/home/0/users/jenkins/mfix.git/model/physical_prop.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
26
27
28
29 SUBROUTINE PHYSICAL_PROP(IER, LEVEL)
30
31 use compar
32 use funits
33 use geometry
34 use indices
35 use mpi_utility
36 use param1
37 use physprop
38
39 use coeff, only: DENSITY
40 use coeff, only: SP_HEAT
41 use coeff, only: PSIZE
42
43 implicit none
44
45
46
47
48 INTEGER, intent(inout) :: IER
49 INTEGER, intent(in) :: LEVEL
50
51
52
53
54
55
56
57
58 INTEGER :: Err_l(0:numPEs-1)
59 INTEGER :: Err_g(0:numPEs-1)
60
61
62 = 0
63
64
65
66 if(LEVEL == 0) then
67 if(DENSITY(0)) CALL PHYSICAL_PROP_ROg
68 if(any(DENSITY(1:DIMENSION_M))) CALL PHYSICAL_PROP_ROs
69
70
71
72 elseif(LEVEL == 1) then
73 if(SP_HEAT(0)) CALL PHYSICAL_PROP_CPg
74
75 if(any(SP_HEAT(1:DIMENSION_M))) CALL PHYSICAL_PROP_CPs
76 if(any(PSIZE(1:DIMENSION_M))) CALL PHYSICAL_PROP_Dp
77
78
79
80
81
82 elseif(LEVEL == 2) then
83 if(DENSITY(0)) CALL PHYSICAL_PROP_ROg
84 if(SP_HEAT(0)) CALL PHYSICAL_PROP_CPg
85 if(any(DENSITY(1:DIMENSION_M))) CALL PHYSICAL_PROP_ROs
86 if(any(SP_HEAT(1:DIMENSION_M))) CALL PHYSICAL_PROP_CPs
87 if(any(PSIZE(1:DIMENSION_M))) CALL PHYSICAL_PROP_Dp
88 endif
89
90
91
92
93 CALL global_all_sum(Err_l, Err_g)
94 IER = maxval(Err_g)
95 if(IER == 0) return
96
97
98
99
100
101
102 IF(IER == 901 .OR. IER == 902) then
103 if(myPE == PE_IO) then
104 write(*,2000) IER
105 write(UNIT_LOG,2000) IER
106 endif
107 CALL MFIX_EXIT(myPE)
108 ENDIF
109
110 return
111
112 2000 FORMAT(/1X,70('*')/' From: PHYSICAL_PROP',/' Fatal Error 2000:', &
113 ' calc_CpoR reporetd an invalid temperature: 0x0', I3/, &
114 'See Cp.log for details. Calling MFIX_EXIT.',/1X,70('*')/)
115
116 contains
117
118
119
120
121
122
123
124
125
126
127 SUBROUTINE PHYSICAL_PROP_ROg
128
129
130
131
132
133 use fldvar, only: X_g
134
135 use fldvar, only: T_g
136
137 use fldvar, only: RO_g
138
139 use fldvar, only: P_g
140
141 use fldvar, only: EP_g
142
143 use fldvar, only: ROP_g
144
145 use toleranc, only: OMW_MAX
146
147 use run, only: REPORT_NEG_DENSITY
148 use functions
149
150 use eos, only: EOSG
151
152 implicit none
153
154
155
156
157 DOUBLE PRECISION :: MW
158
159
160 INTEGER :: IJK
161
162
163 LOGICAL wHeader
164
165
166 = .TRUE.
167
168
169 IF(.NOT.database_read) call read_database0(IER)
170
171 IJK_LP: DO IJK = IJKSTART3, IJKEND3
172 IF(WALL_AT(IJK)) cycle IJK_LP
173 IF (MW_AVG == UNDEFINED) THEN
174
175 = SUM(X_G(IJK,:NMAX(0))/MW_G(:NMAX(0)))
176 MW = ONE/MAX(MW,OMW_MAX)
177 MW_MIX_G(IJK) = MW
178
179 (IJK) = EOSG(MW,P_G(IJK),T_G(IJK))
180 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
181 ELSE
182 RO_G(IJK) = EOSG(MW_AVG,P_G(IJK),T_G(IJK))
183 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
184 ENDIF
185
186 IF(RO_G(IJK) < ZERO) THEN
187 Err_l(myPE) = 100
188 IF(REPORT_NEG_DENSITY)CALL ROgErr_LOG(IJK, wHeader)
189 ENDIF
190 ENDDO IJK_LP
191
192
193 RETURN
194 END SUBROUTINE PHYSICAL_PROP_ROg
195
196
197
198
199
200
201
202
203
204
205
206 SUBROUTINE PHYSICAL_PROP_ROs
207
208
209
210
211 use fldvar, only: X_s
212
213 use fldvar, only: ROP_s, RO_s
214
215 use physprop, only: BASE_ROs
216
217 use physprop, only: X_S0
218
219 use physprop, only: INERT_SPECIES
220
221 use run, only: SOLVE_ROs
222
223 use run, only: REPORT_NEG_DENSITY
224
225 use toleranc, only: DIL_EP_s
226 use functions
227
228 use eos, only: EOSS
229
230 implicit none
231
232
233
234
235 INTEGER :: IJK
236 INTEGER :: M
237
238 INTEGER :: IIS
239
240 LOGICAL wHeader
241
242 DOUBLE PRECISION :: minROPs
243
244 M_LP: DO M=1, MMAX
245 IF(.NOT.SOLVE_ROs(M)) cycle M_LP
246
247 = .TRUE.
248
249 = INERT_SPECIES(M)
250
251 = BASE_ROs(M)*(10.0d0*DIL_EP_s)
252
253
254 IJK_LP: DO IJK = IJKSTART3, IJKEND3
255 IF(WALL_AT(IJK)) cycle IJK_LP
256 IF(ROP_s(IJK,M) > minROPs) THEN
257 RO_S(IJK,M) = EOSS(BASE_ROs(M), X_s0(M,IIS), &
258 X_s(IJK,M,IIS))
259 ELSE
260 RO_s(IJK,M) = BASE_ROs(M)
261 ENDIF
262
263
264 IF(RO_S(IJK,M) <= ZERO) THEN
265 Err_l(myPE) = 101
266 IF(REPORT_NEG_DENSITY) CALL ROsErr_LOG(IJK, M, wHeader)
267 ENDIF
268 ENDDO IJK_LP
269 ENDDO M_LP
270
271
272 RETURN
273 END SUBROUTINE PHYSICAL_PROP_ROs
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288 SUBROUTINE PHYSICAL_PROP_CPg
289
290
291
292
293 use constant, only: RGAS => GAS_CONST_cal
294
295 use toleranc, only: OMW_MAX
296
297 use fldvar, only: X_g
298
299 use fldvar, only: T_g
300
301 use run, only: UNITS
302 use functions
303
304 use read_thermochemical, only: calc_CpoR
305
306 implicit none
307
308
309
310
311 DOUBLE PRECISION :: lCp
312
313 INTEGER :: lCP_Err
314 INTEGER :: gCP_Err
315
316 INTEGER :: IJK
317 INTEGER :: N
318
319
320
321
322
323 IF(.NOT.database_read) CALL read_database0(IER)
324
325 gCP_Err = 0
326 lCP_Err = 0
327
328 IJK_LP: DO IJK = IJKSTART3, IJKEND3
329 IF(WALL_AT(IJK)) CYCLE IJK_LP
330
331 (IJK) = ZERO
332 DO N = 1, NMAX(0)
333 lCp = calc_CpoR(T_G(IJK), 0, N, lCP_Err)
334 C_PG(IJK) = C_PG(IJK) + X_g(IJK,N) * lCp * RGAS / MW_g(N)
335 gCP_Err = max(gCP_Err, lCP_Err)
336 ENDDO
337 ENDDO IJK_LP
338
339
340
341 IF(UNITS == 'SI') C_PG = 4.183925d3 * C_PG
342
343
344 IF(gCP_Err /= 0) Err_l(myPE) = 800 + gCP_Err
345
346 RETURN
347 END SUBROUTINE PHYSICAL_PROP_CPg
348
349
350
351
352
353
354
355
356
357
358
359
360
361 SUBROUTINE PHYSICAL_PROP_CPs
362
363
364 use constant, only: RGAS => GAS_CONST_cal
365
366 use run, only: UNITS
367 use toleranc, only: OMW_MAX
368 use fldvar, only: T_s
369 use fldvar, only: X_s
370 use functions
371 use read_thermochemical, only: calc_CpoR
372
373 implicit none
374
375 DOUBLE PRECISION :: lCp
376
377
378 INTEGER :: IJK
379 INTEGER :: M
380 INTEGER :: N
381
382
383 INTEGER :: lCP_Err
384
385
386
387 IF(.NOT.database_read) CALL read_database0(IER)
388
389 lCP_Err = 0
390
391 M_LP: DO M=1, MMAX
392 IJK_LP: DO IJK = IJKSTART3, IJKEND3
393 IF(WALL_AT(IJK)) CYCLE IJK_LP
394
395 (IJK, M) = ZERO
396
397 DO N = 1, NMAX(M)
398 lCp = calc_CpoR(T_S(IJK,M), M, N, lCP_Err)
399 C_PS(IJK,M) = C_PS(IJK,M) + X_s(IJK,M,N) * &
400 (lCp * RGAS / MW_s(M,N))
401 ENDDO
402
403 ENDDO IJK_LP
404 ENDDO M_LP
405
406
407
408 IF(UNITS == 'SI') C_PS = 4.183925d3 * C_PS
409
410 END SUBROUTINE PHYSICAL_PROP_CPs
411
412
413
414
415
416
417
418
419
420
421 SUBROUTINE PHYSICAL_PROP_Dp
422
423 use run, only: CALL_DQMOM
424
425 use fldvar, only: scalar
426 use scalars, only: phase4scalar
427
428 use fldvar, only: ROP_s
429 use fldvar, only: RO_s
430 use fldvar, only: D_p, EP_S
431 use functions
432
433 implicit none
434
435
436 INTEGER :: IJK
437 INTEGER :: M
438
439
440 INTEGER :: lM
441
442 IF(.NOT.CALL_DQMOM) return
443
444 M_LP: DO M=1, MMAX
445
446 lM = phase4scalar(M)
447
448 IF(.NOT.PSIZE(M)) CYCLE M_LP
449 IJK_LP: DO IJK = IJKSTART3, IJKEND3
450 IF(WALL_AT(IJK)) CYCLE IJK_LP
451
452 IF(EP_s(IJK,lM) > small_number) D_p(IJK,M)= Scalar(IJK,lM)
453
454 ENDDO IJK_LP
455 ENDDO M_LP
456
457 return
458 END SUBROUTINE PHYSICAL_PROP_Dp
459
460
461
462
463
464
465
466
467
468
469
470
471
472 SUBROUTINE ROgErr_LOG(IJK, tHeader)
473
474
475 use run, only: TIME
476
477 use fldvar, only: T_g
478
479 use fldvar, only: RO_g
480
481 use fldvar, only: P_g
482 use cutcell
483
484 INTEGER, intent(in) :: IJK
485 LOGICAL, intent(inout) :: tHeader
486
487 LOGICAL :: lExists
488 CHARACTER(LEN=32) :: lFile
489 INTEGER, parameter :: lUnit = 4868
490 LOGICAL, save :: fHeader = .TRUE.
491
492
493 lFile = '';
494 if(numPEs > 1) then
495 write(lFile,"('ROgErr_',I4.4,'.log')") myPE
496 else
497 write(lFile,"('ROgErr.log')")
498 endif
499 inquire(file=trim(lFile),exist=lExists)
500 if(lExists) then
501 open(lUnit,file=trim(adjustl(lFile)), &
502 status='old', position='append')
503 else
504 open(lUnit,file=trim(adjustl(lFile)), status='new')
505 endif
506
507 if(fHeader) then
508 write(lUnit,1000)
509 fHeader = .FALSE.
510 endif
511
512 if(tHeader) then
513 write(lUnit,"(/2x,'Simulation time: ',g12.5)") TIME
514 tHeader = .FALSE.
515 endif
516
517 write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
518 write(lUnit,"(6x,A,1X,g12.5)",ADVANCE='NO') 'RO_g:', RO_g(IJK)
519 write(lUnit,"(2x,A,1X,g12.5)",ADVANCE='NO') 'P_g:', P_g(IJK)
520 write(lUnit,"(2x,A,1X,g12.5)") 'T_g:', T_g(IJK)
521 if(CARTESIAN_GRID) then
522 write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
523 write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
524 write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
525 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
526 endif
527
528 close(lUnit)
529
530 RETURN
531
532 1000 FORMAT(2X,'One or more cells have reported a negative gas', &
533 ' density (RO_g(IJK)). If',/2x,'this is a persistent issue,', &
534 ' lower UR_FAC(1) in mfix.dat.')
535
536 1001 FORMAT(/4X,'IJK: ',I8,7X,'I: ',I4,' J: ',I4,' K: ',I4)
537
538 END SUBROUTINE ROgErr_LOG
539
540
541
542
543
544
545
546
547
548
549
550
551 SUBROUTINE ROsErr_LOG(IJK, M, tHeader)
552
553
554
555
556 use run, only: TIME
557
558 use fldvar, only: X_s
559
560 use fldvar, only: RO_s
561
562 use physprop, only: BASE_ROs
563
564 use physprop, only: X_S0
565
566 use physprop, only: INERT_SPECIES
567
568
569 use cutcell
570
571
572
573
574
575 INTEGER, intent(in) :: IJK, M
576
577 LOGICAL, intent(inout) :: tHeader
578
579 INTEGER :: N
580
581 LOGICAL :: lExists
582 CHARACTER(LEN=32) :: lFile
583 INTEGER, parameter :: lUnit = 4868
584 LOGICAL, save :: fHeader = .TRUE.
585
586
587 lFile = '';
588 if(numPEs > 1) then
589 write(lFile,"('ROsErr_',I4.4,'.log')") myPE
590 else
591 write(lFile,"('ROsErr.log')")
592 endif
593 inquire(file=trim(lFile),exist=lExists)
594 if(lExists) then
595 open(lUnit,file=trim(adjustl(lFile)), &
596 status='old', position='append')
597 else
598 open(lUnit,file=trim(adjustl(lFile)), status='new')
599 endif
600
601 if(fHeader) then
602 write(lUnit,1000)
603 fHeader = .FALSE.
604 endif
605
606 if(tHeader) then
607 write(lUnit,"(/2x,'Simulation time: ',g12.5,' Phase: ',I2)")&
608 TIME, M
609 tHeader = .FALSE.
610 endif
611
612 N = INERT_SPECIES(M)
613 write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
614 write(lUnit,"(6x,A,1X,g12.5)",advance='no') 'RO_s:', RO_s(IJK,M)
615 write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'Base:', BASE_ROs(M)
616 write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'X_s0:', X_s0(M,N)
617 write(lUnit,"(2x,A,1X,g12.5)") 'X_s:', X_s(IJK,M,N)
618
619 if(CARTESIAN_GRID) then
620 write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
621 write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
622 write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
623 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
624 endif
625
626 close(lUnit)
627
628 RETURN
629
630 1000 FORMAT(2X,'One or more cells have reported a negative gas', &
631 ' density (RO_g(IJK)). If',/2x,'this is a persistent issue,', &
632 ' lower UR_FAC(1) in mfix.dat.')
633
634 1001 FORMAT( 4X,'IJK: ',I8,7X,'I: ',I4,' J: ',I4,' K: ',I4)
635
636 END SUBROUTINE ROsErr_LOG
637
638 END SUBROUTINE PHYSICAL_PROP
639