time_step.f

Go to the documentation of this file.

! -*- f90 -*-
      MODULE step

      USE exit, ONLY: mfix_exit
      USE main, ONLY: nit_total, dncheck, exit_signal, ncheck
      USE run, ONLY: ier

      CONTAINS

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

      !f2py threadsafe
      USE check, only: check_mass_balance
      USE compar, only: mype
      USE dashboard, only: run_status, write_dashboard
      USE discretelement, only: des_continuum_coupled, des_continuum_hybrid, discrete_element
      USE error_manager, only: err_msg
      USE error_manager, only: flush_err_msg
      USE output, only: res_dt
      USE output_man, only: output_manager
      USE param1, only: small_number, undefined
      USE qmom_kinetic_equation, only: qmomk
      USE run, only: auto_restart, automatic_restart, call_dqmom, call_usr, chk_batchq_end
      USE run, only: cn_on, dem_solids, dt, dt_min, dt_prev, ghd_2007, kt_type_enum
      USE run, only: nstep, nsteprst, odt, pic_solids, run_type, time, tstop, units, use_dt_prev
      USE stiff_chem, only: stiff_chemistry, stiff_chem_solver
      USE toleranc, only: max_inlet_vel
      USE utilities, only: max_vel_inlet

      IMPLICIT NONE

! Terminate MFIX normally before batch queue terminates.
      IF (chk_batchq_end) CALL check_batch_queue_end(exit_signal)

      IF (call_usr) CALL usr1

! Remove solids from cells containing very small quantities of solids
      IF(.NOT.(discrete_element .OR. qmomk) .OR. &
           des_continuum_hybrid) THEN
         IF(kt_type_enum == ghd_2007) THEN
            CALL adjust_eps_ghd
         ELSE
            CALL adjust_eps
         ENDIF
      ENDIF

! sof modification: uncomment code below and modify MARK_PHASE_4_COR to
! use previous MFIX algorithm. Nov 22 2010.
! Mark the phase whose continuity will be solved and used to correct
! void/volume fraction in calc_vol_fr (see subroutine for details)
!      CALL MARK_PHASE_4_COR (PHASE_4_P_G, PHASE_4_P_S, DO_CONT, MCP,&
!           DO_P_S, SWITCH_4_P_G, SWITCH_4_P_S, IER)

! Set wall boundary conditions and transient flow b.c.'s
      CALL set_bc1
! include here so they are set before calculations rely on value
! (e.g., calc_mu_g, calc_mu_s)
      CALL set_ep_factors


      CALL output_manager(exit_signal, .false.)

! Update previous-time-step values of field variables
      CALL update_old

! Calculate coefficients
      CALL calc_coeff_all (0, ier)

! Calculate the stress tensor trace and cross terms for all phases.
      CALL calc_trd_and_tau()

! Calculate additional solids phase momentum source terms
      IF (.NOT.discrete_element .OR. des_continuum_hybrid) THEN
         CALL calc_explicit_mom_source_s
      ENDIF

! Check rates and sums of mass fractions every NLOG time steps
      IF (nstep == ncheck) THEN
         IF (dncheck < 256) dncheck = dncheck*2
         ncheck = ncheck + dncheck
! Upate the reaction rates for checking
         CALL calc_rrate(ier)
         CALL check_data_30
      ENDIF

! Double the timestep for 2nd order accurate time implementation
      IF ((cn_on.AND.nstep>1.AND.run_type == 'NEW') .OR. &
           (cn_on.AND.run_type /= 'NEW' .AND. nstep >= (nsteprst+1))) THEN
         dt = 0.5d0*dt
         odt = odt * 2.0d0
      ENDIF

! Check for maximum velocity at inlet to avoid convergence problems
      max_inlet_vel = max_vel_inlet()

      END SUBROUTINE time_step_init


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

      USE check, only: check_mass_balance
      USE compar, only: mype
      USE dashboard, only: run_status, write_dashboard
      USE error_manager, only: err_msg
      USE error_manager, only: flush_err_msg
      USE iterate, only: nit
      USE leqsol, only: solver_statistics, report_solver_stats
      USE output, only: res_dt
      USE output_man, only: output_manager
      USE param1, only: small_number, undefined
      USE qmom_kinetic_equation, only: qmomk
      USE run, only: auto_restart, automatic_restart, call_dqmom, call_usr, chk_batchq_end
      USE run, only: cn_on, dem_solids, dt, dt_min, dt_prev, ghd_2007, kt_type_enum
      USE run, only: nstep, nsteprst, odt, pic_solids, run_type, time, tstop, units, use_dt_prev, steady_state
      USE stiff_chem, only: stiff_chemistry, stiff_chem_solver
      use discretelement, only: des_explicitly_coupled
      IMPLICIT NONE

      IF(dt < dt_min) THEN

         WRITE(err_msg,1100)
         CALL flush_err_msg(header=.false., footer=.false.)
1100 FORMAT('DT < DT_MIN.  Recovery not possible!')

         IF(write_dashboard) THEN
            run_status = 'DT < DT_MIN.  Recovery not possible!'
            CALL update_dashboard(nit,0.0d0,'    ')
         ENDIF
         CALL mfix_exit(mype)
      ENDIF

! Stiff Chemistry Solver.
      IF(stiff_chemistry) THEN
         IF(des_explicitly_coupled) CALL rxns_gs_gas1
         CALL stiff_chem_solver(dt, ier)
      ENDIF


! Check over mass and elemental balances.  This routine is not active by default.
! Edit the routine and specify a reporting interval to activate it.
      CALL check_mass_balance (1)

! Other solids model implementations
      IF (dem_solids) CALL des_time_march
      IF (pic_solids) CALL pic_time_march
      IF (qmomk) CALL qmomk_time_march
      IF (call_dqmom) CALL usr_dqmom

! Advance the time step and continue
      IF((cn_on.AND.nstep>1.AND.run_type == 'NEW') .OR. &
           (cn_on.AND.run_type /= 'NEW' .AND. nstep >= (nsteprst+1))) THEN
! Double the timestep for 2nd order accurate time implementation
         dt = 2.d0*dt
         odt = odt * 0.5d0
! Perform the explicit extrapolation for CN implementation
         CALL cn_extrapol
      ENDIF

      IF (.NOT. steady_state) THEN
         IF(use_dt_prev) THEN
            time = time + dt_prev
         ELSE
            time = time + dt
         ENDIF
         use_dt_prev = .false.
         nstep = nstep + 1
      ENDIF

      nit_total = nit_total+nit

      IF(solver_statistics) CALL report_solver_stats(nit_total, nstep)

! write (*,"('Compute the Courant number')")
! call get_stats(IER)

      FLUSH (6)

      CALL output_manager(exit_signal, .true.)

      END SUBROUTINE time_step_end

      END MODULE step