File: RELATIVE:/../../../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 compar
38           USE cont
39           USE cutcell
40           USE dashboard
41           USE discretelement
42           USE fldvar
43           USE funits
44           USE geometry
45           USE indices
46           USE leqsol
47           USE machine, only: start_log, end_log
48           USE mms, only: USE_MMS
49           USE mpi_utility
50           USE output
51           USE param
52           USE param1
53           USE pgcor
54           USE physprop
55           USE pscor
56           USE qmom_kinetic_equation
57           USE residual
58           USE run
59           USE scalars
60           USE time_cpu
61           USE toleranc
62           USE visc_g
63           USE vtk
64           USE interactive, only: CHECK_INTERACT_ITER
65     
66           use error_manager
67     
68           IMPLICIT NONE
69     !-----------------------------------------------
70     ! Dummy arguments
71     !-----------------------------------------------
72     ! Error index
73           INTEGER, INTENT(INOUT) :: IER
74     ! Number of iterations
75           INTEGER, INTENT(OUT) :: NIT
76     !-----------------------------------------------
77     ! Local variables
78     !-----------------------------------------------
79     ! current cpu time used
80           DOUBLE PRECISION :: CPU_NOW
81     ! cpu time left
82           DOUBLE PRECISION :: TLEFT
83     ! flag indicating convergence status with MUSTIT = 0,1,2 implying
84     ! complete convergence, non-covergence and divergence respectively
85           INTEGER :: MUSTIT
86     ! Normalization factor for gas & solids pressure residual
87           DOUBLE PRECISION :: NORMg, NORMs
88     ! Set normalization factor for gas and solids pressure residual
89           LOGICAL :: SETg, SETs
90     ! gas & solids pressure residual
91           DOUBLE PRECISION :: RESg, RESs
92     ! Weight of solids in the reactor
93           DOUBLE PRECISION :: SMASS
94     ! Heat loss from the reactor
95           DOUBLE PRECISION :: HLOSS
96     ! phase index
97           INTEGER :: M
98     ! average velocity
99           DOUBLE PRECISION :: Vavg
100           DOUBLE PRECISION :: errorpercent(0:MMAX)
101           CHARACTER(LEN=4) :: TUNIT
102     
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     
112     !-----------------------------------------------
113     ! Include statement functions
114     !-----------------------------------------------
115     !-----------------------------------------------
116     
117     ! initializations
118           DT_prev = DT
119           NIT = 0
120           MUSTIT = 0
121           RESG = ZERO
122           RESS = ZERO
123     
124           IF (NORM_G == UNDEFINED) THEN
125              NORMG = ONE
126              SETG = .FALSE.
127           ELSE
128              NORMG = NORM_G
129              SETG = .TRUE.
130           ENDIF
131     
132           IF (NORM_S == UNDEFINED) THEN
133              NORMS = ONE
134              SETS = .FALSE.
135           ELSE
136              NORMS = NORM_S
137              SETS = .TRUE.
138           ENDIF
139     
140           LEQ_ADJUST = .FALSE.
141     
142     
143     ! Initialize residuals
144           CALL INIT_RESID ()
145     
146     
147     ! Initialize the routine for holding gas mass flux constant with cyclic bc
148           IF(CYCLIC) CALL GoalSeekMassFlux(NIT, MUSTIT, .false.)
149     
150     ! CPU time left
151           IF (FULL_LOG) THEN
152              TLEFT = (TSTOP - TIME)*CPUOS
153              CALL GET_TUNIT (TLEFT, TUNIT)
154     
155              IF (DT == UNDEFINED) THEN
156                 CALL GET_SMASS (SMASS)
157                 IF(myPE.eq.PE_IO) THEN
158                    WRITE (*, '(/A,G12.5, A,F9.3,1X,A)') &
159                       ' Starting solids mass = ', SMASS
160                 ENDIF
161              ELSE
162                 IF(myPE.eq.PE_IO) THEN
163                    IF ((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
164                        (CN_ON.AND.RUN_TYPE /= 'NEW' .AND.&
165                         NSTEP >= (NSTEPRST+1))) THEN
166                       WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)')&
167                          ' Time = ', TIME, '  Dt = ', 2.D0*DT
168                    ELSE
169                       WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)') &
170                          ' Time = ', TIME, '  Dt = ', DT
171                    ENDIF
172                 ENDIF
173              ENDIF   ! if/else(dt==undefined)
174           ENDIF   ! if(full_log)
175     
176           CALL CALC_RESID_MB(0, errorpercent)
177     
178     ! Calculate the face values of densities and mass fluxes for the first
179     ! solve_vel_star call.
180           CALL CONV_ROP()
181           CALL CALC_MFLUX ()
182           CALL SET_BC1
183           CALL SET_EP_FACTORS
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()) goto 1000
220     
221     ! Diffusion coefficient and source terms for user-defined scalars
222           IF(NScalar /= 0) CALL SCALAR_PROP()
223     
224     ! Diffusion coefficient and source terms for K & Epsilon Eq.
225           IF(K_Epsilon) CALL K_Epsilon_PROP()
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()
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()) 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)
258                       CALL SOLVE_EPP (NORMS, RESS, IER)
259                       CALL CORRECT_1 ()
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()) 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)
295           ENDIF
296     
297     ! Calculate the face values of densities.
298           CALL CONV_ROP()
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 ()
305           ENDIF
306     
307     ! Recalculate densities.
308           CALL PHYSICAL_PROP(IER, 0)
309           IF (IER_MANAGER()) 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 ()
316     
317     ! Calculate the face values of mass fluxes
318           CALL CALC_MFLUX ()
319           CALL SET_BC1
320           CALL SET_EP_FACTORS
321     
322     ! JFD: modification for cartesian grid implementation
323           IF(CARTESIAN_GRID) CALL CG_SET_OUTFLOW
324     
325     ! Solve energy equations
326           IF (ENERGY_EQ) THEN
327              CALL SOLVE_ENERGY_EQ (IER)
328              IF (IER_MANAGER()) goto 1000
329           ENDIF
330     
331     ! Solve granular energy equation
332           IF (GRANULAR_ENERGY) THEN
333              IF(.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID) THEN
334                 CALL SOLVE_GRANULAR_ENERGY (IER)
335                 IF (IER_MANAGER()) goto 1000
336              ENDIF
337           ENDIF
338     
339     ! Solve species mass balance equations.
340           CALL SOLVE_SPECIES_EQ (IER)
341           IF (IER_MANAGER()) goto 1000
342     
343     ! Solve other scalar transport equations
344           IF(NScalar /= 0) CALL SOLVE_Scalar_EQ (IER)
345     
346     ! Solve K & Epsilon transport equations
347           IF(K_Epsilon) CALL SOLVE_K_Epsilon_EQ (IER)
348     
349     
350     ! User-defined linear equation solver parameters may be adjusted after
351     ! the first iteration
352           IF (.NOT.CYCLIC) LEQ_ADJUST = .TRUE.
353     
354     
355     ! Check for convergence
356           CALL ACCUM_RESID ! Accumulating residuals from all the processors
357           RESG = RESID(RESID_P,0)
358           RESS = RESID(RESID_P,1)
359           CALL CALC_RESID_MB(1, errorpercent)
360           CALL CHECK_CONVERGENCE (NIT, errorpercent(0), MUSTIT)
361     
362           IF(CYCLIC)THEN
363             IF(MUSTIT==0 .OR. NIT >= MAX_NIT) &
364                CALL GoalSeekMassFlux(NIT, MUSTIT, .true.)
365             IF(AUTOMATIC_RESTART) RETURN
366           ENDIF
367     
368     
369     !  If not converged continue iterations; else exit subroutine.
370      1000 CONTINUE
371     !-----------------------------------------------------------------
372     
373     ! Display residuals
374           CALL DISPLAY_RESID (NIT)
375     
376     ! Determine course of simulation: converge, non-converge, diverge?
377           IF (MUSTIT == 0) THEN
378     ! ---------------------------------------------------------------->>>
379              IF (DT==UNDEFINED .AND. NIT==1) GOTO 50   !Iterations converged
380     
381     ! Perform checks and dump to screen every NLOG time steps
382              IF (MOD(NSTEP,NLOG) == 0) THEN
383                 CALL CPU_TIME (CPU_NOW)
384                 CPUOS = (CPU_NOW - CPU_NLOG)/(TIME - TIME_NLOG)
385                 CPU_NLOG = CPU_NOW
386                 TIME_NLOG = TIME
387                 CPU_NOW = CPU_NOW - CPU0
388     
389                 CALL CALC_RESID_MB(1, errorpercent)
390                 CALL GET_SMASS (SMASS)
391                 IF (ENERGY_EQ) CALL GET_HLOSS (HLOSS)
392     
393                 CALL START_LOG
394                 IF (ENERGY_EQ) THEN
395                    WRITE(ERR_MSG,5000)TIME, DT, NIT, SMASS, HLOSS, CPU_NOW
396                    CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
397                 ELSE
398                    WRITE(ERR_MSG,5001) TIME, DT, NIT, SMASS, CPU_NOW
399                    CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
400                 ENDIF
401     
402      5000 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',I3,' Sm=',G12.5, &
403              ' Hl=',G12.5,T84,'CPU=',F8.0,' s')
404     
405      5001 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',I3,' Sm=',G12.5, &
406              T84,'CPU=',F8.0,' s')
407     
408                 IF(DMP_LOG)WRITE (UNIT_LOG, 5002) (errorpercent(M), M=0,MMAX)
409                 IF (FULL_LOG .and. myPE.eq.PE_IO) &
410                    WRITE (*, 5002) (errorpercent(M), M=0,MMAX)
411     
412      5002 FORMAT(3X,'MbError%(0,MMAX):', 5(1X,G11.4))
413     
414     
415     
416                 IF (.NOT.FULL_LOG) THEN
417                    TLEFT = (TSTOP - TIME)*CPUOS
418                    CALL GET_TUNIT (TLEFT, TUNIT)
419                    IF(DMP_LOG)WRITE (UNIT_LOG, '(46X,A,F9.3,1X,A)')
420                 ENDIF
421     
422                 IF (CYCLIC_X .OR. CYCLIC_Y .OR. CYCLIC_Z) THEN
423                    IF (DO_I) THEN
424                      Vavg = VAVG_U_G()
425                      IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'U_g = ', Vavg
426                    ENDIF
427                    IF (DO_J) THEN
428                      Vavg = VAVG_V_G()
429                      IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'V_g = ',  Vavg
430                    ENDIF
431                    IF (DO_K) THEN
432                      Vavg = VAVG_W_G()
433                      IF(DMP_LOG)WRITE (UNIT_LOG, 5050) 'W_g = ', Vavg
434                    ENDIF
435                    DO M = 1, SMAX
436                       IF (DO_I) Then
437                         Vavg = VAVG_U_S(M)
438                         IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'U_s(', M, ') = ', Vavg
439                       ENDIF
440                       IF (DO_J) Then
441                         Vavg = VAVG_V_S(M)
442                         IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'V_s(', M, ') = ', Vavg
443                       ENDIF
444                       IF (DO_K) Then
445                         Vavg = VAVG_W_S(M)
446                         IF(DMP_LOG)WRITE (UNIT_LOG, 5060) 'W_s(', M, ') = ', Vavg
447                       ENDIF
448                    ENDDO
449                 ENDIF   ! end if cyclic_x, cyclic_y or cyclic_z
450     
451                 CALL END_LOG
452              ENDIF   ! end IF (MOD(NSTEP,NLOG) == 0)
453     
454     ! JFD: modification for cartesian grid implementation
455              IF(WRITE_DASHBOARD) THEN
456                 RUN_STATUS = 'In Progress...'
457                 N_DASHBOARD = N_DASHBOARD + 1
458                 IF(MOD(N_DASHBOARD,F_DASHBOARD)==0) THEN
459                    TLEFT = (TSTOP - TIME)*CPUOS
460                    CALL GET_TUNIT (TLEFT, TUNIT)
461                    CALL UPDATE_DASHBOARD(NIT,TLEFT,TUNIT)
462                 ENDIF
463              ENDIF
464     
465              IER = 0
466              RETURN   ! for if mustit =0 (converged)
467     ! end converged: go back to time_march
468     ! ----------------------------------------------------------------<<<
469     
470     ! diverged
471           ELSEIF (MUSTIT==2 .AND. DT/=UNDEFINED) THEN
472     ! ---------------------------------------------------------------->>>
473              IF (FULL_LOG) THEN
474                 CALL START_LOG
475                 CALL CALC_RESID_MB(1, errorpercent)
476     
477                 IF(DMP_LOG) WRITE(UNIT_LOG,5200) TIME, DT, NIT, &
478                    errorpercent(0), trim(adjustl(lMsg))
479                 CALL END_LOG
480     
481                 IF (myPE.EQ.PE_IO) WRITE(*,5200) TIME, DT, NIT, &
482                    errorpercent(0), trim(adjustl(lMsg))
483              ENDIF
484     
485     
486     ! JFD: modification for cartesian grid implementation
487              IF(WRITE_DASHBOARD) THEN
488                 RUN_STATUS = 'Diverged/stalled...'
489                 N_DASHBOARD = N_DASHBOARD + 1
490                 IF(MOD(N_DASHBOARD,F_DASHBOARD)==0) THEN
491                    TLEFT = (TSTOP - TIME)*CPUOS
492                    CALL GET_TUNIT (TLEFT, TUNIT)
493                    CALL UPDATE_DASHBOARD(NIT,TLEFT,TUNIT)
494                 ENDIF
495              ENDIF
496     
497              IER = 1
498              RETURN  ! for if mustit =2 (diverged)
499           ENDIF
500     ! end diverged: go back to time_march, decrease time step, try again
501     ! ----------------------------------------------------------------<<<
502     
503     ! not converged (mustit = 1, !=0,2 )
504     ! ---------------------------------------------------------------->>>
505           IF(INTERACTIVE_MODE .AND. INTERACTIVE_NITS/=UNDEFINED_I) THEN
506              CALL CHECK_INTERACT_ITER(MUSTIT)
507              IF(MUSTIT == 1) THEN
508                 GOTO 50
509              ELSE
510                 GOTO 1000
511              ENDIF
512           ELSEIF(NIT < MAX_NIT) THEN
513              MUSTIT = 0
514              GOTO 50
515           ENDIF ! continue iterate
516     ! ----------------------------------------------------------------<<<
517     
518           CALL GET_SMASS (SMASS)
519           IF (myPE.eq.PE_IO) WRITE(UNIT_OUT, 5100) TIME, DT, NIT, SMASS
520           CALL START_LOG
521           IF(DMP_LOG) WRITE(UNIT_LOG, 5100) TIME, DT, NIT, SMASS
522           CALL END_LOG
523     
524     ! SOF: MFIX will not go the next time step if MAX_NIT is reached,
525     ! instead it will decrease the time step. (IER changed from 0 to 1)
526           IER = 1
527           RETURN
528     
529      5050 FORMAT(5X,'Average ',A,G12.5)
530      5060 FORMAT(5X,'Average ',A,I2,A,G12.5)
531      5100 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT>',I3,' Sm= ',G12.5, 'MbErr%=', G11.4)
532      5200 FORMAT(1X,'t=',F11.4,' Dt=',G11.4,' NIT=',&
533           I3,'MbErr%=', G11.4, ': ',A,' :-(')
534     
535           contains
536     
537     !----------------------------------------------------------------------!
538     ! Function: IER_Manager                                                !
539     !                                                                      !
540     ! Purpose: Identify and account for errors from called subroutines.    !
541     !          Returns .TRUE. for lErr >= 100, otherwise .FALSE.           !
542     !                                                                      !
543     ! Reserved Error Blocks:                                               !
544     !                                                                      !
545     ! [ 100,  109]: PHYSICAL_PROP                                          !
546     ! [ 110,  119]: CALC_VOL_FR                                            !
547     ! [ 120,  129]: SOLVE_ENERGY_EQ                                        !
548     ! [ 130,  139]: SOLVE_SPECIES_EQ                                       !
549     ! [ 140,  149]: SOLVE_GRANULAR_ENERGY                                  !
550     !                                                                      !
551     !----------------------------------------------------------------------!
552           LOGICAL FUNCTION IER_MANAGER()
553     
554     ! Default case: do nothing.
555           IF(IER < 100) THEN
556              IER_MANAGER = .FALSE.
557              return
558           ENDIF
559     
560     ! Errors with an index greater than 100 will force an exit from iterate
561     ! and in turn, reduce the step-size, and restart the time-step.
562           IER_MANAGER = .TRUE.
563           MUSTIT = 2
564     
565     ! Errors reported from PHYSICAL_PROP
566     !```````````````````````````````````````````````````````````````````````
567           IF(IER <  110) THEN
568              IF(IER ==  100) THEN
569                 lMsg = 'Negative gas density detected'
570              ELSEIF(IER ==  101) THEN
571                 lMsg = 'Negative solids density detected'
572              ELSE
573                 lMsg = 'UCE in PHYSICAL_PROP'
574              ENDIF
575     
576     
577     ! Errors reported from CALC_VOL_FR
578     !```````````````````````````````````````````````````````````````````````
579           ELSEIF(IER <  120) THEN
580              IF(IER ==  110) THEN
581                 lMsg = 'Negative void fraction detected'
582              ELSE
583                 lMsg = 'UCE in CALC_VOL_FR'
584              ENDIF
585     
586     
587     ! Errors reported from SOLVE_ENERGY_EQ
588     !```````````````````````````````````````````````````````````````````````
589           ELSEIF(IER <  130) THEN
590              IF(IER ==  120) THEN
591                 lMsg = 'Energy Equation diverged'
592              ELSE
593                 lMsg = 'UCE in SOLVE_ENERGY_EQ'
594              ENDIF
595     
596     
597     ! Errors reported from SOLVE_SPECIES_EQ
598     !```````````````````````````````````````````````````````````````````````
599           ELSEIF(IER <  140) THEN
600              IF(IER ==  130) THEN
601                 lMsg = 'Species Equation diverged'
602              ELSE
603                 lMsg = 'UCE in SOLVE_SPECIES_EQ'
604              ENDIF
605     
606     
607     ! Errors reported from SOLVE_GRANULAR_ENERGY
608     !```````````````````````````````````````````````````````````````````````
609           ELSEIF(IER <  150) THEN
610              IF(IER ==  140) THEN
611                 lMsg = 'Granular Energy Eq diverged'
612              ELSE
613                 lMsg = 'UCE in SOLVE_GRANULAR_ENERGY'
614              ENDIF
615     
616     ! Unclassified Errors
617     !```````````````````````````````````````````````````````````````````````
618           ELSE
619              lMsg = 'Run diverged/stalled with UCE'
620           ENDIF
621     
622     
623           IF(DT == UNDEFINED) IER_MANAGER = .FALSE.
624     
625           return
626           END FUNCTION IER_MANAGER
627     
628     
629     
630           END SUBROUTINE ITERATE
631     
632     
633     
634     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
635     !  Purpose:
636     !
637     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
638           SUBROUTINE GET_TUNIT(TLEFT, TUNIT)
639     
640     !-----------------------------------------------
641     ! Modules
642     !-----------------------------------------------
643           IMPLICIT NONE
644     !-----------------------------------------------
645     ! Dummy arguments
646     !-----------------------------------------------
647           DOUBLE PRECISION, INTENT(INOUT) :: TLEFT
648           CHARACTER(LEN=4) :: TUNIT
649     !-----------------------------------------------
650     
651           IF (TLEFT < 3600.0d0) THEN
652              TUNIT = 's'
653           ELSE
654              TLEFT = TLEFT/3600.0d0
655              TUNIT = 'h'
656              IF (TLEFT >= 24.) THEN
657                 TLEFT = TLEFT/24.0d0
658                 TUNIT = 'days'
659              ENDIF
660           ENDIF
661     
662           RETURN
663           END SUBROUTINE GET_TUNIT
664     
665     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
666     !  Purpose:  In the following subroutine the mass flux across a periodic
667     !            domain with pressure drop is held constant at a
668     !            user-specified value.  This module is activated only if
669     !            the user specifies a value for the keyword flux_g in the
670     !            mfix.dat file.
671     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
672     
673           subroutine GoalSeekMassFlux(NIT, MUSTIT, doit)
674     
675     !-----------------------------------------------
676     ! Modules
677     !-----------------------------------------------
678           USE bc
679           USE geometry
680           USE constant
681           USE compar
682           USE run
683           USE time_cpu
684           USE utilities, ONLY: mfix_isnan
685           IMPLICIT NONE
686     !-----------------------------------------------
687     ! Dummy arguments
688     !-----------------------------------------------
689           INTEGER, INTENT(INOUT) :: NIT, MUSTIT
690           LOGICAL, INTENT(IN) :: doit
691     !-----------------------------------------------
692     ! Local Variables
693     !-----------------------------------------------
694           INTEGER, PARAMETER :: MAXOUTIT = 500
695           DOUBLE PRECISION, PARAMETER          :: omega = 0.9
696           DOUBLE PRECISION, PARAMETER          :: TOL = 1E-03
697           INTEGER, SAVE :: OUTIT
698           LOGICAL, SAVE :: firstPass = .true.
699     
700           DOUBLE PRECISION, SAVE  :: mdot_n, mdot_nm1, delp_n, delp_nm1, err
701           DOUBLE PRECISION        :: mdot_0, delp_xyz
702     
703     !-----------------------------------------------
704     ! Functions
705     !-----------------------------------------------
706           DOUBLE PRECISION, EXTERNAL :: VAVG_Flux_U_G, VAVG_Flux_V_G, &
707                                         VAVG_Flux_W_G
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 (mfix_isnan(mdot_n) .OR. mfix_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