File: /nfs/home/0/users/jenkins/mfix.git/model/iterate.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: ITERATE                                                 C
4     !  Purpose: This module controls the iterations for solving equations  C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 12-APR-96  C
7     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !  Revision Number: 1                                                  C
10     !  Purpose: To incorporate DES flags so that the solids                C
11     !  calculations are not done using Continuum when doing DES            C
12     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  C
13     !                                                                      C
14     !  Revision Number: 2                                                  C
15     !  Purpose: To incorporate Cartesian grid modifications                C
16     !  and utilization of the dashboard                                    C
17     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
18     !
19     !  Revision Number: 3                                                  C
20     !  Purpose: Incorporation of QMOM for the solution of the particle     C
21     !  kinetic equation                                                    C
22     !  Author: Alberto Passalacqua - Fox Research Group   Date: 02-Dec-09  C
23     !                                                                      C
24     !  Literature/Document References:                                     C
25     !                                                                      C
26     !  Variables referenced:                                               C
27     !  Variables modified:                                                 C
28     !  Local variables:                                                    C
29     !                                                                      C
30     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
31     
32           SUBROUTINE ITERATE(IER, NIT)
33     
34     !-----------------------------------------------
35     ! Modules
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     ! Dummy arguments
68     !-----------------------------------------------
69     ! Error index
70           INTEGER, INTENT(INOUT) :: IER
71     ! Number of iterations
72           INTEGER, INTENT(OUT) :: NIT
73     !-----------------------------------------------
74     ! Local variables
75     !-----------------------------------------------
76     ! current cpu time used
77           DOUBLE PRECISION :: CPU_NOW
78     ! cpu time left
79           DOUBLE PRECISION :: TLEFT
80     ! flag indicating convergence status with MUSTIT = 0,1,2 implying
81     ! complete convergence, non-covergence and divergence respectively
82           INTEGER :: MUSTIT
83     ! Normalization factor for gas & solids pressure residual
84           DOUBLE PRECISION :: NORMg, NORMs
85     ! Set normalization factor for gas and solids pressure residual
86           LOGICAL :: SETg, SETs
87     ! gas & solids pressure residual
88           DOUBLE PRECISION :: RESg, RESs
89     ! Weight of solids in the reactor
90           DOUBLE PRECISION :: SMASS
91     ! Heat loss from the reactor
92           DOUBLE PRECISION :: HLOSS
93     ! phase index
94           INTEGER :: M
95     ! average velocity
96           DOUBLE PRECISION :: Vavg
97           DOUBLE PRECISION :: errorpercent(0:MMAX)
98           LOGICAL :: ABORT_IER
99           CHARACTER(LEN=4) :: TUNIT
100     
101     ! Flag indicating which error message to print when run diverges.
102           INTEGER :: lErrMsg
103     ! Error Message
104           CHARACTER(LEN=32) :: lMsg
105     
106     !-----------------------------------------------
107     ! External functions
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     ! Include statement functions
115     !-----------------------------------------------
116     !-----------------------------------------------
117     
118     ! initializations
119           DT_prev = 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     ! Initialize residuals
145           CALL INIT_RESID (IER)
146     
147     
148     ! Initialize the routine for holding gas mass flux constant with cyclic bc
149           IF(CYCLIC) CALL GoalSeekMassFlux(NIT, MUSTIT, .false.)
150     
151     ! CPU time left
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   ! if/else(dt==undefined)
175           ENDIF   ! if(full_log)
176     
177           CALL CALC_RESID_MB(0, errorpercent)
178     
179     ! Calculate the face values of densities and mass fluxes for the first
180     ! solve_vel_star call.
181           CALL CONV_ROP(IER)
182           CALL CALC_MFLUX (IER)
183           CALL SET_BC1
184     
185     ! JFD: modification for cartesian grid implementation
186           IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
187     
188     ! Default/Generic Error message
189           lMsg = 'Run diverged/stalled'
190     
191     ! Begin iterations
192     !-----------------------------------------------------------------
193        50 CONTINUE
194           MUSTIT = 0
195           NIT = NIT + 1
196           PHIP_OUT_ITER=NIT ! To record the output of phip
197     ! mechanism to set the normalization factor for the correction
198     ! after the first iteration to the corresponding residual found
199     ! in the first iteration
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     ! Call user-defined subroutine to set quantities that need to be updated
214     ! every iteration
215           IF (CALL_USR) CALL USR2
216     
217     ! Calculate coefficients, excluding density and reactions.
218           CALL CALC_COEFF(IER, 1)
219           IF (IER_MANAGER(IER)) goto 1000
220     
221     ! Diffusion coefficient and source terms for user-defined scalars
222           IF(NScalar /= 0) CALL SCALAR_PROP(IER)
223     
224     ! Diffusion coefficient and source terms for K & Epsilon Eq.
225           IF(K_Epsilon) CALL K_Epsilon_PROP(IER)
226     
227     ! Update the stress tensor trace and cross terms each subiteration
228     ! for MMS cases.
229           IF(USE_MMS) CALL CALC_TRD_AND_TAU(IER)
230     
231     ! Solve starred velocity components
232           CALL SOLVE_VEL_STAR(IER)
233     
234     ! Correct the centerline velocity for cylindrical simulations.
235           IF(CYLINDRICAL) CALL RADIAL_VEL_CORRECTION
236     
237     ! Calculate densities.
238           CALL PHYSICAL_PROP(IER, 0)
239           IF (IER_MANAGER(IER)) goto 1000
240     
241     ! Calculate chemical reactions.
242           CALL CALC_RRATE(IER)
243     
244     ! Solve solids volume fraction correction equation for close-packed
245     ! solids phases
246           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
247              DES_CONTINUUM_HYBRID) THEN
248              IF (MMAX > 0) THEN
249     ! MMS:  Solve gas continuity only.
250                 IF(USE_MMS) THEN
251                    CALL SOLVE_CONTINUITY(0,IER)
252     ! Regular, non-MMS cases.
253                 ELSE
254                    IF(MMAX == 1 .AND. MCP /= UNDEFINED_I)THEN
255     ! if second phase (m=1) can overpack (e.g., bubbles) then solve its
256     ! continuity equation
257                       CALL CALC_K_CP (K_CP, IER)
258                       CALL SOLVE_EPP (NORMS, RESS, IER)
259                       CALL CORRECT_1 (IER)
260                    ELSE
261     
262     ! If one chooses to revert back to old mark_phase_4_cor wherein the
263     ! continuity of the gas phase can get marked to be solved then this
264     ! loop should start at 0.
265                       DO M=1,SMAX ! mmax -> smax for GHD theory
266     ! Volume fraction correction technique for one of the solids phase
267     ! is not implemented.  This will only slow down convergence.
268     !                      IF (M .EQ. MCP) THEN
269     !                       CALL CALC_K_CP (K_CP, IER)
270     !                       CALL SOLVE_EPP (NORMS, RESS, IER)
271     !                       CALL CORRECT_1 (IER)
272     !                    ELSE
273                             CALL SOLVE_CONTINUITY(M,IER)
274     !                    ENDIF
275                       ENDDO
276                    ENDIF   ! end if/else (mmax==1 .and. mcp /= undefined)
277                 ENDIF ! end if/else (MMS)
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  ! endif (mmax >0)
285     
286           ENDIF  ! end if (.not.discrete_element)
287     
288     
289     ! Calculate P_star in cells where solids continuity equation is
290     ! solved
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     ! Calculate the face values of densities.
298           CALL CONV_ROP(IER)
299     
300           IF (RO_G0 /= ZERO) THEN
301     ! Solve fluid pressure correction equation
302              CALL SOLVE_PP_G (NORMG, RESG, IER)
303     ! Correct pressure, velocities, and density
304              CALL CORRECT_0 (IER)
305           ENDIF
306     
307     ! Recalculate densities.
308           CALL PHYSICAL_PROP(IER, 0)
309           IF (IER_MANAGER(IER)) goto 1000
310     
311     ! Update wall velocities:
312     ! modified by sof to force wall functions so even when NSW or FSW are
313     ! declared, default wall BC will still be treated as NSW and no wall
314     ! functions will be used
315           IF(.NOT. K_EPSILON) CALL SET_WALL_BC (IER)
316     
317     ! Calculate the face values of mass fluxes
318           CALL CALC_MFLUX (IER)
319           CALL SET_BC1
320     
321     ! JFD: modification for cartesian grid implementation
322           IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
323     
324     ! Solve energy equations
325           IF (ENERGY_EQ) THEN
326              CALL SOLVE_ENERGY_EQ (IER)
327              IF (IER_MANAGER(IER)) goto 1000
328           ENDIF
329     
330     ! Solve granular energy equation
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     ! Solve species mass balance equations.
339           CALL SOLVE_SPECIES_EQ (IER)
340           IF (IER_MANAGER(IER)) goto 1000
341     
342     ! Solve other scalar transport equations
343           IF(NScalar /= 0) CALL SOLVE_Scalar_EQ (IER)
344     
345     ! Solve K & Epsilon transport equations
346           IF(K_Epsilon) CALL SOLVE_K_Epsilon_EQ (IER)
347     
348     
349     ! User-defined linear equation solver parameters may be adjusted after
350     ! the first iteration
351           IF (.NOT.CYCLIC) LEQ_ADJUST = .TRUE.
352     
353     
354     ! Check for convergence
355           CALL ACCUM_RESID ! Accumulating residuals from all the processors
356           RESG = 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     !  If not converged continue iterations; else exit subroutine.
369      1000 CONTINUE
370     !-----------------------------------------------------------------
371     
372     ! Display residuals
373           IF (FULL_LOG) CALL DISPLAY_RESID (NIT, IER)
374     
375     ! Determine course of simulation: converge, non-converge, diverge?
376           IF (MUSTIT == 0) THEN
377     ! ---------------------------------------------------------------->>>
378              IF (DT==UNDEFINED .AND. NIT==1) GOTO 50   !Iterations converged
379     
380     ! Perform checks and dump to screen every NLOG time steps
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   ! end if cyclic_x, cyclic_y or cyclic_z
442     
443                 CALL END_LOG
444              ENDIF   ! end IF (MOD(NSTEP,NLOG) == 0)
445     
446     ! JFD: modification for cartesian grid implementation
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   ! for if mustit =0 (converged)
459     ! end converged: go back to time_march
460     ! ----------------------------------------------------------------<<<
461     
462     ! diverged
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     ! JFD: modification for cartesian grid implementation
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  ! for if mustit =2 (diverged)
491           ENDIF
492     ! end diverged: go back to time_march, decrease time step, try again
493     ! ----------------------------------------------------------------<<<
494     
495     ! not converged (mustit = 1, !=0,2 )
496     ! ---------------------------------------------------------------->>>
497           IF (NIT < MAX_NIT) THEN
498              MUSTIT = 0
499              GOTO 50
500           ENDIF ! continue iterate
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     ! SOF: MFIX will not go the next time step if MAX_NIT is reached,
513     ! instead it will decrease the time step. (IER changed from 0 to 1)
514           IER = 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     ! Function: IER_Manager                                                !
534     !                                                                      !
535     ! Purpose: Identify and account for errors from called subroutines.    !
536     !          Returns .TRUE. for lErr >= 100, otherwise .FALSE.           !
537     !                                                                      !
538     ! Reserved Error Blocks:                                               !
539     !                                                                      !
540     ! [ 100,  109]: PHYSICAL_PROP                                          !
541     ! [ 110,  119]: CALC_VOL_FR                                            !
542     ! [ 120,  129]: SOLVE_ENERGY_EQ                                        !
543     ! [ 130,  139]: SOLVE_SPECIES_EQ                                       !
544     ! [ 140,  149]: SOLVE_GRANULAR_ENERGY                                  !
545     !                                                                      !
546     !----------------------------------------------------------------------!
547           LOGICAL FUNCTION IER_MANAGER(lErr)
548     
549           INTEGER, intent(in) :: lErr
550     
551     ! Default case: do nothing.
552           IF(IER < 100) THEN
553              IER_MANAGER = .FALSE.
554              return
555           ENDIF
556     
557     ! Errors with an index greater than 100 will force an exit from iterate
558     ! and in turn, reduce the step-size, and restart the time-step.
559           IER_MANAGER = .TRUE.
560           MUSTIT = 2
561     
562     ! Errors reported from PHYSICAL_PROP
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     ! Errors reported from CALC_VOL_FR
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     ! Errors reported from SOLVE_ENERGY_EQ
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     ! Errors reported from SOLVE_SPECIES_EQ
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     ! Errors reported from SOLVE_GRANULAR_ENERGY
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     ! Unclassified Errors
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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
632     !  Purpose:
633     !
634     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
635           SUBROUTINE GET_TUNIT(TLEFT, TUNIT)
636     
637     !-----------------------------------------------
638     ! Modules
639     !-----------------------------------------------
640           IMPLICIT NONE
641     !-----------------------------------------------
642     ! Dummy arguments
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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
665     !  Purpose:  In the following subroutine the mass flux across a periodic
666     !            domain with pressure drop is held constant at a
667     !            user-specified value.  This module is activated only if
668     !            the user specifies a value for the keyword flux_g in the
669     !            mfix.dat file.
670     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
671     
672           subroutine GoalSeekMassFlux(NIT, MUSTIT, doit)
673     
674     !-----------------------------------------------
675     ! Modules
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     ! Dummy arguments
686     !-----------------------------------------------
687           INTEGER, INTENT(INOUT) :: NIT, MUSTIT
688           LOGICAL, INTENT(IN) :: doit
689     !-----------------------------------------------
690     ! Local Variables
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     ! Functions
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           ! calculate the average gas mass flux and error
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     ! correct delp
759           if(.not.firstPass)then
760     !        delp_xyz = delp_n - omega * (delp_n - delp_nm1) * (mdot_n - mdot_0) &
761     !                          / (mdot_n - mdot_nm1)
762     ! Fail-Safe Newton's method (below) works better than the regular
763     ! Newton method (above)
764     
765              delp_xyz = 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