File: /nfs/home/0/users/jenkins/mfix.git/model/iterate.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
30
31
32 SUBROUTINE ITERATE(IER, NIT)
33
34
35
36
37 USE param
38 USE param1
39 USE toleranc
40 USE run
41 USE physprop
42 USE geometry
43 USE fldvar
44 USE output
45 USE indices
46 USE funits
47 USE time_cpu
48 USE pscor
49 USE leqsol
50 USE visc_g
51 USE pgcor
52 USE cont
53 USE scalars
54 USE compar
55 USE mpi_utility
56 USE discretelement
57 USE residual
58 USE cutcell
59 USE vtk
60 USE dashboard
61 USE qmom_kinetic_equation
62 USE stiff_chem, only : STIFF_CHEMISTRY
63 USE rxns, only : USE_RRATES, NO_OF_RXNS
64 USE mms, only: USE_MMS
65 IMPLICIT NONE
66
67
68
69
70 INTEGER, INTENT(INOUT) :: IER
71
72 INTEGER, INTENT(OUT) :: NIT
73
74
75
76
77 DOUBLE PRECISION :: CPU_NOW
78
79 DOUBLE PRECISION :: TLEFT
80
81
82 INTEGER :: MUSTIT
83
84 DOUBLE PRECISION :: NORMg, NORMs
85
86 LOGICAL :: SETg, SETs
87
88 DOUBLE PRECISION :: RESg, RESs
89
90 DOUBLE PRECISION :: SMASS
91
92 DOUBLE PRECISION :: HLOSS
93
94 INTEGER :: M
95
96 DOUBLE PRECISION :: Vavg
97 DOUBLE PRECISION :: errorpercent(0:MMAX)
98 LOGICAL :: ABORT_IER
99 CHARACTER(LEN=4) :: TUNIT
100
101
102 INTEGER :: lErrMsg
103
104 CHARACTER(LEN=32) :: lMsg
105
106
107
108
109 DOUBLE PRECISION, EXTERNAL :: VAVG_U_G, VAVG_V_G, VAVG_W_G, &
110 VAVG_U_S, VAVG_V_S, VAVG_W_S
111 DOUBLE PRECISION :: WALL_TIME
112
113
114
115
116
117
118
119 = DT
120 NIT = 0
121 MUSTIT = 0
122 RESG = ZERO
123 RESS = ZERO
124
125 IF (NORM_G == UNDEFINED) THEN
126 NORMG = ONE
127 SETG = .FALSE.
128 ELSE
129 NORMG = NORM_G
130 SETG = .TRUE.
131 ENDIF
132
133 IF (NORM_S == UNDEFINED) THEN
134 NORMS = ONE
135 SETS = .FALSE.
136 ELSE
137 NORMS = NORM_S
138 SETS = .TRUE.
139 ENDIF
140
141 LEQ_ADJUST = .FALSE.
142
143
144
145 CALL INIT_RESID (IER)
146
147
148
149 IF(CYCLIC) CALL GoalSeekMassFlux(NIT, MUSTIT, .false.)
150
151
152 IF (FULL_LOG) THEN
153 TLEFT = (TSTOP - TIME)*CPUOS
154 CALL GET_TUNIT (TLEFT, TUNIT)
155
156 IF (DT == UNDEFINED) THEN
157 CALL GET_SMASS (SMASS)
158 IF(myPE.eq.PE_IO) THEN
159 WRITE (*, '(/A,G12.5, A,F9.3,1X,A)') &
160 ' Starting solids mass = ', SMASS
161 ENDIF
162 ELSE
163 IF(myPE.eq.PE_IO) THEN
164 IF ((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
165 (CN_ON.AND.RUN_TYPE /= 'NEW' .AND.&
166 NSTEP >= (NSTEPRST+1))) THEN
167 WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)')&
168 ' Time = ', TIME, ' Dt = ', 2.D0*DT
169 ELSE
170 WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)') &
171 ' Time = ', TIME, ' Dt = ', DT
172 ENDIF
173 ENDIF
174 ENDIF
175 ENDIF
176
177 CALL CALC_RESID_MB(0, errorpercent)
178
179
180
181 CALL CONV_ROP(IER)
182 CALL CALC_MFLUX (IER)
183 CALL SET_BC1
184
185
186 IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
187
188
189 = 'Run diverged/stalled'
190
191
192
193 CONTINUE
194 MUSTIT = 0
195 NIT = NIT + 1
196 PHIP_OUT_ITER=NIT
197
198
199
200 IF (.NOT.SETG) THEN
201 IF (RESG > SMALL_NUMBER) THEN
202 NORMG = RESG
203 SETG = .TRUE.
204 ENDIF
205 ENDIF
206 IF (.NOT.SETS) THEN
207 IF (RESS > SMALL_NUMBER) THEN
208 NORMS = RESS
209 SETS = .TRUE.
210 ENDIF
211 ENDIF
212
213
214
215 IF (CALL_USR) CALL USR2
216
217
218 CALL CALC_COEFF(IER, 1)
219 IF (IER_MANAGER(IER)) goto 1000
220
221
222 IF(NScalar /= 0) CALL SCALAR_PROP(IER)
223
224
225 IF(K_Epsilon) CALL K_Epsilon_PROP(IER)
226
227
228
229 IF(USE_MMS) CALL CALC_TRD_AND_TAU(IER)
230
231
232 CALL SOLVE_VEL_STAR(IER)
233
234
235 IF(CYLINDRICAL) CALL RADIAL_VEL_CORRECTION
236
237
238 CALL PHYSICAL_PROP(IER, 0)
239 IF (IER_MANAGER(IER)) goto 1000
240
241
242 CALL CALC_RRATE(IER)
243
244
245
246 IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
247 DES_CONTINUUM_HYBRID) THEN
248 IF (MMAX > 0) THEN
249
250 IF(USE_MMS) THEN
251 CALL SOLVE_CONTINUITY(0,IER)
252
253 ELSE
254 IF(MMAX == 1 .AND. MCP /= UNDEFINED_I)THEN
255
256
257 CALL CALC_K_CP (K_CP, IER)
258 CALL SOLVE_EPP (NORMS, RESS, IER)
259 CALL CORRECT_1 (IER)
260 ELSE
261
262
263
264
265 DO M=1,SMAX
266
267
268
269
270
271
272
273 CALL SOLVE_CONTINUITY(M,IER)
274
275 ENDDO
276 ENDIF
277 ENDIF
278
279 IF(KT_TYPE_ENUM == GHD_2007) CALL ADJUST_EPS_GHD
280
281 CALL CALC_VOL_FR (P_STAR, RO_G, ROP_G, EP_G, ROP_S, IER)
282 IF (IER_MANAGER(IER)) goto 1000
283
284 ENDIF
285
286 ENDIF
287
288
289
290
291 IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
292 DES_CONTINUUM_HYBRID) THEN
293 IF (MMAX > 0 .AND. .NOT.FRICTION) &
294 CALL CALC_P_STAR (EP_G, P_STAR, IER)
295 ENDIF
296
297
298 CALL CONV_ROP(IER)
299
300 IF (RO_G0 /= ZERO) THEN
301
302 CALL SOLVE_PP_G (NORMG, RESG, IER)
303
304 CALL CORRECT_0 (IER)
305 ENDIF
306
307
308 CALL PHYSICAL_PROP(IER, 0)
309 IF (IER_MANAGER(IER)) goto 1000
310
311
312
313
314
315 IF(.NOT. K_EPSILON) CALL SET_WALL_BC (IER)
316
317
318 CALL CALC_MFLUX (IER)
319 CALL SET_BC1
320
321
322 IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
323
324
325 IF (ENERGY_EQ) THEN
326 CALL SOLVE_ENERGY_EQ (IER)
327 IF (IER_MANAGER(IER)) goto 1000
328 ENDIF
329
330
331 IF (GRANULAR_ENERGY) THEN
332 IF(.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID) THEN
333 CALL SOLVE_GRANULAR_ENERGY (IER)
334 IF (IER_MANAGER(IER)) goto 1000
335 ENDIF
336 ENDIF
337
338
339 CALL SOLVE_SPECIES_EQ (IER)
340 IF (IER_MANAGER(IER)) goto 1000
341
342
343 IF(NScalar /= 0) CALL SOLVE_Scalar_EQ (IER)
344
345
346 IF(K_Epsilon) CALL SOLVE_K_Epsilon_EQ (IER)
347
348
349
350
351 IF (.NOT.CYCLIC) LEQ_ADJUST = .TRUE.
352
353
354
355 CALL ACCUM_RESID
356 = RESID(RESID_P,0)
357 RESS = RESID(RESID_P,1)
358 CALL CALC_RESID_MB(1, errorpercent)
359 CALL CHECK_CONVERGENCE (NIT, errorpercent(0), MUSTIT, IER)
360
361 IF(CYCLIC)THEN
362 IF(MUSTIT==0 .OR. NIT >= MAX_NIT) &
363 CALL GoalSeekMassFlux(NIT, MUSTIT, .true.)
364 IF(AUTOMATIC_RESTART) RETURN
365 ENDIF
366
367
368
369 CONTINUE
370
371
372
373 IF (FULL_LOG) CALL DISPLAY_RESID (NIT, IER)
374
375
376 IF (MUSTIT == 0) THEN
377
378 IF (DT==UNDEFINED .AND. NIT==1) GOTO 50
379
380
381 IF (MOD(NSTEP,NLOG) == 0) THEN
382 CALL CPU_TIME (CPU_NOW)
383 CPUOS = (CPU_NOW - CPU_NLOG)/(TIME - TIME_NLOG)
384 CPU_NLOG = CPU_NOW
385 TIME_NLOG = TIME
386 CPU_NOW = CPU_NOW - CPU0
387
388 CALL CALC_RESID_MB(1, errorpercent)
389 CALL GET_SMASS (SMASS)
390 IF (ENERGY_EQ) CALL GET_HLOSS (HLOSS)
391
392 CALL START_LOG
393 IF (ENERGY_EQ) THEN
394 IF(DMP_LOG)WRITE (UNIT_LOG, 5000) TIME, DT, NIT, SMASS,&
395 HLOSS, CPU_NOW
396 IF(FULL_LOG.and.myPE.eq.PE_IO) &
397 WRITE(*,5000)TIME,DT,NIT,SMASS, HLOSS,CPU_NOW
398 ELSE
399 IF(DMP_LOG)WRITE (UNIT_LOG, 5001) TIME, DT, NIT, &
400 SMASS, CPU_NOW
401 IF (FULL_LOG .and. myPE.eq.PE_IO) &
402 WRITE (*, 5001) TIME, DT, NIT, SMASS, CPU_NOW
403 ENDIF
404
405 IF(DMP_LOG)WRITE (UNIT_LOG, 5002) (errorpercent(M), M=0,MMAX)
406 IF (FULL_LOG .and. myPE.eq.PE_IO) &
407 WRITE (*, 5002) (errorpercent(M), M=0,MMAX)
408 IF (.NOT.FULL_LOG) THEN
409 TLEFT = (TSTOP - TIME)*CPUOS
410 CALL GET_TUNIT (TLEFT, TUNIT)
411 IF(DMP_LOG)WRITE (UNIT_LOG, '(46X,A,F9.3,1X,A)')
412 ENDIF
413
414 IF (CYCLIC_X .OR. CYCLIC_Y .OR. CYCLIC_Z) THEN
415 IF (DO_I) THEN
416 Vavg = VAVG_U_G()
417 IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'U_g = ', Vavg
418 ENDIF
419 IF (DO_J) THEN
420 Vavg = VAVG_V_G()
421 IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'V_g = ', Vavg
422 ENDIF
423 IF (DO_K) THEN
424 Vavg = VAVG_W_G()
425 IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'W_g = ', Vavg
426 ENDIF
427 DO M = 1, SMAX
428 IF (DO_I) Then
429 Vavg = VAVG_U_S(M)
430 IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'U_s(', M, ') = ', Vavg
431 ENDIF
432 IF (DO_J) Then
433 Vavg = VAVG_V_S(M)
434 IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'V_s(', M, ') = ', Vavg
435 ENDIF
436 IF (DO_K) Then
437 Vavg = VAVG_W_S(M)
438 IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'W_s(', M, ') = ', Vavg
439 ENDIF
440 ENDDO
441 ENDIF
442
443 CALL END_LOG
444 ENDIF
445
446
447 IF(WRITE_DASHBOARD) THEN
448 RUN_STATUS = 'In Progress...'
449 N_DASHBOARD = N_DASHBOARD + 1
450 IF(MOD(N_DASHBOARD,F_DASHBOARD)==0) THEN
451 TLEFT = (TSTOP - TIME)*CPUOS
452 CALL GET_TUNIT (TLEFT, TUNIT)
453 CALL UPDATE_DASHBOARD(NIT,TLEFT,TUNIT)
454 ENDIF
455 ENDIF
456
457 IER = 0
458 RETURN
459
460
461
462
463 ELSEIF (MUSTIT==2 .AND. DT/=UNDEFINED) THEN
464
465 IF (FULL_LOG) THEN
466 CALL START_LOG
467 CALL CALC_RESID_MB(1, errorpercent)
468
469 IF(DMP_LOG) WRITE(UNIT_LOG,5200) TIME, DT, NIT, &
470 errorpercent(0), trim(adjustl(lMsg))
471 CALL END_LOG
472
473 IF (myPE.EQ.PE_IO) WRITE(*,5200) TIME, DT, NIT, &
474 errorpercent(0), trim(adjustl(lMsg))
475 ENDIF
476
477
478
479 IF(WRITE_DASHBOARD) THEN
480 RUN_STATUS = 'Diverged/stalled...'
481 N_DASHBOARD = N_DASHBOARD + 1
482 IF(MOD(N_DASHBOARD,F_DASHBOARD)==0) THEN
483 TLEFT = (TSTOP - TIME)*CPUOS
484 CALL GET_TUNIT (TLEFT, TUNIT)
485 CALL UPDATE_DASHBOARD(NIT,TLEFT,TUNIT)
486 ENDIF
487 ENDIF
488
489 IER = 1
490 RETURN
491 ENDIF
492
493
494
495
496
497 IF (NIT < MAX_NIT) THEN
498 MUSTIT = 0
499 GOTO 50
500 ENDIF
501
502
503
504
505
506 CALL GET_SMASS (SMASS)
507 IF (myPE.eq.PE_IO) WRITE(UNIT_OUT, 5100) TIME, DT, NIT, SMASS
508 CALL START_LOG
509 IF(DMP_LOG) WRITE(UNIT_LOG, 5100) TIME, DT, NIT, SMASS
510 CALL END_LOG
511
512
513
514 = 1
515 RETURN
516
517
518 5000 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',I3,' Sm=',G12.5,' Hl=',G12.5,&
519 T84,'CPU=',F8.0,' s')
520 5001 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',I3,' Sm=',G12.5, T84,'CPU=',F8.0,' s')
521 5002 FORMAT(3X,'MbError%(0,MMAX):', 5(1X,G11.4))
522 5050 FORMAT(5X,'Average ',A,G12.5)
523 5060 FORMAT(5X,'Average ',A,I2,A,G12.5)
524 5100 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT>',I3,' Sm= ',G12.5, 'MbErr%=', G11.4)
525 5200 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',&
526 I3,'MbErr%=', G11.4, ': ',A,' :-(')
527 6000 FORMAT(1X,A)
528
529
530 contains
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547 LOGICAL FUNCTION IER_MANAGER(lErr)
548
549 INTEGER, intent(in) :: lErr
550
551
552 IF(IER < 100) THEN
553 IER_MANAGER = .FALSE.
554 return
555 ENDIF
556
557
558
559 = .TRUE.
560 MUSTIT = 2
561
562
563
564 IF(IER < 110) THEN
565 IF(IER == 100) THEN
566 lMsg = 'Negative gas density detected'
567 ELSEIF(IER == 101) THEN
568 lMsg = 'Negative solids density detected'
569 ELSE
570 lMsg = 'UCE in PHYSICAL_PROP'
571 ENDIF
572
573
574
575
576 ELSEIF(IER < 120) THEN
577 IF(IER == 110) THEN
578 lMsg = 'Negative void fraction detected'
579 ELSE
580 lMsg = 'UCE in CALC_VOL_FR'
581 ENDIF
582
583
584
585
586 ELSEIF(IER < 130) THEN
587 IF(IER == 120) THEN
588 lMsg = 'Energy Equation diverged'
589 ELSE
590 lMsg = 'UCE in SOLVE_ENERGY_EQ'
591 ENDIF
592
593
594
595
596 ELSEIF(IER < 140) THEN
597 IF(IER == 130) THEN
598 lMsg = 'Species Equation diverged'
599 ELSE
600 lMsg = 'UCE in SOLVE_SPECIES_EQ'
601 ENDIF
602
603
604
605
606 ELSEIF(IER < 150) THEN
607 IF(IER == 140) THEN
608 lMsg = 'Granular Energy Eq diverged'
609 ELSE
610 lMsg = 'UCE in SOLVE_GRANULAR_ENERGY'
611 ENDIF
612
613
614
615 ELSE
616 lMsg = 'Run diverged/stalled with UCE'
617 ENDIF
618
619
620 IF(DT == UNDEFINED) IER_MANAGER = .FALSE.
621
622 return
623 END FUNCTION IER_MANAGER
624
625
626
627 END SUBROUTINE ITERATE
628
629
630
631
632
633
634
635 SUBROUTINE GET_TUNIT(TLEFT, TUNIT)
636
637
638
639
640 IMPLICIT NONE
641
642
643
644 DOUBLE PRECISION, INTENT(INOUT) :: TLEFT
645 CHARACTER(LEN=4) :: TUNIT
646
647
648 IF (TLEFT < 3600.0d0) THEN
649 TUNIT = 's'
650 ELSE
651 TLEFT = TLEFT/3600.0d0
652 TUNIT = 'h'
653 IF (TLEFT >= 24.) THEN
654 TLEFT = TLEFT/24.0d0
655 TUNIT = 'days'
656 ENDIF
657 ENDIF
658
659 RETURN
660 END SUBROUTINE GET_TUNIT
661
662
663
664
665
666
667
668
669
670
671
672 subroutine GoalSeekMassFlux(NIT, MUSTIT, doit)
673
674
675
676
677 USE bc
678 USE geometry
679 USE constant
680 USE compar
681 USE run
682 USE time_cpu
683 IMPLICIT NONE
684
685
686
687 INTEGER, INTENT(INOUT) :: NIT, MUSTIT
688 LOGICAL, INTENT(IN) :: doit
689
690
691
692 INTEGER, PARAMETER :: MAXOUTIT = 500
693 DOUBLE PRECISION, PARAMETER :: omega = 0.9
694 DOUBLE PRECISION, PARAMETER :: TOL = 1E-03
695 INTEGER, SAVE :: OUTIT
696 LOGICAL, SAVE :: firstPass = .true.
697
698 DOUBLE PRECISION, SAVE :: mdot_n, mdot_nm1, delp_n, delp_nm1, err
699 DOUBLE PRECISION :: mdot_0, delp_xyz
700
701 CHARACTER, SAVE :: Direction
702
703
704
705 DOUBLE PRECISION, EXTERNAL :: VAVG_Flux_U_G, VAVG_Flux_V_G, &
706 VAVG_Flux_W_G
707 LOGICAL, EXTERNAL :: IsNan
708
709
710 IF(CYCLIC_X_MF)THEN
711 delp_n = delp_x
712 ELSEIF(CYCLIC_Y_MF)THEN
713 delp_n = delp_y
714 ELSEIF(CYCLIC_Z_MF)THEN
715 delp_n = delp_z
716 ELSE
717 RETURN
718 ENDIF
719
720 IF(.NOT.doit) THEN
721 OUTIT = 0
722 RETURN
723 ENDIF
724
725 OUTIT = OUTIT + 1
726 IF(OUTIT > MAXOUTIT) THEN
727 IF (myPE.EQ.PE_IO) write(*,5400) MAXOUTIT
728 CALL mfix_exit(0)
729 ENDIF
730
731 mdot_0 = Flux_g
732
733
734
735 IF(CYCLIC_X_MF)THEN
736 mdot_n = VAVG_Flux_U_G()
737 ELSEIF(CYCLIC_Y_MF)THEN
738 mdot_n = VAVG_Flux_V_G()
739 ELSEIF(CYCLIC_Z_MF)THEN
740 mdot_n = VAVG_Flux_W_G()
741 ENDIF
742
743 IF (isNan(mdot_n) .OR. isNan(delp_n)) THEN
744 IF (myPE.eq.PE_IO) write(*,*) mdot_n, delp_n, &
745 ' NaN being caught in GoalSeekMassFlux '
746 AUTOMATIC_RESTART = .TRUE.
747 RETURN
748 ENDIF
749
750 err = abs((mdot_n - mdot_0)/mdot_0)
751 IF( err < TOL) THEN
752 MUSTIT = 0
753 ELSE
754 MUSTIT = 1
755 NIT = 1
756 ENDIF
757
758
759 if(.not.firstPass)then
760
761
762
763
764
765 = delp_n - omega * (delp_n - delp_nm1) * &
766 ((mdot_n - mdot_0)/(mdot_nm1 - mdot_0)) / &
767 ((mdot_n - mdot_0)/(mdot_nm1 - mdot_0) - ONE)
768 else
769 firstPass=.false.
770 delp_xyz = delp_n*0.99
771 endif
772
773 IF(MUSTIT == 0) then
774 IF(myPE.eq.PE_IO) Write(*,5500) TIME, OUTIT, delp_xyz, mdot_n
775 ENDIF
776
777 mdot_nm1 = mdot_n
778 delp_nm1 = delp_n
779
780 IF(CYCLIC_X_MF)THEN
781 delp_x = delp_xyz
782 ELSEIF(CYCLIC_Y_MF)THEN
783 delp_y = delp_xyz
784 ELSEIF(CYCLIC_Z_MF)THEN
785 delp_z = delp_xyz
786 ENDIF
787
788
789 RETURN
790
791 5400 FORMAT(/1X,70('*')//' From: GoalSeekMassFlux',/&
792 ' Message: Number of outer iterations exceeded ', I4,/1X,70('*')/)
793 5500 Format(' Time=', G12.5, ' MassFluxIterations=', I4, ' DelP=', &
794 G12.5, ' Gas Flux=', G12.5)
795
796 END SUBROUTINE GoalSeekMassFlux
797
798