File: N:\mfix\model\iterate.f

1     ! -*- f90 -*-
2           MODULE ITERATE
3     
4           USE EXIT, ONLY: MFIX_EXIT
5           USE RUN, ONLY: IER, TUNIT, TIME
6     
7     ! flag indicating convergence status with MUSTIT = 0,1,2 implying
8     ! complete convergence, non-covergence and divergence respectively
9           INTEGER :: MUSTIT
10     
11     ! Number of iterations completed for current timestep
12           INTEGER :: NIT
13     
14     ! User defined maximum number of iterations
15           INTEGER :: MAX_NIT
16     
17           LOGICAL :: CONVERGED, DIVERGED
18     
19     ! cpu time left
20           DOUBLE PRECISION :: TLEFT
21     ! Normalization factor for gas & solids pressure residual
22           DOUBLE PRECISION :: NORMg, NORMs
23     ! Set normalization factor for gas and solids pressure residual
24           LOGICAL :: SETg, SETs
25     ! gas & solids pressure residual
26           DOUBLE PRECISION :: RESg, RESs
27     ! Weight of solids in the reactor
28           DOUBLE PRECISION :: SMASS
29           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: errorpercent
30     ! Error Message
31           CHARACTER(LEN=32) :: lMsg
32     
33           CONTAINS
34     
35     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
36     !                                                                      !
37     !  Subroutine: ITERATE                                                 !
38     !  Author: M. Syamlal                                 Date: 12-APR-96  !
39     !                                                                      !
40     !  Purpose: This module controls the iterations for solving equations  !
41     !                                                                      !
42     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
43           SUBROUTINE ITERATE_INIT
44     
45           USE compar, only: mype, pe_io
46           USE cutcell, only: cartesian_grid
47           USE geometry, only: cyclic
48           USE leqsol, only: leq_adjust
49           USE output, only: full_log
50           USE param1, only: one, small_number, undefined, zero
51           USE run, only: dt, dt_prev, run_type, time, tstop, nstep, nsteprst, cn_on, get_tunit, steady_state
52           USE time_cpu
53           USE toleranc, only: norm_g, norm_s
54     
55           IMPLICIT NONE
56     
57     ! initializations
58           DT_prev = DT
59           NIT = 0
60           MUSTIT = 0
61           CONVERGED = .FALSE.
62           DIVERGED  = .FALSE.
63           RESG = ZERO
64           RESS = ZERO
65     
66           IF(NORM_G == UNDEFINED) THEN
67              NORMG = ONE
68              SETG = .FALSE.
69           ELSE
70              NORMG = NORM_G
71              SETG = .TRUE.
72           ENDIF
73     
74           IF(NORM_S == UNDEFINED) THEN
75              NORMS = ONE
76              SETS = .FALSE.
77           ELSE
78              NORMS = NORM_S
79              SETS = .TRUE.
80           ENDIF
81     
82           LEQ_ADJUST = .FALSE.
83     
84     ! Initialize residuals
85           CALL INIT_RESID ()
86     
87     ! Initialize the routine for holding gas mass flux constant with cyclic bc
88           IF(CYCLIC) CALL GoalSeekMassFlux(NIT, MUSTIT, .false.)
89     
90     ! CPU time left
91           IF (FULL_LOG) THEN
92              TLEFT = (TSTOP - TIME)*CPUOS
93              CALL GET_TUNIT (TLEFT, TUNIT)
94     
95              IF (STEADY_STATE) THEN
96                 CALL GET_SMASS (SMASS)
97                 IF(myPE.eq.PE_IO) THEN
98                    WRITE (*, '(/A,G12.5, A,F9.3,1X,A)') &
99                       ' Starting solids mass = ', SMASS
100                 ENDIF
101              ELSE
102                 IF(myPE.eq.PE_IO) THEN
103                    IF ((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
104                        (CN_ON.AND.RUN_TYPE /= 'NEW' .AND.&
105                         NSTEP >= (NSTEPRST+1))) THEN
106                       WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)')&
107                          ' Time = ', TIME, '  Dt = ', 2.D0*DT
108                    ELSE
109                       WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)') &
110                          ' Time = ', TIME, '  Dt = ', DT
111                    ENDIF
112                 ENDIF
113              ENDIF   ! if/else(steady_state)
114           ENDIF   ! if(full_log)
115     
116           CALL CALC_RESID_MB(0, errorpercent)
117     
118     ! Calculate the face values of densities and mass fluxes for the first
119     ! solve_vel_star call.
120           CALL CONV_ROP()
121           CALL CALC_MFLUX ()
122           CALL SET_BC1
123           CALL SET_EP_FACTORS
124     
125     ! JFD: modification for cartesian grid implementation
126           IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
127     
128     ! Default/Generic Error message
129           lMsg = 'Run diverged/stalled'
130     
131           RETURN
132           END SUBROUTINE ITERATE_INIT
133     
134     
135     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
136     !                                                                      !
137     !  Subroutine: ITERATE                                                 !
138     !  Author: M. Syamlal                                 Date: 12-APR-96  !
139     !                                                                      !
140     !  Purpose: This module controls the iterations for solving equations  !
141     !                                                                      !
142     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
143           SUBROUTINE DO_ITERATION
144     
145           USE cont, only: solve_continuity
146              USE cutcell, only: cartesian_grid
147              USE discretelement, only: discrete_element, des_continuum_hybrid
148              USE fldvar, only: ep_g, ro_g, rop_g, rop_s, p_star
149              USE geometry, only: cyclic, cylindrical
150              USE leqsol, only: leq_adjust
151              USE mms, only: USE_MMS
152              USE param1, only: small_number, zero, undefined, undefined_i
153              USE physprop, only: mmax, ro_g0, smax
154              USE pscor, only: k_cp, mcp
155              USE qmom_kinetic_equation, only: qmomk
156              USE residual, only: resid_p, resid
157              USE run, only: auto_restart, automatic_restart, call_dqmom, call_usr, chk_batchq_end, friction, ghd_2007, granular_energy, k_epsilon, kt_type_enum, phip_out_iter, energy_eq, dt, ier, steady_state
158              USE scalars, only: nscalar
159     
160           IMPLICIT NONE
161     
162           INTEGER :: M
163     
164           PHIP_OUT_ITER=NIT ! To record the output of phip
165     ! mechanism to set the normalization factor for the correction
166     ! after the first iteration to the corresponding residual found
167     ! in the first iteration
168           IF (.NOT.SETG) THEN
169              IF (RESG > SMALL_NUMBER) THEN
170                 NORMG = RESG
171                 SETG = .TRUE.
172              ENDIF
173           ENDIF
174           IF (.NOT.SETS) THEN
175              IF (RESS > SMALL_NUMBER) THEN
176                 NORMS = RESS
177                 SETS = .TRUE.
178              ENDIF
179           ENDIF
180     
181     ! Call user-defined subroutine to set quantities that need to be updated
182     ! every iteration
183           IF (CALL_USR) CALL USR2
184     
185     ! Calculate coefficients, excluding density and reactions.
186           CALL CALC_COEFF(IER, 1)
187           IF (IER_MANAGER()) return
188     
189     ! Diffusion coefficient and source terms for user-defined scalars
190           IF(NScalar /= 0) CALL SCALAR_PROP()
191     
192     ! Diffusion coefficient and source terms for K & Epsilon Eq.
193           IF(K_Epsilon) CALL K_Epsilon_PROP()
194     
195     ! Update the stress tensor trace and cross terms each subiteration
196     ! for MMS cases.
197           IF(USE_MMS) CALL CALC_TRD_AND_TAU()
198     
199     ! Solve starred velocity components
200           CALL SOLVE_VEL_STAR(IER)
201     
202     ! Correct the centerline velocity for cylindrical simulations.
203           IF(CYLINDRICAL) CALL RADIAL_VEL_CORRECTION
204     
205     ! Calculate densities.
206           CALL PHYSICAL_PROP(IER, 0)
207           IF (IER_MANAGER()) return
208     
209     ! Calculate chemical reactions.
210           CALL CALC_RRATE(IER)
211     
212     ! Solve solids volume fraction correction equation for close-packed
213     ! solids phases
214           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
215              DES_CONTINUUM_HYBRID) THEN
216              IF (MMAX > 0) THEN
217     ! MMS:  Solve gas continuity only.
218                 IF(USE_MMS) THEN
219                    CALL SOLVE_CONTINUITY(0,IER)
220     ! Regular, non-MMS cases.
221                 ELSE
222                    IF(MMAX == 1 .AND. MCP /= UNDEFINED_I)THEN
223     ! if second phase (m=1) can overpack (e.g., bubbles) then solve its
224     ! continuity equation
225                       CALL CALC_K_CP (K_CP)
226                       CALL SOLVE_EPP (NORMS, RESS, IER)
227                       CALL CORRECT_1 ()
228                    ELSE
229     
230     ! If one chooses to revert back to old mark_phase_4_cor wherein the
231     ! continuity of the gas phase can get marked to be solved then this
232     ! loop should start at 0.
233                       DO M=1,SMAX ! mmax -> smax for GHD theory
234     ! Volume fraction correction technique for one of the solids phase
235     ! is not implemented.  This will only slow down convergence.
236     !                      IF (M .EQ. MCP) THEN
237     !                       CALL CALC_K_CP (K_CP, IER)
238     !                       CALL SOLVE_EPP (NORMS, RESS, IER)
239     !                       CALL CORRECT_1 (IER)
240     !                    ELSE
241                             CALL SOLVE_CONTINUITY(M,IER)
242     !                    ENDIF
243                       ENDDO
244                    ENDIF   ! end if/else (mmax==1 .and. mcp /= undefined)
245                 ENDIF ! end if/else (MMS)
246     
247                 IF(KT_TYPE_ENUM == GHD_2007) CALL ADJUST_EPS_GHD
248     
249                 CALL CALC_VOL_FR (P_STAR, RO_G, ROP_G, EP_G, ROP_S, IER)
250                 IF (IER_MANAGER()) return
251     
252              ENDIF  ! endif (mmax >0)
253     
254           ENDIF  ! end if (.not.discrete_element)
255     
256     
257     ! Calculate P_star in cells where solids continuity equation is
258     ! solved
259           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
260              DES_CONTINUUM_HYBRID) THEN
261              IF (MMAX > 0 .AND. .NOT.FRICTION) &
262                 CALL CALC_P_STAR (EP_G, P_STAR)
263           ENDIF
264     
265     ! Calculate the face values of densities.
266           CALL CONV_ROP()
267     
268           IF (RO_G0 /= ZERO) THEN
269     ! Solve fluid pressure correction equation
270     #ifndef FLAG_MMS
271              CALL SOLVE_PP_G (NORMG, RESG, IER)
272     #endif
273     ! Correct pressure, velocities, and density
274              CALL CORRECT_0 ()
275           ENDIF
276     
277     ! Recalculate densities.
278           CALL PHYSICAL_PROP(IER, 0)
279           IF (IER_MANAGER()) return
280     
281     ! Update wall velocities:
282     ! modified by sof to force wall functions so even when NSW or FSW are
283     ! declared, default wall BC will still be treated as NSW and no wall
284     ! functions will be used
285           IF(.NOT. K_EPSILON) CALL SET_WALL_BC ()
286     
287     ! Calculate the face values of mass fluxes
288           CALL CALC_MFLUX ()
289           CALL SET_BC1
290           CALL SET_EP_FACTORS
291     
292     ! JFD: modification for cartesian grid implementation
293           IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
294     
295     ! Solve energy equations
296           IF (ENERGY_EQ) THEN
297              CALL SOLVE_ENERGY_EQ (IER)
298              IF (IER_MANAGER()) return
299           ENDIF
300     
301     ! Solve granular energy equation
302           IF (GRANULAR_ENERGY) THEN
303              IF(.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID) THEN
304                 CALL SOLVE_GRANULAR_ENERGY (IER)
305                 IF (IER_MANAGER()) return
306              ENDIF
307           ENDIF
308     
309     ! Solve species mass balance equations.
310           CALL SOLVE_SPECIES_EQ (IER)
311           IF (IER_MANAGER()) return
312     
313     ! Solve other scalar transport equations
314           IF(NScalar /= 0) CALL SOLVE_Scalar_EQ (IER)
315     
316     ! Solve K & Epsilon transport equations
317           IF(K_Epsilon) CALL SOLVE_K_Epsilon_EQ (IER)
318     
319     ! User-defined linear equation solver parameters may be adjusted after
320     ! the first iteration
321           IF (.NOT.CYCLIC) LEQ_ADJUST = .TRUE.
322     
323     ! Check for convergence
324           CALL ACCUM_RESID ! Accumulating residuals from all the processors
325           RESG = RESID(RESID_P,0)
326           RESS = RESID(RESID_P,1)
327           CALL CALC_RESID_MB(1, errorpercent)
328           MUSTIT = 0
329           CALL CHECK_CONVERGENCE (NIT, errorpercent(0), MUSTIT)
330     
331           IF(CYCLIC .AND. (MUSTIT==0 .OR. NIT >= MAX_NIT)) &
332              CALL GoalSeekMassFlux(NIT, MUSTIT, .true.)
333     
334     ! Display residuals
335           CALL DISPLAY_RESID (NIT)
336     
337           CALL END_ITERATION
338     
339           contains
340     
341           SUBROUTINE END_ITERATION
342              IMPLICIT NONE
343     
344              ! Determine course of simulation: converge, non-converge, diverge?
345              IF (MUSTIT == 0) THEN
346                 IF (STEADY_STATE .AND. NIT==1) RETURN   !Iterations converged
347                 CONVERGED = .TRUE.
348                 IER = 0
349              ELSEIF (MUSTIT==2 .AND. .NOT.STEADY_STATE) THEN
350                 DIVERGED = .TRUE.
351                 IER = 1
352              ENDIF
353     
354           END SUBROUTINE END_ITERATION
355     
356     !----------------------------------------------------------------------!
357     ! Function: IER_Manager                                                !
358     !                                                                      !
359     ! Purpose: Identify and account for errors from called subroutines.    !
360     !          Returns .TRUE. for lErr >= 100, otherwise .FALSE.           !
361     !                                                                      !
362     ! Reserved Error Blocks:                                               !
363     !                                                                      !
364     ! [ 100,  109]: PHYSICAL_PROP                                          !
365     ! [ 110,  119]: CALC_VOL_FR                                            !
366     ! [ 120,  129]: SOLVE_ENERGY_EQ                                        !
367     ! [ 130,  139]: SOLVE_SPECIES_EQ                                       !
368     ! [ 140,  149]: SOLVE_GRANULAR_ENERGY                                  !
369     !                                                                      !
370     !----------------------------------------------------------------------!
371           LOGICAL FUNCTION IER_MANAGER()
372     
373     ! Default case: do nothing.
374           IF(IER < 100) THEN
375              IER_MANAGER = .FALSE.
376              return
377           ENDIF
378     
379     ! Errors with an index greater than 100 will force an exit from iterate
380     ! and in turn, reduce the step-size, and restart the time-step.
381           IER_MANAGER = .TRUE.
382           MUSTIT = 2
383     
384     ! Errors reported from PHYSICAL_PROP
385     !```````````````````````````````````````````````````````````````````````
386           IF(IER <  110) THEN
387              IF(IER ==  100) THEN
388                 lMsg = 'Negative gas density detected'
389              ELSEIF(IER ==  101) THEN
390                 lMsg = 'Negative solids density detected'
391              ELSE
392                 lMsg = 'UCE in PHYSICAL_PROP'
393              ENDIF
394     
395     
396     ! Errors reported from CALC_VOL_FR
397     !```````````````````````````````````````````````````````````````````````
398           ELSEIF(IER <  120) THEN
399              IF(IER ==  110) THEN
400                 lMsg = 'Negative void fraction detected'
401              ELSE
402                 lMsg = 'UCE in CALC_VOL_FR'
403              ENDIF
404     
405     
406     ! Errors reported from SOLVE_ENERGY_EQ
407     !```````````````````````````````````````````````````````````````````````
408           ELSEIF(IER <  130) THEN
409              IF(IER ==  120) THEN
410                 lMsg = 'Energy Equation diverged'
411              ELSE
412                 lMsg = 'UCE in SOLVE_ENERGY_EQ'
413              ENDIF
414     
415     
416     ! Errors reported from SOLVE_SPECIES_EQ
417     !```````````````````````````````````````````````````````````````````````
418           ELSEIF(IER <  140) THEN
419              IF(IER ==  130) THEN
420                 lMsg = 'Species Equation diverged'
421              ELSE
422                 lMsg = 'UCE in SOLVE_SPECIES_EQ'
423              ENDIF
424     
425     
426     ! Errors reported from SOLVE_GRANULAR_ENERGY
427     !```````````````````````````````````````````````````````````````````````
428           ELSEIF(IER <  150) THEN
429              IF(IER ==  140) THEN
430                 lMsg = 'Granular Energy Eq diverged'
431              ELSE
432                 lMsg = 'UCE in SOLVE_GRANULAR_ENERGY'
433              ENDIF
434     
435     ! Unclassified Errors
436     !```````````````````````````````````````````````````````````````````````
437           ELSE
438              lMsg = 'Run diverged/stalled with UCE'
439           ENDIF
440     
441     
442           IF(STEADY_STATE) IER_MANAGER = .FALSE.
443     
444           CALL END_ITERATION
445     
446           return
447           END FUNCTION IER_MANAGER
448     
449           END SUBROUTINE DO_ITERATION
450     
451     
452     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
453     !                                                                      !
454     !  Subroutine: ITERATE                                                 !
455     !  Author: M. Syamlal                                 Date: 12-APR-96  !
456     !                                                                      !
457     !  Purpose: This module controls the iterations for solving equations  !
458     !                                                                      !
459     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
460           SUBROUTINE POST_ITERATE
461     
462           USE compar, only: mype, pe_io
463           USE funits, only: dmp_log, unit_log, unit_out
464           USE machine, only: start_log, end_log
465           USE run, only: dt, time
466     
467           IMPLICIT NONE
468     
469           IF (CONVERGED) THEN
470              CALL LOG_CONVERGED
471           ELSEIF (DIVERGED) THEN
472              CALL LOG_DIVERGED
473           ELSE
474              CALL GET_SMASS (SMASS)
475              IF (myPE.eq.PE_IO) WRITE(UNIT_OUT, 5100) TIME, DT, NIT, SMASS
476              CALL START_LOG
477              IF(DMP_LOG) WRITE(UNIT_LOG, 5100) TIME, DT, NIT, SMASS
478              CALL END_LOG
479           ENDIF
480     
481           IF(.NOT.(CONVERGED .OR. DIVERGED)) THEN
482              IER = 1
483           ENDIF
484     
485     5100  FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT>',I3,' Sm= ',G12.5, &
486                'MbErr%=', G11.4)
487     
488           END SUBROUTINE POST_ITERATE
489     
490     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
491     !                                                                      !
492     !  Subroutine: ITERATE                                                 !
493     !  Author: M. Syamlal                                 Date: 12-APR-96  !
494     !                                                                      !
495     !  Purpose: This module controls the iterations for solving equations  !
496     !                                                                      !
497     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
498           SUBROUTINE LOG_DIVERGED
499           USE compar, only: mype, pe_io
500           USE dashboard, only: f_dashboard, n_dashboard, run_status, write_dashboard
501           USE funits, only: dmp_log, unit_log
502           USE machine, only: start_log, end_log
503           USE output, only: full_log
504           USE run, only: dt, time, tstop, get_tunit
505           USE time_cpu, only: cpuos
506     
507           IMPLICIT NONE
508     
509           CHARACTER(LEN=4) :: TUNIT
510     
511           IF (FULL_LOG) THEN
512              CALL START_LOG
513              CALL CALC_RESID_MB(1, errorpercent)
514     
515              IF(DMP_LOG) WRITE(UNIT_LOG,5200) TIME, DT, NIT, &
516                   errorpercent(0), trim(adjustl(lMsg))
517              CALL END_LOG
518     
519              IF (myPE.EQ.PE_IO) WRITE(*,5200) TIME, DT, NIT, &
520                   errorpercent(0), trim(adjustl(lMsg))
521           ENDIF
522     
523           ! JFD: modification for cartesian grid implementation
524           IF(WRITE_DASHBOARD) THEN
525              RUN_STATUS = 'Diverged/stalled...'
526              N_DASHBOARD = N_DASHBOARD + 1
527              IF(MOD(N_DASHBOARD,F_DASHBOARD)==0) THEN
528                 TLEFT = (TSTOP - TIME)*CPUOS
529                 CALL GET_TUNIT (TLEFT, TUNIT)
530                 CALL UPDATE_DASHBOARD(NIT,TLEFT,TUNIT)
531              ENDIF
532           ENDIF
533     5200  FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',&
534                      I3,'MbErr%=', G11.4, ': ',A,' :-(')
535           END SUBROUTINE LOG_DIVERGED
536     
537     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
538     !                                                                      !
539     !  Subroutine: ITERATE                                                 !
540     !  Author: M. Syamlal                                 Date: 12-APR-96  !
541     !                                                                      !
542     !  Purpose: This module controls the iterations for solving equations  !
543     !                                                                      !
544     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
545           SUBROUTINE LOG_CONVERGED
546     
547           USE compar, only: mype, pe_io
548           USE dashboard, only: f_dashboard, n_dashboard, run_status, write_dashboard
549           USE error_manager, only: err_msg
550           USE funits, only: dmp_log, unit_log
551           USE geometry, only: cyclic_x, cyclic_y, cyclic_z
552           USE geometry, only: do_i, do_j, do_k
553           USE machine, only: start_log, end_log
554           USE output, only: full_log, nlog
555           USE physprop, only: mmax, smax
556           USE run, only: dt, energy_eq, time, tstop, nstep, get_tunit
557           USE time_cpu, only: cpu0, cpu_nlog, cpuos, time_nlog
558     
559           IMPLICIT NONE
560     
561           CHARACTER(LEN=4) :: TUNIT
562     ! Perform checks and dump to screen every NLOG time steps
563           IF (MOD(NSTEP,NLOG) == 0) CALL DUMP_TO_SCREEN
564     
565     ! JFD: modification for cartesian grid implementation
566           IF(WRITE_DASHBOARD) THEN
567              RUN_STATUS = 'In Progress...'
568              N_DASHBOARD = N_DASHBOARD + 1
569              IF(MOD(N_DASHBOARD,F_DASHBOARD)==0) THEN
570                 TLEFT = (TSTOP - TIME)*CPUOS
571                 CALL GET_TUNIT (TLEFT, TUNIT)
572                 CALL UPDATE_DASHBOARD(NIT,TLEFT,TUNIT)
573              ENDIF
574           ENDIF
575     
576           CONTAINS
577     
578           SUBROUTINE DUMP_TO_SCREEN
579           USE error_manager, only: flush_err_msg
580           IMPLICIT NONE
581     
582           ! phase index
583           INTEGER :: M
584     ! current cpu time used
585           DOUBLE PRECISION :: CPU_NOW
586     ! Heat loss from the reactor
587           DOUBLE PRECISION :: HLOSS
588     ! average velocity
589           DOUBLE PRECISION :: Vavg
590     
591     !-----------------------------------------------
592     ! External functions
593     !-----------------------------------------------
594           DOUBLE PRECISION, EXTERNAL :: VAVG_U_G, VAVG_V_G, VAVG_W_G, &
595              VAVG_U_S, VAVG_V_S, VAVG_W_S
596     
597           CALL CPU_TIME (CPU_NOW)
598           CPUOS = (CPU_NOW - CPU_NLOG)/(TIME - TIME_NLOG)
599           CPU_NLOG = CPU_NOW
600           TIME_NLOG = TIME
601           CPU_NOW = CPU_NOW - CPU0
602     
603           CALL CALC_RESID_MB(1, errorpercent)
604           CALL GET_SMASS (SMASS)
605           IF (ENERGY_EQ) CALL GET_HLOSS (HLOSS)
606     
607           CALL START_LOG
608           IF (ENERGY_EQ) THEN
609              WRITE(ERR_MSG,5000)TIME, DT, NIT, SMASS, HLOSS, CPU_NOW
610              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
611           ELSE
612              WRITE(ERR_MSG,5001) TIME, DT, NIT, SMASS, CPU_NOW
613              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
614           ENDIF
615     
616      5000 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',I3,' Sm=',G12.5, &
617              ' Hl=',G12.5,T84,'CPU=',F8.0,' s')
618     
619      5001 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',I3,' Sm=',G12.5, &
620              T84,'CPU=',F8.0,' s')
621     
622           IF(DMP_LOG)WRITE (UNIT_LOG, 5002) (errorpercent(M), M=0,MMAX)
623           IF (FULL_LOG .and. myPE.eq.PE_IO) &
624              WRITE (*, 5002) (errorpercent(M), M=0,MMAX)
625     
626      5002 FORMAT(3X,'MbError%(0,MMAX):', 5(1X,G11.4))
627     
628           IF (.NOT.FULL_LOG) THEN
629              TLEFT = (TSTOP - TIME)*CPUOS
630              CALL GET_TUNIT (TLEFT, TUNIT)
631              IF(DMP_LOG)WRITE (UNIT_LOG, '(46X,A,F9.3,1X,A)')
632           ENDIF
633     
634           IF (CYCLIC_X .OR. CYCLIC_Y .OR. CYCLIC_Z) THEN
635              IF (DO_I) THEN
636                Vavg = VAVG_U_G()
637                IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'U_g = ', Vavg
638              ENDIF
639              IF (DO_J) THEN
640                Vavg = VAVG_V_G()
641                IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'V_g = ',  Vavg
642              ENDIF
643              IF (DO_K) THEN
644                Vavg = VAVG_W_G()
645                IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'W_g = ', Vavg
646              ENDIF
647              DO M = 1, SMAX
648                 IF (DO_I) Then
649                   Vavg = VAVG_U_S(M)
650                   IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'U_s(', M, ') = ', Vavg
651                 ENDIF
652                 IF (DO_J) Then
653                   Vavg = VAVG_V_S(M)
654                   IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'V_s(', M, ') = ', Vavg
655                 ENDIF
656                 IF (DO_K) Then
657                   Vavg = VAVG_W_S(M)
658                   IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'W_s(', M, ') = ', Vavg
659                 ENDIF
660              ENDDO
661           ENDIF   ! end if cyclic_x, cyclic_y or cyclic_z
662     
663           CALL END_LOG
664     
665     5050  FORMAT(5X,'Average ',A,G12.5)
666     5060  FORMAT(5X,'Average ',A,I2,A,G12.5)
667     
668           END SUBROUTINE DUMP_TO_SCREEN
669     
670           END SUBROUTINE LOG_CONVERGED
671     
672     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
673     !                                                                      !
674     !  Subroutine: GoalSeekMassFlux                                        !
675     !  Author: M. Syamlal                                 Date: 12-APR-96  !
676     !                                                                      !
677     ! Purpose:  In the following subroutine the mass flux across a periodic!
678     ! domain with pressure drop is held constant at a user-specified value.!
679     ! This module is activated only if the user specifies a value for the  !
680     ! keyword flux_g in the mfix.dat file.                                 !
681     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
682           SUBROUTINE GoalSeekMassFlux(NIT, MUSTIT, doit)
683     
684     !-----------------------------------------------
685     ! Modules
686     !-----------------------------------------------
687           USE bc
688           USE compar
689           USE constant
690           USE exit, only: mfix_exit
691           USE geometry
692           USE param1, only: one
693           USE run
694           USE time_cpu
695           USE utilities, ONLY: mfix_isnan
696           IMPLICIT NONE
697     !-----------------------------------------------
698     ! Dummy arguments
699     !-----------------------------------------------
700           INTEGER, INTENT(INOUT) :: NIT, MUSTIT
701           LOGICAL, INTENT(IN) :: doit
702     !-----------------------------------------------
703     ! Local Variables
704     !-----------------------------------------------
705           INTEGER, PARAMETER :: MAXOUTIT = 500
706           DOUBLE PRECISION, PARAMETER          :: omega = 0.9
707           DOUBLE PRECISION, PARAMETER          :: TOL = 1E-03
708           INTEGER, SAVE :: OUTIT
709           LOGICAL, SAVE :: firstPass = .true.
710     
711           DOUBLE PRECISION, SAVE  :: mdot_n, mdot_nm1, delp_n, delp_nm1, err
712           DOUBLE PRECISION        :: mdot_0, delp_xyz
713     
714     !-----------------------------------------------
715     ! Functions
716     !-----------------------------------------------
717           DOUBLE PRECISION, EXTERNAL :: VAVG_Flux_U_G, VAVG_Flux_V_G, &
718                                         VAVG_Flux_W_G
719     !-----------------------------------------------
720     
721           IF(CYCLIC_X_MF)THEN
722              delp_n = delp_x
723           ELSEIF(CYCLIC_Y_MF)THEN
724              delp_n = delp_y
725           ELSEIF(CYCLIC_Z_MF)THEN
726              delp_n = delp_z
727           ELSE
728              RETURN
729           ENDIF
730     
731           IF(.NOT.doit) THEN
732              OUTIT = 0
733              RETURN
734           ENDIF
735     
736           OUTIT = OUTIT + 1
737           IF(OUTIT > MAXOUTIT) THEN
738              IF (myPE.EQ.PE_IO) write(*,5400) MAXOUTIT
739              CALL mfix_exit(0)
740           ENDIF
741     
742           mdot_0 = Flux_g
743     
744     
745           ! calculate the average gas mass flux and error
746           IF(CYCLIC_X_MF)THEN
747             mdot_n = VAVG_Flux_U_G()
748           ELSEIF(CYCLIC_Y_MF)THEN
749             mdot_n = VAVG_Flux_V_G()
750           ELSEIF(CYCLIC_Z_MF)THEN
751             mdot_n = VAVG_Flux_W_G()
752           ENDIF
753     
754           IF (mfix_isnan(mdot_n) .OR. mfix_isnan(delp_n)) THEN
755              IF (myPE.eq.PE_IO) write(*,*) mdot_n, delp_n, &
756                 ' NaN being caught in GoalSeekMassFlux '
757              RETURN
758           ENDIF
759     
760           err = abs((mdot_n - mdot_0)/mdot_0)
761           IF( err < TOL) THEN
762              MUSTIT = 0
763           ELSE
764             MUSTIT = 1
765             NIT = 1
766           ENDIF
767     
768     ! correct delp
769           if(.not.firstPass)then
770     !        delp_xyz = delp_n - omega * (delp_n - delp_nm1) * (mdot_n - mdot_0) &
771     !                          / (mdot_n - mdot_nm1)
772     ! Fail-Safe Newton's method (below) works better than the regular
773     ! Newton method (above)
774     
775              delp_xyz = delp_n - omega * (delp_n - delp_nm1) * &
776                          ((mdot_n - mdot_0)/(mdot_nm1 - mdot_0)) / &
777                          ((mdot_n - mdot_0)/(mdot_nm1 - mdot_0) - ONE)
778           else
779              firstPass=.false.
780              delp_xyz = delp_n*0.99
781           endif
782     
783           IF(MUSTIT == 0) then
784             IF(myPE.eq.PE_IO) Write(*,5500) TIME, OUTIT, delp_xyz, mdot_n
785           ENDIF
786     
787           mdot_nm1 = mdot_n
788           delp_nm1 = delp_n
789     
790           IF(CYCLIC_X_MF)THEN
791             delp_x = delp_xyz
792           ELSEIF(CYCLIC_Y_MF)THEN
793             delp_y = delp_xyz
794           ELSEIF(CYCLIC_Z_MF)THEN
795             delp_z = delp_xyz
796           ENDIF
797     
798     
799           RETURN
800     
801     5400 FORMAT(/1X,70('*')//' From: GoalSeekMassFlux',/&
802           ' Message: Number of outer iterations exceeded ', I4,/1X,70('*')/)
803     5500  Format('  Time=', G12.5, ' MassFluxIterations=', I4, ' DelP=', &
804           G12.5, ' Gas Flux=', G12.5)
805     
806           END SUBROUTINE GoalSeekMassFlux
807     
808     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
809     !                                                                      !
810     !  Module name: ADJUST_DT()                                            !
811     !  Author: M. Syamlal                                 Date: FEB-10-97  !
812     !                                                                      !
813     !  Purpose: Automatically adjust time step.                            !
814     !                                                                      !
815     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
816           LOGICAL FUNCTION ADJUSTDT()
817     
818     ! Global Variables:
819     !---------------------------------------------------------------------//
820     ! User defined type of run: new or restart
821           use run, only: RUN_TYPE
822     ! Integer flag: 0=Good, 100=initialize, otherwise bad.
823           use run, only: IER
824     ! User defined: min, max DT and adjustment factor
825           use run, only: DT_MIN, DT_MAX, DT_FAC
826     ! Flag: Use stored DT value for advancing TIME
827           use run, only: USE_DT_PREV
828     ! Flag: 2nd order time implementation
829           use run, only: CN_ON
830     ! Flag: Continue to run at DT_MIN
831           use run, only: PERSISTENT_MODE
832     ! The current number of time steps (value at restart).
833           use run, only: NSTEP, NSTEPRST
834     ! Current DT (1/DT) and direction of last change (+/-)
835           use run, only: DT, oDT, DT_DIR, STEADY_STATE
836     
837     ! Global Parameters:
838     !---------------------------------------------------------------------//
839           use param1, only: ZERO, ONE, UNDEFINED
840     
841     ! Module procedures:
842     !---------------------------------------------------------------------//
843           use error_manager
844     
845           IMPLICIT NONE
846     
847     ! Dummy Arguments:
848     
849     ! Local Variables:
850     !---------------------------------------------------------------------//
851     ! Number of steps in between DT adjustments.
852           INTEGER, PARAMETER :: STEPS_MIN = 5
853     ! Number of time steps since last DT adjustment
854           INTEGER, SAVE :: STEPS_TOT=0
855     ! number of iterations since last DT adjustment
856           INTEGER, SAVE :: NIT_TOT=0
857     ! Iterations per second for last dt
858           DOUBLE PRECISION, SAVE :: NIToS=0.0
859     ! Current number of iterations per second
860           DOUBLE PRECISION :: NITOS_NEW
861     ! Flag to half/double the current time step
862           LOGICAL :: CN_ADJUST_DT
863     !......................................................................!
864     
865     ! Initialize the function result.
866           ADJUSTDT = .FALSE.
867           USE_DT_PREV = .FALSE.
868     
869     ! Steady-state simulation.
870           IF (STEADY_STATE .OR. DT<ZERO) RETURN
871     
872     ! Local flag for adjusting the time step for CN implementation.
873           CN_ADJUST_DT = CN_ON .AND. ((RUN_TYPE=='NEW' .AND. NSTEP>1) .OR. &
874              (RUN_TYPE/='NEW' .AND. NSTEP >= (NSTEPRST+1)))
875     
876     ! Iterate successfully converged.
877     !---------------------------------------------------------------------//
878           IF(IER == 0) THEN
879     
880     ! Set back the timestep to original size which was halved previously for
881     ! 2nd order accurate time implementation.
882              IF(CN_ADJUST_DT) DT = 2.0D0*DT
883     
884     ! Calculate a new DT every STEPS_MIN time steps.
885              IF(STEPS_TOT >= STEPS_MIN) THEN
886                 NITOS_NEW = DBLE(NIT_TOT)/(STEPS_TOT*DT)
887                 IF (NITOS_NEW > NITOS) DT_DIR = DT_DIR*(-1)
888                 STEPS_TOT = 0
889                 NITOS = NITOS_NEW
890                 NIT_TOT = 0
891                 IF (DT_DIR > 0) THEN
892                    IF(NIT < MAX_NIT) DT = MIN(DT_MAX,DT/DT_FAC)
893                 ELSE
894                    DT = DT*DT_FAC
895                    IF(PERSISTENT_MODE) DT = max(DT, DT_MIN)
896                 ENDIF
897     
898     ! DT was modified. Use the stored DT should be used to update TIME.
899                 USE_DT_PREV = .TRUE.
900     
901     ! Write the convergence stats to the screen/log file.
902                 WRITE(ERR_MSG,"('DT=',g11.4,3x,'NIT/s=',A)")  &
903                    DT, trim(iVal(nint(NITOS)))
904                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., &
905                    FOOTER=.FALSE., LOG=.FALSE.)
906     
907              ELSE
908                 STEPS_TOT = STEPS_TOT + 1
909                 NIT_TOT = NIT_TOT + NIT
910              ENDIF
911     ! No need to iterate again
912              ADJUSTDT = .FALSE.
913     ! Cut the timestep into half for 2nd order accurate time implementation.
914              IF(CN_ADJUST_DT) DT = 0.5d0*DT
915     
916     ! Iterate failed to converge.
917     !---------------------------------------------------------------------//
918           ELSE
919     
920     ! Clear the error flag.
921              IER = 0
922     
923     ! Reset the timestep to original size which was halved previously for
924     ! 2nd order accurate time implementation.
925              IF(CN_ADJUST_DT) DT = 2.0d0*DT
926     
927     ! Reset counters.
928              DT_DIR = -1
929              STEPS_TOT = 0
930              NITOS = 0.
931              NIT_TOT = 0
932     
933     ! Reduce the step size.
934              DT = DT*DT_FAC
935     
936     ! The step size has decreased to the minimum.
937              IF (DT_FAC >= ONE) THEN
938     
939                 IF(PERSISTENT_MODE) THEN
940                    ADJUSTDT = .FALSE.
941                 ELSE
942                    WRITE(ERR_MSG,"(3X,A)") &
943                       'DT_FAC >= 1. Recovery not possible!'
944                    CALL FLUSH_ERR_MSG(ABORT=.TRUE., &
945                       HEADER=.FALSE., FOOTER=.FALSE.)
946                 ENDIF
947     
948              ELSEIF (DT > DT_MIN) THEN
949     
950                 WRITE(ERR_MSG,"(3X,'Recovered: Dt=',G12.5,' :-)')") DT
951                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
952     
953                 CALL RESET_NEW
954     
955     ! Iterate again with new dt
956                 ADJUSTDT = .TRUE.
957     
958     ! Cut the timestep for 2nd order accurate time implementation.
959                 IF(CN_ADJUST_DT) DT = 0.5d0*DT
960     
961     ! Set the return flag stop iterating.
962              ELSE
963     
964     ! Prevent DT from dropping below DT_MIN.
965                 IF(PERSISTENT_MODE) DT = max(DT, DT_MIN)
966                 ADJUSTDT = .FALSE.
967              ENDIF
968     
969           ENDIF
970     
971     ! Update ONE/DT variable.
972           ODT = ONE/DT
973     
974           RETURN
975           END FUNCTION ADJUSTDT
976     
977           END MODULE ITERATE
978