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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  SUBROUTINE: TIME_MARCH                                              !
4     !  Author: M. Syamlal                                 Date: 21-JAN-92  !
5     !                                                                      !
6     !  Purpose: Controlling module for time marching and finding the       !
7     !           solution of equations from TIME to TSTOP at intervals of   !
8     !           DT, updating the b.c.'s, and creating output.              !
9     !                                                                      !
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
11     
12           SUBROUTINE TIME_MARCH
13     
14     !-----------------------------------------------
15     ! Modules
16     !-----------------------------------------------
17           USE param
18           USE param1
19           USE run
20           USE physprop
21           USE fldvar
22           USE geometry
23           USE pgcor
24           USE pscor
25           USE cont
26           USE tau_g
27           USE tau_s
28           USE visc_g
29           USE visc_s
30           USE funits
31           USE vshear
32           USE scalars
33           USE toleranc
34           USE drag
35           USE rxns, only: nRR
36           USE compar
37           USE time_cpu
38           USE discretelement
39           USE leqsol, only: SOLVER_STATISTICS, REPORT_SOLVER_STATS
40           use mpi_utility
41           USE cdist
42           USE MFIX_netcdf
43           USE cutcell
44           USE vtk
45           USE qmom_kinetic_equation
46           USE dashboard
47           USE indices
48           USE bc
49           USE coeff
50           USE stiff_chem, only : STIFF_CHEMISTRY, STIFF_CHEM_SOLVER
51           USE vtp
52     
53           use output, only: RES_DT, NLOG
54     
55     
56           IMPLICIT NONE
57     !-----------------------------------------------
58     ! Local variables
59     !-----------------------------------------------
60     ! Flag to indicate one pass through iterate for steady
61     ! state conditions.
62           LOGICAL :: FINISH
63     
64     ! Loop indices
65           INTEGER :: L, M , I, IJK, N
66     ! Error index
67           INTEGER :: IER
68     ! Number of iterations
69           INTEGER :: NIT, NIT_TOTAL
70     ! used for activating check_data_30
71           INTEGER :: NCHECK, DNCHECK
72     ! dummy logical variable for initializing adjust_dt
73           LOGICAL :: dummy
74     
75     ! AEOLUS : stop trigger mechanism to terminate MFIX normally before
76     ! batch queue terminates
77           DOUBLE PRECISION :: CPU_STOP
78           LOGICAL :: eofBATCHQ
79     
80     !-----------------------------------------------
81     ! External functions
82     !-----------------------------------------------
83     ! use function MAX_VEL_INLET to compute max. velocity at inlet
84           DOUBLE PRECISION :: MAX_VEL_INLET
85     ! use function vavg_v_g to catch NaN's
86           DOUBLE PRECISION :: VAVG_U_G, VAVG_V_G, VAVG_W_G, X_vavg
87     
88           LOGICAL , EXTERNAL :: ADJUST_DT
89     !-----------------------------------------------
90     
91           IF(AUTOMATIC_RESTART) RETURN
92     
93     
94           FINISH  = .FALSE.
95           NCHECK  = NSTEP
96           DNCHECK = 1
97           CPU_IO  = ZERO
98           NIT_TOTAL = 0
99           eofBATCHQ = .FALSE.
100     
101     
102           CALL INIT_OUTPUT_VARS
103     
104     ! Parse residual strings
105           CALL PARSE_RESID_STRING (IER)
106     
107     ! Call user-defined subroutine to set constants, check data, etc.
108           IF (CALL_USR) CALL USR0
109     
110           CALL RRATES_INIT(IER)
111     
112     ! Calculate all the coefficients once before entering the time loop
113           CALL INIT_COEFF(IER)
114     
115           DO M=1, MMAX
116              CALL ZERO_ARRAY (F_gs(1,M), IER)
117           ENDDO
118     
119     ! Remove undefined values at wall cells for scalars
120           CALL UNDEF_2_0 (ROP_G, IER)
121           DO M = 1, MMAX
122              CALL UNDEF_2_0 (ROP_S(1,M), IER)
123           ENDDO
124     
125     ! Initialize d's and e's to zero
126           DO M = 0, MMAX
127              CALL ZERO_ARRAY (D_E(1,M), IER)
128              CALL ZERO_ARRAY (D_N(1,M), IER)
129              CALL ZERO_ARRAY (D_T(1,M), IER)
130           ENDDO
131           CALL ZERO_ARRAY (E_E, IER)
132           CALL ZERO_ARRAY (E_N, IER)
133           CALL ZERO_ARRAY (E_T, IER)
134     
135     ! Initialize adjust_ur
136           dummy = ADJUST_DT(100, 0)
137     
138     ! calculate shear velocities if periodic shear BCs are used
139           IF(SHEAR) CALL CAL_D(V_sh)
140     
141     ! Initialize check_mass_balance.  This routine is not active by default.
142     ! Specify a reporting interval (hard-wired in the routine) to activate
143     ! the routine.
144           Call check_mass_balance (0)
145     
146     ! sof modification: now it's only needed to do this once before time-loop
147     ! Mark the phase whose continuity will be solved and used to correct
148     ! void/volume fraction in calc_vol_fr (see subroutine for details)
149           CALL MARK_PHASE_4_COR (PHASE_4_P_G, PHASE_4_P_S, DO_CONT, MCP,&
150               DO_P_S, SWITCH_4_P_G, SWITCH_4_P_S, IER)
151     
152     ! uncoupled discrete element simulations do not need to be within
153     ! the two fluid model time-loop
154           IF(DISCRETE_ELEMENT.AND.(.NOT.DES_CONTINUUM_COUPLED))  THEN
155              IF(WRITE_VTK_FILES) THEN
156                 DO L = 1, DIMENSION_VTK
157                    ! CALL WRITE_VTP_FILE(L)
158                 ENDDO
159              ENDIF
160              CALL DES_TIME_MARCH
161              CALL CPU_TIME(CPU_STOP)
162              CPU_STOP = CPU_STOP - CPU00
163              IF(myPE.EQ.PE_IO) &
164                 write(*,"('Elapsed CPU time = ',E15.6,' sec')") CPU_STOP
165              CALL PARALLEL_FIN
166              STOP
167           ENDIF
168     
169     
170     ! The TIME loop begins here.............................................
171      100  CONTINUE
172     
173     ! AEOLUS: stop trigger mechanism to terminate MFIX normally before batch
174     ! queue terminates
175           IF (CHK_BATCHQ_END) CALL CHECK_BATCH_QUEUE_END
176     
177           IF (CALL_USR) CALL USR1
178     
179     ! Remove solids from cells containing very small quantities of solids
180           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
181              DES_CONTINUUM_HYBRID) THEN
182              IF(KT_TYPE_ENUM == GHD_2007) THEN
183                 CALL ADJUST_EPS_GHD
184              ELSE
185                 CALL ADJUST_EPS
186              ENDIF
187           ENDIF
188     
189     ! sof modification: uncomment code below and modify MARK_PHASE_4_COR to
190     ! use previous MFIX algorithm. Nov 22 2010.
191     ! Mark the phase whose continuity will be solved and used to correct
192     ! void/volume fraction in calc_vol_fr (see subroutine for details)
193     !      CALL MARK_PHASE_4_COR (PHASE_4_P_G, PHASE_4_P_S, DO_CONT, MCP,&
194     !           DO_P_S, SWITCH_4_P_G, SWITCH_4_P_S, IER)
195     
196     ! Set wall boundary conditions and transient flow b.c.'s
197           CALL SET_BC1
198     
199           CALL OUTPUT_MANAGER(eofBATCHQ, FINISH)
200     
201           IF (DT == UNDEFINED) THEN
202              IF (FINISH) THEN
203                 RETURN
204              ELSE
205                 FINISH = .TRUE.
206              ENDIF
207     
208     ! Mechanism to terminate MFIX normally before batch queue terminates.
209           ELSEIF (TIME + 0.1d0*DT >= TSTOP .OR. eofBATCHQ) THEN
210              IF(SOLVER_STATISTICS) &
211                 CALL REPORT_SOLVER_STATS(NIT_TOTAL, NSTEP)
212              RETURN
213           ENDIF
214     
215     ! Update previous-time-step values of field variables
216           CALL UPDATE_OLD
217     
218     ! Calculate coefficients
219           CALL CALC_COEFF_ALL (0, IER)
220     
221     ! Calculate the stress tensor trace and cross terms for all phases.
222           CALL CALC_TRD_AND_TAU(IER)
223     
224     ! Calculate additional solid phase momentum source terms
225     ! that arise from kinetic theory constitutive relations
226           IF (.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID) THEN
227              CALL CALC_KTMOMSOURCE_U_S (IER)
228              CALL CALC_KTMOMSOURCE_V_S (IER)
229              CALL CALC_KTMOMSOURCE_W_S (IER)
230           ENDIF
231     
232     ! Check rates and sums of mass fractions every NLOG time steps
233           IF (NSTEP == NCHECK) THEN
234              IF (DNCHECK < 256) DNCHECK = DNCHECK*2
235              NCHECK = NCHECK + DNCHECK
236     ! Upate the reaction rates for checking
237              CALL CALC_RRATE(IER)
238              CALL CHECK_DATA_30
239           ENDIF
240     
241     ! Double the timestep for 2nd order accurate time implementation
242           IF ((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
243               (CN_ON.AND.RUN_TYPE /= 'NEW' .AND. NSTEP >= (NSTEPRST+1))) THEN
244              DT = 0.5d0*DT
245              ODT = ODT * 2.0d0
246           ENDIF
247     
248     ! Check for maximum velocity at inlet to avoid convergence problems
249           MAX_INLET_VEL = 100.0d0*MAX_VEL_INLET()
250     ! if no inlet velocity is specified, use an upper limit defined in
251     ! toleranc_mod.f
252           IF(MAX_INLET_VEL <= SMALL_NUMBER) THEN
253              MAX_INLET_VEL = MAX_ALLOWED_VEL
254              IF (UNITS == 'SI') MAX_INLET_VEL = 1D-2 * MAX_ALLOWED_VEL
255           ENDIF
256     ! Scale the value using a user defined scale factor
257           MAX_INLET_VEL = MAX_INLET_VEL * MAX_INLET_VEL_FAC
258     
259     ! Advance the solution in time by iteratively solving the equations
260      150  CALL ITERATE (IER, NIT)
261           IF(AUTOMATIC_RESTART) RETURN
262     
263     ! Just to Check for NaN's, Uncomment the following lines and also lines
264     ! of code in  VAVG_U_G, VAVG_V_G, VAVG_W_G to use.
265     !      X_vavg = VAVG_U_G ()
266     !      X_vavg = VAVG_V_G ()
267     !      X_vavg = VAVG_W_G ()
268     !      IF(AUTOMATIC_RESTART) RETURN
269     
270           DO WHILE (ADJUST_DT(IER,NIT))
271              CALL ITERATE (IER, NIT)
272           ENDDO
273     
274     
275           IF(DT.LT.DT_MIN) THEN
276              IF(TIME.LE.RES_DT .AND. AUTO_RESTART) THEN
277                 WRITE(ERR_MSG,1000)
278                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
279                 AUTO_RESTART = .FALSE.
280              ENDIF
281     
282      1000 FORMAT('Automatic restart not possible as Total Time < RES_DT')
283     
284              IF(AUTO_RESTART) THEN
285                 AUTOMATIC_RESTART = .TRUE.
286                 RETURN
287              ELSE
288                 IF(WRITE_DASHBOARD) THEN
289                    RUN_STATUS = 'DT < DT_MIN.  Recovery not possible!'
290                    CALL UPDATE_DASHBOARD(NIT,0.0d0,'    ')
291                 ENDIF
292                 CALL MFIX_EXIT(MyPE)
293              ENDIF
294           ENDIF
295     
296     ! Stiff Chemistry Solver.
297           IF(STIFF_CHEMISTRY) THEN
298              CALL STIFF_CHEM_SOLVER(DT, IER)
299              IF(IER /= 0) THEN
300                 dummy = ADJUST_DT(IER, NIT)
301                 GOTO 150
302              ENDIF
303           ENDIF
304     
305     ! Check over mass and elemental balances.  This routine is not active by default.
306     ! Edit the routine and specify a reporting interval to activate it.
307           CALL CHECK_MASS_BALANCE (1)
308     
309     ! Othe solids model implementations
310           IF (DEM_SOLIDS) CALL DES_TIME_MARCH
311           IF (PIC_SOLIDS) CALL PIC_TIME_MARCH
312           IF (QMOMK) CALL QMOMK_TIME_MARCH
313           IF (CALL_DQMOM) CALL USR_DQMOM
314     
315     ! Advance the time step and continue
316           IF((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
317              (CN_ON.AND.RUN_TYPE /= 'NEW' .AND. NSTEP >= (NSTEPRST+1))) THEN
318     ! Double the timestep for 2nd order accurate time implementation
319              DT = 2.d0*DT
320              ODT = ODT * 0.5d0
321     ! Perform the explicit extrapolation for CN implementation
322              CALL CN_EXTRAPOL
323           ENDIF
324     
325           IF (DT /= UNDEFINED) THEN
326              IF(use_DT_prev) THEN
327                 TIME = TIME + DT_PREV
328              ELSE
329                 TIME = TIME + DT
330              ENDIF
331              use_DT_prev = .FALSE.
332              NSTEP = NSTEP + 1
333           ENDIF
334     
335           NIT_TOTAL = NIT_TOTAL+NIT
336     
337     ! write (*,"('Compute the Courant number')")
338     ! call get_stats(IER)
339     
340           FLUSH (6)
341     ! The TIME loop ends here....................................................
342           GOTO 100
343     
344           IF(SOLVER_STATISTICS) CALL REPORT_SOLVER_STATS(NIT_TOTAL, NSTEP)
345     
346           contains
347     
348     !----------------------------------------------------------------------!
349     !                                                                      !
350     !  Subroutine: CHECK_BATCH_QUEUE_END                                   !
351     !  Author: A.Gel                                      Date:            !
352     !                                                                      !
353     !  Purpose:                                                            !
354     !                                                                      !
355     !----------------------------------------------------------------------!
356           SUBROUTINE CHECK_BATCH_QUEUE_END
357     
358           use time_cpu, only: WALL_START
359     
360           use error_manager
361     
362     ! Logical flags for hault cases.
363           LOGICAL :: USER_HAULT, WALL_HAULT
364     ! Elapsed wall time, and fancy formatted buffer/batch queue times.
365           DOUBLE PRECISION :: WALL_STOP, FANCY_BUFF, FANCY_BATCH
366     ! Time units for formatted output.
367           CHARACTER(LEN=4) :: WT_UNIT, BF_UNIT, BC_UNIT
368     ! External function
369           DOUBLE PRECISION :: WALL_TIME
370     
371     ! Calculate the current elapsed wall time.
372           WALL_STOP = WALL_TIME()
373           WALL_STOP = WALL_STOP - WALL_START
374     
375     ! Set flags for wall time exceeded or user specified hault.
376           WALL_HAULT = ((WALL_STOP+TERM_BUFFER) >= BATCH_WALLCLOCK)
377           INQUIRE(file="MFIX.STOP", exist=USER_HAULT)
378     
379     ! Report that the max user wall time was reached and exit.
380           IF(WALL_HAULT) THEN
381              CALL GET_TUNIT(WALL_STOP,WT_UNIT)
382              FANCY_BUFF = TERM_BUFFER
383              CALL GET_TUNIT(FANCY_BUFF, BF_UNIT)
384              FANCY_BATCH = BATCH_WALLCLOCK
385              CALL GET_TUNIT(FANCY_BATCH, BC_UNIT)
386              WRITE(ERR_MSG, 1100) WALL_STOP, WT_UNIT, FANCY_BUFF, BF_UNIT, &
387                 FANCY_BATCH, BC_UNIT
388              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
389           ENDIF
390     
391      1100 FORMAT(2/,15('='),' REQUESTED CPU TIME LIMIT REACHED ',('='),/   &
392              'Batch Wall Time:',3X,F9.2,1X,A,/'Elapsed Wall Time: ',F9.2,  &
393              1X,A,/'Term Buffer:',7X,F9.2,A,/15('='),' REQUESTED CPU ',    &
394              'TIME LIMIT REACHED ',('='))
395     
396     ! Report that the hault signal was detected.
397           IF(USER_HAULT) THEN
398              WRITE(ERR_MSG, 1200)
399              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
400           ENDIF
401     
402      1200 FORMAT(2/,19('='),' MFIX STOP SIGNAL DETECTED ',19('='),/'MFIX.',&
403              'STOP file detected in run directory. Terminating MFIX.',/    &
404              'Please DO NOT FORGET to erase the MFIX.STOP file before ',   &
405              'restarting',/19('='),'MFIX STOP SIGNAL DETECTED ',19('='))
406     
407     ! This routine was restructured so all MPI ranks to the same action. As
408     ! a result, broadcasting the BATCHQ flag may not be needed.
409           eofBATCHQ = (WALL_HAULT .OR. USER_HAULT)
410           call bcast (eofBATCHQ,PE_IO)
411     
412           END SUBROUTINE CHECK_BATCH_QUEUE_END
413     
414           END SUBROUTINE TIME_MARCH
415     
416