File: RELATIVE:/../../../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 des_thermo, only: DES_ENERGY_SOURCE
13           use desgrid, only: desgrid_pic
14           use discretelement
15           use error_manager
16           use fldvar, only: EP_g, ROP_g, ROP_s
17           use fldvar, only: EP_g, ROP_g, ROP_s
18           use functions
19           use machine
20           use mpi_funs_des, only: DES_PAR_EXCHANGE
21           use mpi_utility
22           use run, only: ANY_SPECIES_EQ
23           use run, only: ANY_SPECIES_EQ
24           use run, only: CALL_USR
25           use run, only: ENERGY_EQ
26           use run, only: NSTEP
27           use run, only: TIME, TSTOP, DT
28           use sendrecv
29     
30           IMPLICIT NONE
31     !------------------------------------------------
32     ! Local variables
33     !------------------------------------------------
34     ! Total number of particles
35           INTEGER, SAVE :: NP=0
36     
37     ! time step loop counter index
38           INTEGER :: NN
39     
40     ! loop counter index for any initial particle settling incoupled cases
41           INTEGER :: FACTOR
42     
43     ! Temporary variables when des_continuum_coupled is T to track
44     ! changes in solid time step
45           DOUBLE PRECISION :: TMP_DTS, DTSOLID_TMP
46     
47     ! Numbers to calculate wall time spent in DEM calculations.
48           DOUBLE PRECISION :: TMP_WALL
49     
50     ! In case of restarts assign S_TIME from MFIX TIME
51           S_TIME = TIME
52           TMP_DTS = ZERO
53           DTSOLID_TMP = ZERO
54           TMP_WALL = WALL_TIME()
55     
56     ! Initialize time stepping variables for coupled gas/solids simulations.
57           IF(DES_CONTINUUM_COUPLED) THEN
58              IF(DT.GE.DTSOLID) THEN
59                 FACTOR = CEILING(real(DT/DTSOLID))
60              ELSE
61                 FACTOR = 1
62                 DTSOLID_TMP = DTSOLID
63                 DTSOLID = DT
64              ENDIF
65     
66     ! Initialize time stepping variable for pure granular simulations.
67           ELSE
68              FACTOR = CEILING(real((TSTOP-TIME)/DTSOLID))
69              DT = DTSOLID
70              CALL OUTPUT_MANAGER(.FALSE., .FALSE.)
71           ENDIF   ! end if/else (des_continuum_coupled)
72     
73           NP = PIP - IGHOST_CNT
74           CALL GLOBAL_ALL_SUM(NP)
75     
76           IF(DES_CONTINUUM_COUPLED) THEN
77              WRITE(ERR_MSG, 1000) trim(iVal(factor)), trim(iVAL(NP))
78              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE., LOG=.FALSE.)
79           ELSE
80              WRITE(ERR_MSG, 1100) TIME, DTSOLID, trim(iVal(factor))
81              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE., LOG=.FALSE.)
82           ENDIF
83      1000 FORMAT(/'DEM NITs: ',A,3x,'Total PIP: ', A)
84      1100 FORMAT(/'Time: ',g12.5,3x,'DT: ',g12.5,3x,'DEM NITs: ',A)
85     
86           IF(CALL_USR) CALL USR0_DES
87     
88           IF(DES_CONTINUUM_COUPLED) THEN
89              IF(DES_EXPLICITLY_COUPLED) THEN
90                 CALL DRAG_GS_DES1
91              ELSE
92                 IF(ANY_SPECIES_EQ) CALL ZERO_RRATE_DES
93                 IF(ENERGY_EQ) CALL ZERO_ENERGY_SOURCE
94              ENDIF
95              CALL CALC_PG_GRAD
96           ENDIF
97     
98     
99     ! Main DEM time loop
100     !----------------------------------------------------------------->>>
101           DO NN = 1, FACTOR
102     
103              IF(DES_CONTINUUM_COUPLED) THEN
104     ! If the current time in the discrete loop exceeds the current time in
105     ! the continuum simulation, exit the discrete loop
106                 IF(S_TIME.GE.(TIME+DT)) EXIT
107     ! If next time step in the discrete loop will exceed the current time
108     ! in the continuum simulation, modify the discrete time step so final
109     ! time will match
110                 IF((S_TIME+DTSOLID).GT.(TIME+DT)) THEN
111                    TMP_DTS = DTSOLID
112                    DTSOLID = TIME + DT - S_TIME
113                 ENDIF
114              ENDIF
115     
116     ! Calculate forces acting on particles (collisions, drag, etc).
117              CALL CALC_FORCE_DEM
118     ! Calculate or distribute fluid-particle drag force.
119              CALL CALC_DRAG_DES
120     
121     ! Update the old values of particle position and velocity with the new
122     ! values computed
123              IF (DO_OLD) CALL CFUPDATEOLD
124     ! Calculate thermochemical sources (energy and  rates of formation).
125              CALL CALC_THERMO_DES
126     ! Call user functions.
127              IF(CALL_USR) CALL USR1_DES
128     ! Update position and velocities
129              CALL CFNEWVALUES
130     ! Update particle temperatures
131              CALL DES_THERMO_NEWVALUES
132     ! Update particle from reactive chemistry process.
133              CALL DES_REACTION_MODEL
134     
135     ! Set DO_NSEARCH before calling DES_PAR_EXCHANGE.
136              DO_NSEARCH = (NN == 1 .OR. MOD(NN,NEIGHBOR_SEARCH_N) == 0)
137     
138     ! Add/Remove particles to the system via flow BCs.
139              IF(DEM_BCMI > 0) CALL MASS_INFLOW_DEM
140              IF(DEM_BCMO > 0) CALL MASS_OUTFLOW_DEM(DO_NSEARCH)
141     
142     ! Call exchange particles - this will exchange particle crossing
143     ! boundaries as well as updates ghost particles information
144              IF (DO_NSEARCH .OR. (numPEs>1) .OR. DES_PERIODIC_WALLS) THEN
145                 CALL DESGRID_PIC(.TRUE.)
146                 CALL DES_PAR_EXCHANGE
147              ENDIF
148     
149              IF(DO_NSEARCH) CALL NEIGHBOUR
150     
151     ! Explicitly coupled simulations do not need to rebin particles to
152     ! the fluid grid every time step. However, this implies that the
153     ! fluid cell information and interpolation weights become stale.
154              IF(DES_CONTINUUM_COUPLED .AND. &
155                 .NOT.DES_EXPLICITLY_COUPLED) THEN
156     ! Bin particles to fluid grid.
157                 CALL PARTICLES_IN_CELL
158     ! Calculate interpolation weights
159                 CALL CALC_INTERP_WEIGHTS
160     ! Calculate mean fields (EPg).
161                 CALL COMP_MEAN_FIELDS
162              ENDIF
163     
164     ! Update time to reflect changes
165              S_TIME = S_TIME + DTSOLID
166     
167     ! The following section targets data writes for DEM only cases:
168              IF(.NOT.DES_CONTINUUM_COUPLED) THEN
169     ! Keep track of TIME and number of steps for DEM simulations
170                 TIME = S_TIME
171                 NSTEP = NSTEP + 1
172     ! Call the output manager to write RES and SPx data.
173                 CALL OUTPUT_MANAGER(.FALSE., .FALSE.)
174              ENDIF  ! end if (.not.des_continuum_coupled)
175     
176              IF(CALL_USR) CALL USR2_DES
177     
178           ENDDO ! end do NN = 1, FACTOR
179     
180     ! END DEM time loop
181     !-----------------------------------------------------------------<<<
182     
183           IF(CALL_USR) CALL USR3_DES
184     
185     !      CALL DIFFUSE_MEAN_FIELDS
186     !      CALL CALC_EPG_DES
187     
188     ! When coupled, and if needed, reset the discrete time step accordingly
189           IF(DT.LT.DTSOLID_TMP) THEN
190              DTSOLID = DTSOLID_TMP
191           ENDIF
192     
193           IF(TMP_DTS.NE.ZERO) THEN
194              DTSOLID = TMP_DTS
195              TMP_DTS = ZERO
196           ENDIF
197     
198           IF(.NOT.DES_CONTINUUM_COUPLED)THEN
199              WRITE(ERR_MSG,"('<---------- END DES_TIME_MARCH ----------')")
200              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
201           ELSE
202              call send_recv(ep_g,2)
203              call send_recv(rop_g,2)
204              call send_recv(des_u_s,2)
205              call send_recv(des_v_s,2)
206              if(do_K) call send_recv(des_w_s,2)
207              call send_recv(rop_s,2)
208              if(ENERGY_EQ) call send_recv(des_energy_source,2)
209     
210              TMP_WALL = WALL_TIME() - TMP_WALL
211              IF(TMP_WALL > 1.0d-10) THEN
212                 WRITE(ERR_MSG, 9000) trim(iVal(dble(FACTOR)/TMP_WALL))
213              ELSE
214                 WRITE(ERR_MSG, 9000) '+Inf'
215              ENDIF
216              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE., LOG=.FALSE.)
217     
218      9000 FORMAT('    NITs/SEC = ',A)
219     
220           ENDIF
221     
222           RETURN
223           END SUBROUTINE DES_TIME_MARCH
224