File: RELATIVE:/../../../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 MFIX_netcdf
18           USE bc
19           USE cdist
20           USE coeff
21           USE compar
22           USE cont
23           USE cutcell
24           USE dashboard
25           USE discretelement
26           USE drag
27           USE fldvar
28           USE funits
29           USE geometry
30           USE indices
31           USE leqsol, only: SOLVER_STATISTICS, REPORT_SOLVER_STATS
32           USE output, only: RES_DT
33           USE param
34           USE param1
35           USE pgcor
36           USE physprop
37           USE pscor
38           USE qmom_kinetic_equation
39           USE run
40           USE scalars
41           USE stiff_chem, only : STIFF_CHEMISTRY, STIFF_CHEM_SOLVER
42           USE tau_g
43           USE tau_s
44           USE time_cpu
45           USE toleranc
46     ! use function MAX_VEL_INLET to compute max. velocity at inlet
47           USE utilities, ONLY: MAX_VEL_INLET
48           USE visc_g
49           USE visc_s
50           USE vshear
51           USE vtk
52           USE vtp
53           use output, only: RES_DT, NLOG
54           use interactive, only: INTERACT, INIT_INTERACTIVE_MODE
55           use adjust_dt
56     
57           IMPLICIT NONE
58     !-----------------------------------------------
59     ! Local variables
60     !-----------------------------------------------
61     ! Flag to indicate one pass through iterate for steady
62     ! state conditions.
63           LOGICAL :: FINISH
64     
65     ! Loop indices
66           INTEGER :: L, M
67     ! Error index
68           INTEGER :: IER
69     ! Number of iterations
70           INTEGER :: NIT, NIT_TOTAL
71     ! used for activating check_data_30
72           INTEGER :: NCHECK, DNCHECK
73     ! dummy logical variable for initializing adjust_dt
74           LOGICAL :: dummy
75     
76     ! AEOLUS : stop trigger mechanism to terminate MFIX normally before
77     ! batch queue terminates
78           DOUBLE PRECISION :: CPU_STOP
79     
80     ! Flag to save results and cleanly exit.
81           LOGICAL :: EXIT_SIGNAL = .FALSE.
82     
83     ! C Function
84           INTERFACE
85              SUBROUTINE CHECK_SOCKETS() BIND ( C )
86                use, INTRINSIC :: iso_c_binding
87              END SUBROUTINE CHECK_SOCKETS
88           END INTERFACE
89     
90     !-----------------------------------------------
91     ! External functions
92     !-----------------------------------------------
93     ! use function vavg_v_g to catch NaN's
94           DOUBLE PRECISION, EXTERNAL :: VAVG_U_G, VAVG_V_G, VAVG_W_G, X_vavg
95     
96     !-----------------------------------------------
97     
98           IF(AUTOMATIC_RESTART) RETURN
99     
100           FINISH  = .FALSE.
101           NCHECK  = NSTEP
102           DNCHECK = 1
103           CPU_IO  = ZERO
104           NIT_TOTAL = 0
105     
106           CALL INIT_OUTPUT_VARS
107     
108           IF(INTERACTIVE_MODE) CALL INIT_INTERACTIVE_MODE
109     
110     ! Parse residual strings
111           CALL PARSE_RESID_STRING ()
112     
113     ! Call user-defined subroutine to set constants, check data, etc.
114           IF (CALL_USR) CALL USR0
115     
116           CALL RRATES_INIT()
117     
118     ! Calculate all the coefficients once before entering the time loop
119           CALL INIT_COEFF(IER)
120     
121           DO M=1, MMAX
122              CALL ZERO_ARRAY (F_gs(1,M))
123           ENDDO
124     
125     ! Remove undefined values at wall cells for scalars
126           CALL UNDEF_2_0 (ROP_G)
127           DO M = 1, MMAX
128              CALL UNDEF_2_0 (ROP_S(1,M))
129           ENDDO
130     
131     ! Initialize d's and e's to zero
132           DO M = 0, MMAX
133              CALL ZERO_ARRAY (D_E(1,M))
134              CALL ZERO_ARRAY (D_N(1,M))
135              CALL ZERO_ARRAY (D_T(1,M))
136           ENDDO
137           CALL ZERO_ARRAY (E_E)
138           CALL ZERO_ARRAY (E_N)
139           CALL ZERO_ARRAY (E_T)
140     
141     
142     ! calculate shear velocities if periodic shear BCs are used
143           IF(SHEAR) CALL CAL_D(V_sh)
144     
145     ! Initialize check_mass_balance.  This routine is not active by default.
146     ! Specify a reporting interval (hard-wired in the routine) to activate
147     ! the routine.
148           Call check_mass_balance (0)
149     
150     ! sof modification: now it's only needed to do this once before time-loop
151     ! Mark the phase whose continuity will be solved and used to correct
152     ! void/volume fraction in calc_vol_fr (see subroutine for details)
153           CALL MARK_PHASE_4_COR (PHASE_4_P_G, PHASE_4_P_S, DO_CONT, MCP,&
154               DO_P_S, SWITCH_4_P_G, SWITCH_4_P_S)
155     
156     ! uncoupled discrete element simulations do not need to be within
157     ! the two fluid model time-loop
158           IF(DISCRETE_ELEMENT.AND.(.NOT.DES_CONTINUUM_COUPLED))  THEN
159              IF(WRITE_VTK_FILES) THEN
160                 DO L = 1, DIMENSION_VTK
161                    ! CALL WRITE_VTP_FILE(L)
162                 ENDDO
163              ENDIF
164              CALL DES_TIME_MARCH
165              CALL CPU_TIME(CPU_STOP)
166              CPU_STOP = CPU_STOP - CPU00
167              RETURN
168           ENDIF
169     
170     
171     ! The TIME loop begins here.............................................
172      100  CONTINUE
173     
174     
175     #ifdef socket
176           IF(INTERACTIVE_MODE) THEN
177              write(*,*) "INTERACTIVE_MODE: ",INTERACTIVE_MODE
178              CALL CHECK_SOCKETS()
179              CALL INTERACT(EXIT_SIGNAL, ABORT_SIGNAL)
180              IF(ABORT_SIGNAL) RETURN
181           ENDIF
182     #endif
183     
184     ! Terminate MFIX normally before batch queue terminates.
185           IF (CHK_BATCHQ_END) CALL CHECK_BATCH_QUEUE_END(EXIT_SIGNAL)
186     
187           IF (CALL_USR) CALL USR1
188     
189     ! Remove solids from cells containing very small quantities of solids
190           IF(.NOT.(DISCRETE_ELEMENT .OR. QMOMK) .OR. &
191              DES_CONTINUUM_HYBRID) THEN
192              IF(KT_TYPE_ENUM == GHD_2007) THEN
193                 CALL ADJUST_EPS_GHD
194              ELSE
195                 CALL ADJUST_EPS
196              ENDIF
197           ENDIF
198     
199     ! sof modification: uncomment code below and modify MARK_PHASE_4_COR to
200     ! use previous MFIX algorithm. Nov 22 2010.
201     ! Mark the phase whose continuity will be solved and used to correct
202     ! void/volume fraction in calc_vol_fr (see subroutine for details)
203     !      CALL MARK_PHASE_4_COR (PHASE_4_P_G, PHASE_4_P_S, DO_CONT, MCP,&
204     !           DO_P_S, SWITCH_4_P_G, SWITCH_4_P_S, IER)
205     
206     ! Set wall boundary conditions and transient flow b.c.'s
207           CALL SET_BC1
208     
209           CALL OUTPUT_MANAGER(EXIT_SIGNAL, FINISH)
210     
211           IF (DT == UNDEFINED) THEN
212              IF (FINISH) THEN
213                 RETURN
214              ELSE
215                 FINISH = .TRUE.
216              ENDIF
217     
218     ! Mechanism to terminate MFIX normally before batch queue terminates.
219           ELSEIF (TIME + 0.1d0*DT >= TSTOP .OR. EXIT_SIGNAL) THEN
220              IF(SOLVER_STATISTICS) &
221                 CALL REPORT_SOLVER_STATS(NIT_TOTAL, NSTEP)
222              RETURN
223           ENDIF
224     
225     ! Update previous-time-step values of field variables
226           CALL UPDATE_OLD
227     
228     ! Calculate coefficients
229           CALL CALC_COEFF_ALL (0, IER)
230     
231     ! Calculate the stress tensor trace and cross terms for all phases.
232           CALL CALC_TRD_AND_TAU()
233     
234     ! Calculate additional solid phase momentum source terms
235     ! that arise from kinetic theory constitutive relations
236           IF (.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID) THEN
237              CALL CALC_KTMOMSOURCE_U_S ()
238              CALL CALC_KTMOMSOURCE_V_S ()
239              CALL CALC_KTMOMSOURCE_W_S ()
240           ENDIF
241     
242     ! Check rates and sums of mass fractions every NLOG time steps
243           IF (NSTEP == NCHECK) THEN
244              IF (DNCHECK < 256) DNCHECK = DNCHECK*2
245              NCHECK = NCHECK + DNCHECK
246     ! Upate the reaction rates for checking
247              CALL CALC_RRATE(IER)
248              CALL CHECK_DATA_30
249           ENDIF
250     
251     ! Double the timestep for 2nd order accurate time implementation
252           IF ((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
253               (CN_ON.AND.RUN_TYPE /= 'NEW' .AND. NSTEP >= (NSTEPRST+1))) THEN
254              DT = 0.5d0*DT
255              ODT = ODT * 2.0d0
256           ENDIF
257     
258     ! Check for maximum velocity at inlet to avoid convergence problems
259           MAX_INLET_VEL = 100.0d0*MAX_VEL_INLET()
260     ! if no inlet velocity is specified, use an upper limit defined in
261     ! toleranc_mod.f
262           IF(MAX_INLET_VEL <= SMALL_NUMBER) THEN
263              MAX_INLET_VEL = MAX_ALLOWED_VEL
264              IF (UNITS == 'SI') MAX_INLET_VEL = 1D-2 * MAX_ALLOWED_VEL
265           ENDIF
266     ! Scale the value using a user defined scale factor
267           MAX_INLET_VEL = MAX_INLET_VEL * MAX_INLET_VEL_FAC
268     
269     ! Advance the solution in time by iteratively solving the equations
270      150  CALL ITERATE (IER, NIT)
271     
272           IF(AUTOMATIC_RESTART) RETURN
273     
274     ! Just to Check for NaN's, Uncomment the following lines and also lines
275     ! of code in  VAVG_U_G, VAVG_V_G, VAVG_W_G to use.
276     !      X_vavg = VAVG_U_G ()
277     !      X_vavg = VAVG_V_G ()
278     !      X_vavg = VAVG_W_G ()
279     !      IF(AUTOMATIC_RESTART) RETURN
280     
281           DO WHILE (ADJUSTDT(IER,NIT))
282              CALL ITERATE (IER, NIT)
283           ENDDO
284     
285     ! A single to interupt the time step was sent for interactive mode.
286           IF(INTERUPT) GOTO 100
287     
288           IF(DT < DT_MIN) THEN
289              AUTO_RESTART = .FALSE.
290              IF(AUTO_RESTART) THEN
291                 IF(TIME <= RES_DT) THEN
292     
293      1000 FORMAT('Automatic restart not possible as Total Time < RES_DT')
294     
295                    WRITE(ERR_MSG,1000)
296                    CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
297                    AUTOMATIC_RESTART = .TRUE.
298                 ENDIF
299                 RETURN
300     
301              ELSE
302     
303      1100 FORMAT('DT < DT_MIN.  Recovery not possible!')
304     
305                 WRITE(ERR_MSG,1100)
306                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
307     
308                 IF(WRITE_DASHBOARD) THEN
309                    RUN_STATUS = 'DT < DT_MIN.  Recovery not possible!'
310                    CALL UPDATE_DASHBOARD(NIT,0.0d0,'    ')
311                 ENDIF
312                 CALL MFIX_EXIT(MyPE)
313              ENDIF
314           ENDIF
315     
316     
317     ! Stiff Chemistry Solver.
318           IF(STIFF_CHEMISTRY) THEN
319              CALL STIFF_CHEM_SOLVER(DT, IER)
320              IF(IER /= 0) THEN
321                 dummy = ADJUSTDT(IER, NIT)
322                 GOTO 150
323              ENDIF
324           ENDIF
325     
326     ! Check over mass and elemental balances.  This routine is not active by default.
327     ! Edit the routine and specify a reporting interval to activate it.
328           CALL CHECK_MASS_BALANCE (1)
329     
330     ! Other solids model implementations
331           IF (DEM_SOLIDS) CALL DES_TIME_MARCH
332           IF (PIC_SOLIDS) CALL PIC_TIME_MARCH
333           IF (QMOMK) CALL QMOMK_TIME_MARCH
334           IF (CALL_DQMOM) CALL USR_DQMOM
335     
336     ! Advance the time step and continue
337           IF((CN_ON.AND.NSTEP>1.AND.RUN_TYPE == 'NEW') .OR. &
338              (CN_ON.AND.RUN_TYPE /= 'NEW' .AND. NSTEP >= (NSTEPRST+1))) THEN
339     ! Double the timestep for 2nd order accurate time implementation
340              DT = 2.d0*DT
341              ODT = ODT * 0.5d0
342     ! Perform the explicit extrapolation for CN implementation
343              CALL CN_EXTRAPOL
344           ENDIF
345     
346           IF (DT /= UNDEFINED) THEN
347              IF(USE_DT_PREV) THEN
348                 TIME = TIME + DT_PREV
349              ELSE
350                 TIME = TIME + DT
351              ENDIF
352              USE_DT_PREV = .FALSE.
353              NSTEP = NSTEP + 1
354           ENDIF
355     
356           NIT_TOTAL = NIT_TOTAL+NIT
357     
358     ! write (*,"('Compute the Courant number')")
359     ! call get_stats(IER)
360     
361           FLUSH (6)
362     ! The TIME loop ends here....................................................
363           GOTO 100
364     
365           IF(SOLVER_STATISTICS) CALL REPORT_SOLVER_STATS(NIT_TOTAL, NSTEP)
366     
367           END SUBROUTINE TIME_MARCH
368     
369