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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !     Subroutine: DES_TIME_MARCH                                       !
4     !     Author: Jay Boyalakuntla                        Date: 21-Jun-04  !
5     !                                                                      !
6     !     Purpose: Main DEM driver routine                                 !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE DES_TIME_MARCH
10     
11           use des_bc, only: DEM_BCMI, DEM_BCMO
12           use fldvar, only: EP_g, ROP_g, ROP_s
13           use mfix_pic, only: MPPIC
14           use output, only: SPX_DT
15           use run, only: ANY_SPECIES_EQ
16           use run, only: CALL_USR
17           use run, only: ENERGY_EQ
18           use run, only: RUN_TYPE
19           use run, only: ANY_SPECIES_EQ
20     
21           use output, only: SPX_DT
22           USE fldvar, only: EP_g, ROP_g, ROP_s
23           use des_thermo, only: DES_ENERGY_SOURCE
24           use run, only: TIME, TSTOP, DT
25     
26           use mpi_funs_des, only: DES_PAR_EXCHANGE
27     
28           use discretelement
29           use error_manager
30           use functions
31           use mpi_utility
32           use sendrecv
33           use stl
34           use vtp
35           use vtk
36     
37           IMPLICIT NONE
38     !------------------------------------------------
39     ! Local variables
40     !------------------------------------------------
41     ! Total number of particles
42           INTEGER, SAVE :: NP=0
43     
44     ! time step loop counter index
45           INTEGER :: NN
46     
47     ! loop counter index for any initial particle settling incoupled cases
48           INTEGER :: FACTOR
49     
50           INTEGER :: II, JJ, KK, IJK, CELL_ID, I_CELL, J_CELL, K_CELL, COUNT, NF,L
51           INTEGER :: IMINUS1, IPLUS1, JMINUS1, JPLUS1, KMINUS1, KPLUS1, PHASELL, LOC_MIN_PIP
52     
53     ! Local variables to keep track of time when dem restart and des
54     ! write data need to be written when des_continuum_coupled is F
55           DOUBLE PRECISION :: DES_RES_TIME, DES_SPX_TIME, NTIME
56     
57     ! Temporary variables when des_continuum_coupled is T to track
58     ! changes in solid time step
59           DOUBLE PRECISION :: TMP_DTS, DTSOLID_TMP
60     
61     ! Numbers to calculate wall time spent in DEM calculations.
62           DOUBLE PRECISION :: WALL_TIME, TMP_WALL
63     
64     
65     ! In case of restarts assign S_TIME from MFIX TIME
66           S_TIME = TIME
67           TMP_DTS = ZERO
68           DTSOLID_TMP = ZERO
69           TMP_WALL = WALL_TIME()
70     
71     ! Initialize time stepping variables for coupled gas/solids simulations.
72           IF(DES_CONTINUUM_COUPLED) THEN
73              IF(DT.GE.DTSOLID) THEN
74                 FACTOR = CEILING(real(DT/DTSOLID))
75              ELSE
76                 FACTOR = 1
77                 DTSOLID_TMP = DTSOLID
78                 DTSOLID = DT
79              ENDIF
80     
81     ! Initialize time stepping variable for pure granular simulations.
82           ELSE
83              FACTOR = CEILING(real((TSTOP-TIME)/DTSOLID))
84     
85     ! Set the DES_SPX and RES variables.
86              IF(RUN_TYPE .EQ. 'NEW') THEN
87                 DES_SPX_TIME = S_TIME
88                 DES_RES_TIME = S_TIME
89              ELSE
90                 NTIME = (S_TIME+0.1d0*DTSOLID)
91                 DES_SPX_TIME = (INT(NTIME/DES_SPX_DT) + 1)*DES_SPX_DT
92                 DES_RES_TIME = (INT(NTIME/DES_RES_DT) + 1)*DES_RES_DT
93              ENDIF
94           ENDIF   ! end if/else (des_continuum_coupled)
95     
96           NP = PIP - IGHOST_CNT
97           CALL GLOBAL_ALL_SUM(NP)
98     
99           WRITE(ERR_MSG, 1000) trim(iVal(factor)), trim(iVAL(NP))
100           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
101     
102      1000 FORMAT(/'DEM NITs: ',A,3x,'Total PIP: ', A)
103     
104           IF(CALL_USR) CALL USR0_DES
105     
106           IF(DES_CONTINUUM_COUPLED) THEN
107              DES_SPX_DT = SPX_DT(1)
108              CALL CALC_PG_GRAD
109              IF(.NOT.DES_EXPLICITLY_COUPLED) THEN
110                 IF(ANY_SPECIES_EQ) CALL ZERO_RRATE_DES
111                 IF(ENERGY_EQ) CALL ZERO_ENERGY_SOURCE
112              ENDIF
113           ENDIF
114     
115     ! Main DEM time loop
116     !----------------------------------------------------------------->>>
117           DO NN = 1, FACTOR
118     
119              IF(DES_CONTINUUM_COUPLED) THEN
120     ! If the current time in the discrete loop exceeds the current time in
121     ! the continuum simulation, exit the discrete loop
122                 IF(S_TIME.GE.(TIME+DT)) EXIT
123     ! If next time step in the discrete loop will exceed the current time
124     ! in the continuum simulation, modify the discrete time step so final
125     ! time will match
126                 IF((S_TIME+DTSOLID).GT.(TIME+DT)) THEN
127                    TMP_DTS = DTSOLID
128                    DTSOLID = TIME + DT - S_TIME
129                 ENDIF
130              ENDIF
131     
132     ! Calculate forces acting on particles (collisions, drag, etc).
133              CALL CALC_FORCE_DEM
134     ! Calculate or distribute fluid-particle drag force.
135              CALL CALC_DRAG_DES
136     
137     ! Update the old values of particle position and velocity with the new
138     ! values computed
139              IF (DO_OLD) CALL CFUPDATEOLD
140     ! Calculate thermochemical sources (energy and  rates of formation).
141              CALL CALC_THERMO_DES
142     ! Call user functions.
143              IF(CALL_USR) CALL USR1_DES
144     ! Update position and velocities
145              CALL CFNEWVALUES
146     ! Update particle temperatures
147              CALL DES_THERMO_NEWVALUES
148     ! Update particle from reactive chemistry process.
149              CALL DES_REACTION_MODEL
150     
151     ! Set DO_NSEARCH before calling DES_PAR_EXCHANGE.
152              DO_NSEARCH = (NN == 1 .OR. MOD(NN,NEIGHBOR_SEARCH_N) == 0)
153     
154     ! Add/Remove particles to the system via flow BCs.
155              IF(DEM_BCMI > 0) CALL MASS_INFLOW_DEM
156              IF(DEM_BCMO > 0) CALL MASS_OUTFLOW_DEM(DO_NSEARCH)
157     
158     ! Call exchange particles - this will exchange particle crossing
159     ! boundaries as well as updates ghost particles information
160              IF(DO_NSEARCH) THEN
161                 CALL DES_PAR_EXCHANGE
162                 CALL NEIGHBOUR
163              ELSEIF ((numPEs>1) .OR. DES_PERIODIC_WALLS) THEN
164                 CALL DES_PAR_EXCHANGE
165              ENDIF
166     
167     ! Explicitly coupled simulations do not need to rebin particles to
168     ! the fluid grid every time step. However, this implies that the
169     ! fluid cell information and interpolation weights become stale.
170              IF(.NOT.DES_EXPLICITLY_COUPLED) THEN
171     ! Bin particles to fluid grid.
172                 CALL PARTICLES_IN_CELL
173     ! Calculate interpolation weights
174                 CALL CALC_INTERP_WEIGHTS
175     ! Calculate mean fields (EPg).
176                 CALL COMP_MEAN_FIELDS
177              ENDIF
178     
179     ! Update time to reflect changes
180              S_TIME = S_TIME + DTSOLID
181     
182     ! When coupled, all write calls are made in time_march (the continuum
183     ! portion) according to user settings for spx_time and res_time.
184     ! The following section targets data writes for DEM only cases:
185              IF(.NOT.DES_CONTINUUM_COUPLED) THEN
186     ! Keep track of TIME for DEM simulations
187                 TIME = S_TIME
188     ! Calculate the 'next time' output is written.
189                 NTIME = S_TIME + 0.1d0*DTSOLID
190     
191     ! Write data using des_spx_time and des_res_time; note the time will
192     ! reflect current position of particles
193                 IF(PRINT_DES_DATA) THEN
194                    IF((NTIME >= DES_SPX_TIME) .OR. (NTIME >= TSTOP)&
195                       .OR. (NN == FACTOR) ) THEN
196                       DES_SPX_TIME =(INT(NTIME/DES_SPX_DT)+1)*DES_SPX_DT
197                       CALL WRITE_DES_DATA
198                    ENDIF
199                 ENDIF
200     
201     ! write vtp files
202           IF(WRITE_VTK_FILES) THEN
203              DO L = 1, DIMENSION_VTK
204                 IF (VTK_TIME(L)/=UNDEFINED .AND.  TIME+0.1d0*DT>=VTK_TIME(L)) THEN
205                    VTK_TIME(L) = (INT((TIME + 0.1d0*DT)/VTK_DT(L))+1)*VTK_DT(L)
206                    CALL WRITE_VTP_FILE(L)
207                 ENDIF
208              ENDDO
209           ENDIF
210     
211     ! Write out the restart infomration here. Note that there is a call
212     ! to write out the continuum variables.
213                 IF((NTIME >= DES_RES_TIME) .OR. (NTIME >= TSTOP) .OR.      &
214                    (NN == FACTOR)) THEN
215                    DES_RES_TIME = (INT(NTIME/DES_RES_DT)+1)*DES_RES_DT
216                    CALL WRITE_RES0_DES
217                    CALL WRITE_RES1
218                 ENDIF
219              ENDIF  ! end if (.not.des_continuum_coupled)
220     
221     
222              IF(CALL_USR) CALL USR2_DES
223     
224           ENDDO     ! end do NN = 1, FACTOR
225     ! END DEM time loop
226     !-----------------------------------------------------------------<<<
227     
228           IF(CALL_USR) CALL USR3_DES
229     
230     !      CALL DIFFUSE_MEAN_FIELDS
231     !      CALL CALC_EPG_DES
232     
233     ! When coupled, and if needed, reset the discrete time step accordingly
234           IF(DT.LT.DTSOLID_TMP) THEN
235              DTSOLID = DTSOLID_TMP
236           ENDIF
237     
238           IF(TMP_DTS.NE.ZERO) THEN
239              DTSOLID = TMP_DTS
240              TMP_DTS = ZERO
241           ENDIF
242     
243           IF(.NOT.DES_CONTINUUM_COUPLED)THEN
244              WRITE(ERR_MSG,"('<---------- END DES_TIME_MARCH ----------')")
245              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
246           ELSE
247              call send_recv(ep_g,2)
248              call send_recv(rop_g,2)
249              call send_recv(des_u_s,2)
250              call send_recv(des_v_s,2)
251              if(do_K) call send_recv(des_w_s,2)
252              call send_recv(rop_s,2)
253              if(ENERGY_EQ) call send_recv(des_energy_source,2)
254     
255              TMP_WALL = WALL_TIME() - TMP_WALL
256              IF(TMP_WALL > 1.0d-10) THEN
257                 WRITE(ERR_MSG, 9000) trim(iVal(dble(FACTOR)/TMP_WALL))
258              ELSE
259                 WRITE(ERR_MSG, 9000) '+Inf'
260              ENDIF
261              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
262     
263      9000 FORMAT('    NITs/SEC = ',A)
264     
265           ENDIF
266     
267           RETURN
268           END SUBROUTINE DES_TIME_MARCH
269