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

1     ! -*- f90 -*-
2           MODULE STEP
3     
4           USE EXIT, ONLY: MFIX_EXIT
5           USE MAIN, ONLY: NIT_TOTAL, DNCHECK, EXIT_SIGNAL, NCHECK
6           USE RUN, ONLY: IER
7     
8           CONTAINS
9     
10     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
11     !                                                                      !
12     !  Subroutine: TIME_STEP_INIT                                          !
13     !  Author: M. Syamlal                                 Date: 12-APR-96  !
14     !                                                                      !
15     !  Purpose: This module controls the iterations for solving equations  !
16     !                                                                      !
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
18           SUBROUTINE TIME_STEP_INIT
19     
20           !f2py threadsafe
21           USE check, only: check_mass_balance
22           USE compar, only: mype
23           USE dashboard, only: run_status, write_dashboard
24           USE discretelement, only: des_continuum_coupled, des_continuum_hybrid, discrete_element
25           USE error_manager, only: err_msg
26           USE error_manager, only: flush_err_msg
27           USE output, only: res_dt
28           USE output_man, only: output_manager
29           USE param1, only: small_number, undefined
30           USE qmom_kinetic_equation, only: qmomk
31           USE run, only: auto_restart, automatic_restart, call_dqmom, call_usr, chk_batchq_end
32           USE run, only: cn_on, dem_solids, dt, dt_min, dt_prev, ghd_2007, kt_type_enum
33           USE run, only: nstep, nsteprst, odt, pic_solids, run_type, time, tstop, units, use_dt_prev
34           USE stiff_chem, only: stiff_chemistry, stiff_chem_solver
35           USE toleranc, only: max_inlet_vel
36           USE utilities, only: max_vel_inlet
37     
38           IMPLICIT NONE
39     
40     ! Terminate MFIX normally before batch queue terminates.
41           IF (CHK_BATCHQ_END) CALL CHECK_BATCH_QUEUE_END(EXIT_SIGNAL)
42     
43           IF (CALL_USR) CALL USR1
44     
45     ! Remove solids from cells containing very small quantities of solids
46           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
47                DES_CONTINUUM_HYBRID) THEN
48              IF(KT_TYPE_ENUM == GHD_2007) THEN
49                 CALL ADJUST_EPS_GHD
50              ELSE
51                 CALL ADJUST_EPS
52              ENDIF
53           ENDIF
54     
55     ! sof modification: uncomment code below and modify MARK_PHASE_4_COR to
56     ! use previous MFIX algorithm. Nov 22 2010.
57     ! Mark the phase whose continuity will be solved and used to correct
58     ! void/volume fraction in calc_vol_fr (see subroutine for details)
59     !      CALL MARK_PHASE_4_COR (PHASE_4_P_G, PHASE_4_P_S, DO_CONT, MCP,&
60     !           DO_P_S, SWITCH_4_P_G, SWITCH_4_P_S, IER)
61     
62     ! Set wall boundary conditions and transient flow b.c.'s
63           CALL SET_BC1
64     ! include here so they are set before calculations rely on value
65     ! (e.g., calc_mu_g, calc_mu_s)
66           CALL SET_EP_FACTORS
67     
68     
69           CALL OUTPUT_MANAGER(EXIT_SIGNAL, .FALSE.)
70     
71     ! Update previous-time-step values of field variables
72           CALL UPDATE_OLD
73     
74     ! Calculate coefficients
75           CALL CALC_COEFF_ALL (0, IER)
76     
77     ! Calculate the stress tensor trace and cross terms for all phases.
78           CALL CALC_TRD_AND_TAU()
79     
80     ! Calculate additional solids phase momentum source terms
81           IF (.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID) THEN
82              CALL CALC_EXPLICIT_MOM_SOURCE_S
83           ENDIF
84     
85     ! Check rates and sums of mass fractions every NLOG time steps
86           IF (NSTEP == NCHECK) THEN
87              IF (DNCHECK < 256) DNCHECK = DNCHECK*2
88              NCHECK = NCHECK + DNCHECK
89     ! Upate the reaction rates for checking
90              CALL CALC_RRATE(IER)
91              CALL CHECK_DATA_30
92           ENDIF
93     
94     ! Double the timestep for 2nd order accurate time implementation
95           IF ((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
96                (CN_ON.AND.RUN_TYPE /= 'NEW' .AND. NSTEP >= (NSTEPRST+1))) THEN
97              DT = 0.5d0*DT
98              ODT = ODT * 2.0d0
99           ENDIF
100     
101     ! Check for maximum velocity at inlet to avoid convergence problems
102           MAX_INLET_VEL = MAX_VEL_INLET()
103     
104           END SUBROUTINE TIME_STEP_INIT
105     
106     
107     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
108     !                                                                      !
109     !  Subroutine: TIME_STEP_INIT                                          !
110     !  Author: M. Syamlal                                 Date: 12-APR-96  !
111     !                                                                      !
112     !  Purpose: This module controls the iterations for solving equations  !
113     !                                                                      !
114     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
115           SUBROUTINE TIME_STEP_END
116     
117           USE check, only: check_mass_balance
118           USE compar, only: mype
119           USE dashboard, only: run_status, write_dashboard
120           USE error_manager, only: err_msg
121           USE error_manager, only: flush_err_msg
122           USE iterate, only: nit
123           USE leqsol, only: solver_statistics, report_solver_stats
124           USE output, only: res_dt
125           USE output_man, only: output_manager
126           USE param1, only: small_number, undefined
127           USE qmom_kinetic_equation, only: qmomk
128           USE run, only: auto_restart, automatic_restart, call_dqmom, call_usr, chk_batchq_end
129           USE run, only: cn_on, dem_solids, dt, dt_min, dt_prev, ghd_2007, kt_type_enum
130           USE run, only: nstep, nsteprst, odt, pic_solids, run_type, time, tstop, units, use_dt_prev, steady_state
131           USE stiff_chem, only: stiff_chemistry, stiff_chem_solver
132           use discretelement, only: DES_EXPLICITLY_COUPLED
133           IMPLICIT NONE
134     
135           IF(DT < DT_MIN) THEN
136     
137              WRITE(ERR_MSG,1100)
138              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
139     1100 FORMAT('DT < DT_MIN.  Recovery not possible!')
140     
141              IF(WRITE_DASHBOARD) THEN
142                 RUN_STATUS = 'DT < DT_MIN.  Recovery not possible!'
143                 CALL UPDATE_DASHBOARD(NIT,0.0d0,'    ')
144              ENDIF
145              CALL MFIX_EXIT(MyPE)
146           ENDIF
147     
148     ! Stiff Chemistry Solver.
149           IF(STIFF_CHEMISTRY) THEN
150              IF(DES_EXPLICITLY_COUPLED) CALL RXNS_GS_GAS1
151              CALL STIFF_CHEM_SOLVER(DT, IER)
152           ENDIF
153     
154     
155     ! Check over mass and elemental balances.  This routine is not active by default.
156     ! Edit the routine and specify a reporting interval to activate it.
157           CALL CHECK_MASS_BALANCE (1)
158     
159     ! Other solids model implementations
160           IF (DEM_SOLIDS) CALL DES_TIME_MARCH
161           IF (PIC_SOLIDS) CALL PIC_TIME_MARCH
162           IF (QMOMK) CALL QMOMK_TIME_MARCH
163           IF (CALL_DQMOM) CALL USR_DQMOM
164     
165     ! Advance the time step and continue
166           IF((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
167                (CN_ON.AND.RUN_TYPE /= 'NEW' .AND. NSTEP >= (NSTEPRST+1))) THEN
168     ! Double the timestep for 2nd order accurate time implementation
169              DT = 2.d0*DT
170              ODT = ODT * 0.5d0
171     ! Perform the explicit extrapolation for CN implementation
172              CALL CN_EXTRAPOL
173           ENDIF
174     
175           IF (.NOT. STEADY_STATE) THEN
176              IF(USE_DT_PREV) THEN
177                 TIME = TIME + DT_PREV
178              ELSE
179                 TIME = TIME + DT
180              ENDIF
181              USE_DT_PREV = .FALSE.
182              NSTEP = NSTEP + 1
183           ENDIF
184     
185           NIT_TOTAL = NIT_TOTAL+NIT
186     
187           IF(SOLVER_STATISTICS) CALL REPORT_SOLVER_STATS(NIT_TOTAL, NSTEP)
188     
189     ! write (*,"('Compute the Courant number')")
190     ! call get_stats(IER)
191     
192           FLUSH (6)
193     
194           CALL OUTPUT_MANAGER(EXIT_SIGNAL, .TRUE.)
195     
196           END SUBROUTINE TIME_STEP_END
197     
198           END MODULE STEP
199