iterate.f

Go to the documentation of this file.

! -*- f90 -*-
      MODULE iterate

      USE exit, ONLY: mfix_exit
      USE run, ONLY: ier, tunit, time

! flag indicating convergence status with MUSTIT = 0,1,2 implying
! complete convergence, non-covergence and divergence respectively
      INTEGER :: mustit

! Number of iterations completed for current timestep
      INTEGER :: nit

! User defined maximum number of iterations
      INTEGER :: max_nit

      LOGICAL :: converged, diverged

! cpu time left
      DOUBLE PRECISION :: tleft
! Normalization factor for gas & solids pressure residual
      DOUBLE PRECISION :: normg, norms
! Set normalization factor for gas and solids pressure residual
      LOGICAL :: setg, sets
! gas & solids pressure residual
      DOUBLE PRECISION :: resg, ress
! Weight of solids in the reactor
      DOUBLE PRECISION :: smass
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: errorpercent
! Error Message
      CHARACTER(LEN=32) :: lmsg

      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ITERATE                                                 !
!  Author: M. Syamlal                                 Date: 12-APR-96  !
!                                                                      !
!  Purpose: This module controls the iterations for solving equations  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE iterate_init

      USE compar, only: mype, pe_io
      USE cutcell, only: cartesian_grid
      USE geometry, only: cyclic
      USE leqsol, only: leq_adjust
      USE output, only: full_log
      USE param1, only: one, small_number, undefined, zero
      USE run, only: dt, dt_prev, run_type, time, tstop, nstep, nsteprst, cn_on, get_tunit, steady_state
      USE time_cpu
      USE toleranc, only: norm_g, norm_s

      IMPLICIT NONE

! initializations
      dt_prev = dt
      nit = 0
      mustit = 0
      converged = .false.
      diverged  = .false.
      resg = zero
      ress = zero

      IF(norm_g == undefined) THEN
         normg = one
         setg = .false.
      ELSE
         normg = norm_g
         setg = .true.
      ENDIF

      IF(norm_s == undefined) THEN
         norms = one
         sets = .false.
      ELSE
         norms = norm_s
         sets = .true.
      ENDIF

      leq_adjust = .false.

! Initialize residuals
      CALL init_resid ()

! Initialize the routine for holding gas mass flux constant with cyclic bc
      IF(cyclic) CALL goalseekmassflux(nit, mustit, .false.)

! CPU time left
      IF (full_log) THEN
         tleft = (tstop - time)*cpuos
         CALL get_tunit (tleft, tunit)

         IF (steady_state) THEN
            CALL get_smass (smass)
            IF(mype.eq.pe_io) THEN
               WRITE (*, '(/A,G12.5, A,F9.3,1X,A)') &
                  ' Starting solids mass = ', smass
            ENDIF
         ELSE
            IF(mype.eq.pe_io) THEN
               IF ((cn_on.AND.nstep>1.AND.run_type == 'NEW') .OR. &
                   (cn_on.AND.run_type /= 'NEW' .AND.&
                    nstep >= (nsteprst+1))) THEN
                  WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)')&
                     ' Time = ', time, '  Dt = ', 2.d0*dt
               ELSE
                  WRITE (*, '(/A,G12.5, A,G12.5, A,F9.3,1X,A)') &
                     ' Time = ', time, '  Dt = ', dt
               ENDIF
            ENDIF
         ENDIF   ! if/else(steady_state)
      ENDIF   ! if(full_log)

      CALL calc_resid_mb(0, errorpercent)

! Calculate the face values of densities and mass fluxes for the first
! solve_vel_star call.
      CALL conv_rop()
      CALL calc_mflux ()
      CALL set_bc1
      CALL set_ep_factors

! JFD: modification for cartesian grid implementation
      IF(cartesian_grid) CALL cg_set_outflow

! Default/Generic Error message
      lmsg = 'Run diverged/stalled'

      RETURN
      END SUBROUTINE iterate_init


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ITERATE                                                 !
!  Author: M. Syamlal                                 Date: 12-APR-96  !
!                                                                      !
!  Purpose: This module controls the iterations for solving equations  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE do_iteration

      USE cont, only: solve_continuity
         USE cutcell, only: cartesian_grid
         USE discretelement, only: discrete_element, des_continuum_hybrid
         USE fldvar, only: ep_g, ro_g, rop_g, rop_s, p_star
         USE geometry, only: cyclic, cylindrical
         USE leqsol, only: leq_adjust
         USE mms, only: use_mms
         USE param1, only: small_number, zero, undefined, undefined_i
         USE physprop, only: mmax, ro_g0, smax
         USE pscor, only: k_cp, mcp
         USE qmom_kinetic_equation, only: qmomk
         USE residual, only: resid_p, resid
         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
         USE scalars, only: nscalar

      IMPLICIT NONE

      INTEGER :: M

      phip_out_iter=nit ! To record the output of phip
! mechanism to set the normalization factor for the correction
! after the first iteration to the corresponding residual found
! in the first iteration
      IF (.NOT.setg) THEN
         IF (resg > small_number) THEN
            normg = resg
            setg = .true.
         ENDIF
      ENDIF
      IF (.NOT.sets) THEN
         IF (ress > small_number) THEN
            norms = ress
            sets = .true.
         ENDIF
      ENDIF

! Call user-defined subroutine to set quantities that need to be updated
! every iteration
      IF (call_usr) CALL usr2

! Calculate coefficients, excluding density and reactions.
      CALL calc_coeff(ier, 1)
      IF (ier_manager()) return

! Diffusion coefficient and source terms for user-defined scalars
      IF(nscalar /= 0) CALL scalar_prop()

! Diffusion coefficient and source terms for K & Epsilon Eq.
      IF(k_epsilon) CALL k_epsilon_prop()

! Update the stress tensor trace and cross terms each subiteration
! for MMS cases.
      IF(use_mms) CALL calc_trd_and_tau()

! Solve starred velocity components
      CALL solve_vel_star(ier)

! Correct the centerline velocity for cylindrical simulations.
      IF(cylindrical) CALL radial_vel_correction

! Calculate densities.
      CALL physical_prop(ier, 0)
      IF (ier_manager()) return

! Calculate chemical reactions.
      CALL calc_rrate(ier)

! Solve solids volume fraction correction equation for close-packed
! solids phases
      IF(.NOT.(discrete_element .OR. qmomk) .OR. &
         des_continuum_hybrid) THEN
         IF (mmax > 0) THEN
! MMS:  Solve gas continuity only.
            IF(use_mms) THEN
               CALL solve_continuity(0,ier)
! Regular, non-MMS cases.
            ELSE
               IF(mmax == 1 .AND. mcp /= undefined_i)THEN
! if second phase (m=1) can overpack (e.g., bubbles) then solve its
! continuity equation
                  CALL calc_k_cp (k_cp)
                  CALL solve_epp (norms, ress, ier)
                  CALL correct_1 ()
               ELSE

! If one chooses to revert back to old mark_phase_4_cor wherein the
! continuity of the gas phase can get marked to be solved then this
! loop should start at 0.
                  DO m=1,smax ! mmax -> smax for GHD theory
! Volume fraction correction technique for one of the solids phase
! is not implemented.  This will only slow down convergence.
!                      IF (M .EQ. MCP) THEN
!                       CALL CALC_K_CP (K_CP, IER)
!                       CALL SOLVE_EPP (NORMS, RESS, IER)
!                       CALL CORRECT_1 (IER)
!                    ELSE
                        CALL solve_continuity(m,ier)
!                    ENDIF
                  ENDDO
               ENDIF   ! end if/else (mmax==1 .and. mcp /= undefined)
            ENDIF ! end if/else (MMS)

            IF(kt_type_enum == ghd_2007) CALL adjust_eps_ghd

            CALL calc_vol_fr (p_star, ro_g, rop_g, ep_g, rop_s, ier)
            IF (ier_manager()) return

         ENDIF  ! endif (mmax >0)

      ENDIF  ! end if (.not.discrete_element)


! Calculate P_star in cells where solids continuity equation is
! solved
      IF(.NOT.(discrete_element .OR. qmomk) .OR. &
         des_continuum_hybrid) THEN
         IF (mmax > 0 .AND. .NOT.friction) &
            CALL calc_p_star (ep_g, p_star)
      ENDIF

! Calculate the face values of densities.
      CALL conv_rop()

      IF (ro_g0 /= zero) THEN
! Solve fluid pressure correction equation
#ifndef FLAG_MMS
         CALL solve_pp_g (normg, resg, ier)
#endif
! Correct pressure, velocities, and density
         CALL correct_0 ()
      ENDIF

! Recalculate densities.
      CALL physical_prop(ier, 0)
      IF (ier_manager()) return

! Update wall velocities:
! modified by sof to force wall functions so even when NSW or FSW are
! declared, default wall BC will still be treated as NSW and no wall
! functions will be used
      IF(.NOT. k_epsilon) CALL set_wall_bc ()

! Calculate the face values of mass fluxes
      CALL calc_mflux ()
      CALL set_bc1
      CALL set_ep_factors

! JFD: modification for cartesian grid implementation
      IF(cartesian_grid) CALL cg_set_outflow

! Solve energy equations
      IF (energy_eq) THEN
         CALL solve_energy_eq (ier)
         IF (ier_manager()) return
      ENDIF

! Solve granular energy equation
      IF (granular_energy) THEN
         IF(.NOT.discrete_element .OR. des_continuum_hybrid) THEN
            CALL solve_granular_energy (ier)
            IF (ier_manager()) return
         ENDIF
      ENDIF

! Solve species mass balance equations.
      CALL solve_species_eq (ier)
      IF (ier_manager()) return

! Solve other scalar transport equations
      IF(nscalar /= 0) CALL solve_scalar_eq (ier)

! Solve K & Epsilon transport equations
      IF(k_epsilon) CALL solve_k_epsilon_eq (ier)

! User-defined linear equation solver parameters may be adjusted after
! the first iteration
      IF (.NOT.cyclic) leq_adjust = .true.

! Check for convergence
      CALL accum_resid ! Accumulating residuals from all the processors
      resg = resid(resid_p,0)
      ress = resid(resid_p,1)
      CALL calc_resid_mb(1, errorpercent)
      mustit = 0
      CALL check_convergence (nit, errorpercent(0), mustit)

      IF(cyclic .AND. (mustit==0 .OR. nit >= max_nit)) &
         CALL goalseekmassflux(nit, mustit, .true.)

! Display residuals
      CALL display_resid (nit)

      CALL end_iteration

      contains

      SUBROUTINE end_iteration
         IMPLICIT NONE

         ! Determine course of simulation: converge, non-converge, diverge?
         IF (mustit == 0) THEN
            IF (steady_state .AND. nit==1) RETURN   !Iterations converged
            converged = .true.
            ier = 0
         ELSEIF (mustit==2 .AND. .NOT.steady_state) THEN
            diverged = .true.
            ier = 1
         ENDIF

      END SUBROUTINE end_iteration

!----------------------------------------------------------------------!
! Function: IER_Manager                                                !
!                                                                      !
! Purpose: Identify and account for errors from called subroutines.    !
!          Returns .TRUE. for lErr >= 100, otherwise .FALSE.           !
!                                                                      !
! Reserved Error Blocks:                                               !
!                                                                      !
! [ 100,  109]: PHYSICAL_PROP                                          !
! [ 110,  119]: CALC_VOL_FR                                            !
! [ 120,  129]: SOLVE_ENERGY_EQ                                        !
! [ 130,  139]: SOLVE_SPECIES_EQ                                       !
! [ 140,  149]: SOLVE_GRANULAR_ENERGY                                  !
!                                                                      !
!----------------------------------------------------------------------!
      LOGICAL FUNCTION ier_manager()

! Default case: do nothing.
      IF(ier < 100) THEN
         ier_manager = .false.
         return
      ENDIF

! Errors with an index greater than 100 will force an exit from iterate
! and in turn, reduce the step-size, and restart the time-step.
      ier_manager = .true.
      mustit = 2

! Errors reported from PHYSICAL_PROP
!```````````````````````````````````````````````````````````````````````
      IF(ier <  110) THEN
         IF(ier ==  100) THEN
            lmsg = 'Negative gas density detected'
         ELSEIF(ier ==  101) THEN
            lmsg = 'Negative solids density detected'
         ELSE
            lmsg = 'UCE in PHYSICAL_PROP'
         ENDIF


! Errors reported from CALC_VOL_FR
!```````````````````````````````````````````````````````````````````````
      ELSEIF(ier <  120) THEN
         IF(ier ==  110) THEN
            lmsg = 'Negative void fraction detected'
         ELSE
            lmsg = 'UCE in CALC_VOL_FR'
         ENDIF


! Errors reported from SOLVE_ENERGY_EQ
!```````````````````````````````````````````````````````````````````````
      ELSEIF(ier <  130) THEN
         IF(ier ==  120) THEN
            lmsg = 'Energy Equation diverged'
         ELSE
            lmsg = 'UCE in SOLVE_ENERGY_EQ'
         ENDIF


! Errors reported from SOLVE_SPECIES_EQ
!```````````````````````````````````````````````````````````````````````
      ELSEIF(ier <  140) THEN
         IF(ier ==  130) THEN
            lmsg = 'Species Equation diverged'
         ELSE
            lmsg = 'UCE in SOLVE_SPECIES_EQ'
         ENDIF


! Errors reported from SOLVE_GRANULAR_ENERGY
!```````````````````````````````````````````````````````````````````````
      ELSEIF(ier <  150) THEN
         IF(ier ==  140) THEN
            lmsg = 'Granular Energy Eq diverged'
         ELSE
            lmsg = 'UCE in SOLVE_GRANULAR_ENERGY'
         ENDIF

! Unclassified Errors
!```````````````````````````````````````````````````````````````````````
      ELSE
         lmsg = 'Run diverged/stalled with UCE'
      ENDIF


      IF(steady_state) ier_manager = .false.

      CALL end_iteration

      return
      END FUNCTION ier_manager

      END SUBROUTINE do_iteration


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ITERATE                                                 !
!  Author: M. Syamlal                                 Date: 12-APR-96  !
!                                                                      !
!  Purpose: This module controls the iterations for solving equations  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE post_iterate

      USE compar, only: mype, pe_io
      USE funits, only: dmp_log, unit_log, unit_out
      USE machine, only: start_log, end_log
      USE run, only: dt, time

      IMPLICIT NONE

      IF (converged) THEN
         CALL log_converged
      ELSEIF (diverged) THEN
         CALL log_diverged
      ELSE
         CALL get_smass (smass)
         IF (mype.eq.pe_io) WRITE(unit_out, 5100) time, dt, nit, smass
         CALL start_log
         IF(dmp_log) WRITE(unit_log, 5100) time, dt, nit, smass
         CALL end_log
      ENDIF

      IF(.NOT.(converged .OR. diverged)) THEN
         ier = 1
      ENDIF

5100  FORMAT(1x,'t=',f11.4,' Dt=',g11.4,' NIT>',i3,' Sm= ',g12.5, &
           'MbErr%=', g11.4)

      END SUBROUTINE post_iterate

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ITERATE                                                 !
!  Author: M. Syamlal                                 Date: 12-APR-96  !
!                                                                      !
!  Purpose: This module controls the iterations for solving equations  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE log_diverged
      USE compar, only: mype, pe_io
      USE dashboard, only: f_dashboard, n_dashboard, run_status, write_dashboard
      USE funits, only: dmp_log, unit_log
      USE machine, only: start_log, end_log
      USE output, only: full_log
      USE run, only: dt, time, tstop, get_tunit
      USE time_cpu, only: cpuos

      IMPLICIT NONE

      CHARACTER(LEN=4) :: TUNIT

      IF (full_log) THEN
         CALL start_log
         CALL calc_resid_mb(1, errorpercent)

         IF(dmp_log) WRITE(unit_log,5200) time, dt, nit, &
              errorpercent(0), trim(adjustl(lmsg))
         CALL end_log

         IF (mype.EQ.pe_io) WRITE(*,5200) time, dt, nit, &
              errorpercent(0), trim(adjustl(lmsg))
      ENDIF

      ! JFD: modification for cartesian grid implementation
      IF(write_dashboard) THEN
         run_status = 'Diverged/stalled...'
         n_dashboard = n_dashboard + 1
         IF(mod(n_dashboard,f_dashboard)==0) THEN
            tleft = (tstop - time)*cpuos
            CALL get_tunit (tleft, tunit)
            CALL update_dashboard(nit,tleft,tunit)
         ENDIF
      ENDIF
5200  FORMAT(1x,'t=',f11.4,' Dt=',g11.4,' NIT=',&
                 i3,'MbErr%=', g11.4, ': ',a,' :-(')
      END SUBROUTINE log_diverged

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ITERATE                                                 !
!  Author: M. Syamlal                                 Date: 12-APR-96  !
!                                                                      !
!  Purpose: This module controls the iterations for solving equations  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE log_converged

      USE compar, only: mype, pe_io
      USE dashboard, only: f_dashboard, n_dashboard, run_status, write_dashboard
      USE error_manager, only: err_msg
      USE funits, only: dmp_log, unit_log
      USE geometry, only: cyclic_x, cyclic_y, cyclic_z
      USE geometry, only: do_i, do_j, do_k
      USE machine, only: start_log, end_log
      USE output, only: full_log, nlog
      USE physprop, only: mmax, smax
      USE run, only: dt, energy_eq, time, tstop, nstep, get_tunit
      USE time_cpu, only: cpu0, cpu_nlog, cpuos, time_nlog

      IMPLICIT NONE

      CHARACTER(LEN=4) :: TUNIT
! Perform checks and dump to screen every NLOG time steps
      IF (mod(nstep,nlog) == 0) CALL dump_to_screen

! JFD: modification for cartesian grid implementation
      IF(write_dashboard) THEN
         run_status = 'In Progress...'
         n_dashboard = n_dashboard + 1
         IF(mod(n_dashboard,f_dashboard)==0) THEN
            tleft = (tstop - time)*cpuos
            CALL get_tunit (tleft, tunit)
            CALL update_dashboard(nit,tleft,tunit)
         ENDIF
      ENDIF

      CONTAINS

      SUBROUTINE dump_to_screen
      USE error_manager, only: flush_err_msg
      IMPLICIT NONE

      ! phase index
      INTEGER :: M
! current cpu time used
      DOUBLE PRECISION :: CPU_NOW
! Heat loss from the reactor
      DOUBLE PRECISION :: HLOSS
! average velocity
      DOUBLE PRECISION :: Vavg

!-----------------------------------------------
! External functions
!-----------------------------------------------
      DOUBLE PRECISION, EXTERNAL :: VAVG_U_G, VAVG_V_G, VAVG_W_G, &
         VAVG_U_S, VAVG_V_S, VAVG_W_S

      CALL cpu_time (cpu_now)
      cpuos = (cpu_now - cpu_nlog)/(time - time_nlog)
      cpu_nlog = cpu_now
      time_nlog = time
      cpu_now = cpu_now - cpu0

      CALL calc_resid_mb(1, errorpercent)
      CALL get_smass (smass)
      IF (energy_eq) CALL get_hloss (hloss)

      CALL start_log
      IF (energy_eq) THEN
         WRITE(err_msg,5000)time, dt, nit, smass, hloss, cpu_now
         CALL flush_err_msg(header=.false., footer=.false.)
      ELSE
         WRITE(err_msg,5001) time, dt, nit, smass, cpu_now
         CALL flush_err_msg(header=.false., footer=.false.)
      ENDIF

 5000 FORMAT(1x,'t=',f11.4,' Dt=',g11.4,' NIT=',i3,' Sm=',g12.5, &
         ' Hl=',g12.5,t84,'CPU=',f8.0,' s')

 5001 FORMAT(1x,'t=',f11.4,' Dt=',g11.4,' NIT=',i3,' Sm=',g12.5, &
         t84,'CPU=',f8.0,' s')

      IF(dmp_log)WRITE (unit_log, 5002) (errorpercent(m), m=0,mmax)
      IF (full_log .and. mype.eq.pe_io) &
         WRITE (*, 5002) (errorpercent(m), m=0,mmax)

 5002 FORMAT(3x,'MbError%(0,MMAX):', 5(1x,g11.4))

      IF (.NOT.full_log) THEN
         tleft = (tstop - time)*cpuos
         CALL get_tunit (tleft, tunit)
         IF(dmp_log)WRITE (unit_log, '(46X,A,F9.3,1X,A)')
      ENDIF

      IF (cyclic_x .OR. cyclic_y .OR. cyclic_z) THEN
         IF (do_i) THEN
           vavg = vavg_u_g()
           IF(dmp_log)WRITE (unit_log, 5050) 'U_g = ', vavg
         ENDIF
         IF (do_j) THEN
           vavg = vavg_v_g()
           IF(dmp_log)WRITE (unit_log, 5050) 'V_g = ',  vavg
         ENDIF
         IF (do_k) THEN
           vavg = vavg_w_g()
           IF(dmp_log)WRITE (unit_log, 5050) 'W_g = ', vavg
         ENDIF
         DO m = 1, smax
            IF (do_i) Then
              vavg = vavg_u_s(m)
              IF(dmp_log)WRITE (unit_log, 5060) 'U_s(', m, ') = ', vavg
            ENDIF
            IF (do_j) Then
              vavg = vavg_v_s(m)
              IF(dmp_log)WRITE (unit_log, 5060) 'V_s(', m, ') = ', vavg
            ENDIF
            IF (do_k) Then
              vavg = vavg_w_s(m)
              IF(dmp_log)WRITE (unit_log, 5060) 'W_s(', m, ') = ', vavg
            ENDIF
         ENDDO
      ENDIF   ! end if cyclic_x, cyclic_y or cyclic_z

      CALL end_log

5050  FORMAT(5x,'Average ',a,g12.5)
5060  FORMAT(5x,'Average ',a,i2,a,g12.5)

      END SUBROUTINE dump_to_screen

      END SUBROUTINE log_converged

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GoalSeekMassFlux                                        !
!  Author: M. Syamlal                                 Date: 12-APR-96  !
!                                                                      !
! Purpose:  In the following subroutine the mass flux across a periodic!
! domain with pressure drop is held constant at a user-specified value.!
! This module is activated only if the user specifies a value for the  !
! keyword flux_g in the mfix.dat file.                                 !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE goalseekmassflux(NIT, MUSTIT, doit)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE bc
      USE compar
      USE constant
      USE exit, only: mfix_exit
      USE geometry
      USE param1, only: one
      USE run
      USE time_cpu
      USE utilities, ONLY: mfix_isnan
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
      INTEGER, INTENT(INOUT) :: NIT, MUSTIT
      LOGICAL, INTENT(IN) :: doit
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
      INTEGER, PARAMETER :: MAXOUTIT = 500
      DOUBLE PRECISION, PARAMETER          :: omega = 0.9
      DOUBLE PRECISION, PARAMETER          :: TOL = 1e-03
      INTEGER, SAVE :: OUTIT
      LOGICAL, SAVE :: firstPass = .true.

      DOUBLE PRECISION, SAVE  :: mdot_n, mdot_nm1, delp_n, delp_nm1, err
      DOUBLE PRECISION        :: mdot_0, delp_xyz

!-----------------------------------------------
! Functions
!-----------------------------------------------
      DOUBLE PRECISION, EXTERNAL :: VAVG_Flux_U_G, VAVG_Flux_V_G, &
                                    VAVG_Flux_W_G
!-----------------------------------------------

      IF(cyclic_x_mf)THEN
         delp_n = delp_x
      ELSEIF(cyclic_y_mf)THEN
         delp_n = delp_y
      ELSEIF(cyclic_z_mf)THEN
         delp_n = delp_z
      ELSE
         RETURN
      ENDIF

      IF(.NOT.doit) THEN
         outit = 0
         RETURN
      ENDIF

      outit = outit + 1
      IF(outit > maxoutit) THEN
         IF (mype.EQ.pe_io) write(*,5400) maxoutit
         CALL mfix_exit(0)
      ENDIF

      mdot_0 = flux_g


      ! calculate the average gas mass flux and error
      IF(cyclic_x_mf)THEN
        mdot_n = vavg_flux_u_g()
      ELSEIF(cyclic_y_mf)THEN
        mdot_n = vavg_flux_v_g()
      ELSEIF(cyclic_z_mf)THEN
        mdot_n = vavg_flux_w_g()
      ENDIF

      IF (mfix_isnan(mdot_n) .OR. mfix_isnan(delp_n)) THEN
         IF (mype.eq.pe_io) write(*,*) mdot_n, delp_n, &
            ' NaN being caught in GoalSeekMassFlux '
         RETURN
      ENDIF

      err = abs((mdot_n - mdot_0)/mdot_0)
      IF( err < tol) THEN
         mustit = 0
      ELSE
        mustit = 1
        nit = 1
      ENDIF

! correct delp
      if(.not.firstpass)then
!        delp_xyz = delp_n - omega * (delp_n - delp_nm1) * (mdot_n - mdot_0) &
!                          / (mdot_n - mdot_nm1)
! Fail-Safe Newton's method (below) works better than the regular
! Newton method (above)

         delp_xyz = delp_n - omega * (delp_n - delp_nm1) * &
                     ((mdot_n - mdot_0)/(mdot_nm1 - mdot_0)) / &
                     ((mdot_n - mdot_0)/(mdot_nm1 - mdot_0) - one)
      else
         firstpass=.false.
         delp_xyz = delp_n*0.99
      endif

      IF(mustit == 0) then
        IF(mype.eq.pe_io) Write(*,5500) time, outit, delp_xyz, mdot_n
      ENDIF

      mdot_nm1 = mdot_n
      delp_nm1 = delp_n

      IF(cyclic_x_mf)THEN
        delp_x = delp_xyz
      ELSEIF(cyclic_y_mf)THEN
        delp_y = delp_xyz
      ELSEIF(cyclic_z_mf)THEN
        delp_z = delp_xyz
      ENDIF


      RETURN

5400 FORMAT(/1x,70('*')//' From: GoalSeekMassFlux',/&
      ' Message: Number of outer iterations exceeded ', i4,/1x,70('*')/)
5500  Format('  Time=', g12.5, ' MassFluxIterations=', i4, ' DelP=', &
      g12.5, ' Gas Flux=', g12.5)

      END SUBROUTINE goalseekmassflux

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: ADJUST_DT()                                            !
!  Author: M. Syamlal                                 Date: FEB-10-97  !
!                                                                      !
!  Purpose: Automatically adjust time step.                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      LOGICAL FUNCTION adjustdt()

! Global Variables:
!---------------------------------------------------------------------//
! User defined type of run: new or restart
      use run, only: run_type
! Integer flag: 0=Good, 100=initialize, otherwise bad.
      use run, only: ier
! User defined: min, max DT and adjustment factor
      use run, only: dt_min, dt_max, dt_fac
! Flag: Use stored DT value for advancing TIME
      use run, only: use_dt_prev
! Flag: 2nd order time implementation
      use run, only: cn_on
! Flag: Continue to run at DT_MIN
      use run, only: persistent_mode
! The current number of time steps (value at restart).
      use run, only: nstep, nsteprst
! Current DT (1/DT) and direction of last change (+/-)
      use run, only: dt, odt, dt_dir, steady_state

! Global Parameters:
!---------------------------------------------------------------------//
      use param1, only: zero, one, undefined

! Module procedures:
!---------------------------------------------------------------------//
      use error_manager

      IMPLICIT NONE

! Dummy Arguments:

! Local Variables:
!---------------------------------------------------------------------//
! Number of steps in between DT adjustments.
      INTEGER, PARAMETER :: STEPS_MIN = 5
! Number of time steps since last DT adjustment
      INTEGER, SAVE :: STEPS_TOT=0
! number of iterations since last DT adjustment
      INTEGER, SAVE :: NIT_TOT=0
! Iterations per second for last dt
      DOUBLE PRECISION, SAVE :: NIToS=0.0
! Current number of iterations per second
      DOUBLE PRECISION :: NITOS_NEW
! Flag to half/double the current time step
      LOGICAL :: CN_ADJUST_DT
!......................................................................!

! Initialize the function result.
      adjustdt = .false.
      use_dt_prev = .false.

! Steady-state simulation.
      IF (steady_state .OR. dt<zero) RETURN

! Local flag for adjusting the time step for CN implementation.
      cn_adjust_dt = cn_on .AND. ((run_type=='NEW' .AND. nstep>1) .OR. &
         (run_type/='NEW' .AND. nstep >= (nsteprst+1)))

! Iterate successfully converged.
!---------------------------------------------------------------------//
      IF(ier == 0) THEN

! Set back the timestep to original size which was halved previously for
! 2nd order accurate time implementation.
         IF(cn_adjust_dt) dt = 2.0d0*dt

! Calculate a new DT every STEPS_MIN time steps.
         IF(steps_tot >= steps_min) THEN
            nitos_new = dble(nit_tot)/(steps_tot*dt)
            IF (nitos_new > nitos) dt_dir = dt_dir*(-1)
            steps_tot = 0
            nitos = nitos_new
            nit_tot = 0
            IF (dt_dir > 0) THEN
               IF(nit < max_nit) dt = min(dt_max,dt/dt_fac)
            ELSE
               dt = dt*dt_fac
               IF(persistent_mode) dt = max(dt, dt_min)
            ENDIF

! DT was modified. Use the stored DT should be used to update TIME.
            use_dt_prev = .true.

! Write the convergence stats to the screen/log file.
            WRITE(err_msg,"('DT=',g11.4,3x,'NIT/s=',A)")  &
               dt, trim(ival(nint(nitos)))
            CALL flush_err_msg(header=.false., &
               footer=.false., log=.false.)

         ELSE
            steps_tot = steps_tot + 1
            nit_tot = nit_tot + nit
         ENDIF
! No need to iterate again
         adjustdt = .false.
! Cut the timestep into half for 2nd order accurate time implementation.
         IF(cn_adjust_dt) dt = 0.5d0*dt

! Iterate failed to converge.
!---------------------------------------------------------------------//
      ELSE

! Clear the error flag.
         ier = 0

! Reset the timestep to original size which was halved previously for
! 2nd order accurate time implementation.
         IF(cn_adjust_dt) dt = 2.0d0*dt

! Reset counters.
         dt_dir = -1
         steps_tot = 0
         nitos = 0.
         nit_tot = 0

! Reduce the step size.
         dt = dt*dt_fac

! The step size has decreased to the minimum.
         IF (dt_fac >= one) THEN

            IF(persistent_mode) THEN
               adjustdt = .false.
            ELSE
               WRITE(err_msg,"(3X,A)") &
                  'DT_FAC >= 1. Recovery not possible!'
               CALL flush_err_msg(abort=.true., &
                  header=.false., footer=.false.)
            ENDIF

         ELSEIF (dt > dt_min) THEN

            WRITE(err_msg,"(3X,'Recovered: Dt=',G12.5,' :-)')") dt
            CALL flush_err_msg(header=.false., footer=.false.)

            CALL reset_new

! Iterate again with new dt
            adjustdt = .true.

! Cut the timestep for 2nd order accurate time implementation.
            IF(cn_adjust_dt) dt = 0.5d0*dt

! Set the return flag stop iterating.
         ELSE

! Prevent DT from dropping below DT_MIN.
            IF(persistent_mode) dt = max(dt, dt_min)
            adjustdt = .false.
         ENDIF

      ENDIF

! Update ONE/DT variable.
      odt = one/dt

      RETURN
      END FUNCTION adjustdt

      END MODULE iterate