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 ! Modules 12 !---------------------------------------------------------------------// 13 use calc_collision_wall, only: calc_dem_thermo_with_wall_stl 14 use des_bc, only: DEM_BCMI, DEM_BCMO 15 use des_thermo, only: CALC_RADT_DES 16 use desgrid, only: desgrid_pic 17 use discretelement 18 use error_manager 19 use functions 20 use machine 21 use mpi_funs_des, only: DESMPI_SEND_RECV_FIELD_VARS 22 use mpi_funs_des, only: DES_PAR_EXCHANGE 23 use mpi_utility 24 use output_man, only: OUTPUT_MANAGER 25 use run, only: ANY_SPECIES_EQ 26 use run, only: CALL_USR 27 use run, only: ENERGY_EQ 28 use run, only: NSTEP 29 use run, only: TIME, TSTOP, DT 30 use sendrecv 31 IMPLICIT NONE 32 33 ! Local variables 34 !---------------------------------------------------------------------// 35 ! Total number of particles 36 INTEGER, SAVE :: NP=0 37 38 !type(sap_t) :: sap 39 40 ! time step loop counter index 41 INTEGER :: NN,ii,nnn 42 ! loop counter index for any initial particle settling incoupled cases 43 INTEGER :: FACTOR 44 ! Temporary variables when des_continuum_coupled is T to track 45 ! changes in solid time step 46 DOUBLE PRECISION :: DTSOLID_TMP 47 ! Numbers to calculate wall time spent in DEM calculations. 48 DOUBLE PRECISION :: TMP_WALL 49 50 !......................................................................! 51 52 ! In case of restarts assign S_TIME from MFIX TIME 53 S_TIME = TIME 54 DTSOLID_TMP = DTSOLID 55 TMP_WALL = WALL_TIME() 56 57 ! Initialize time stepping variables for coupled gas/solids simulations. 58 IF(DES_CONTINUUM_COUPLED) THEN 59 IF(DT.GE.DTSOLID) THEN 60 FACTOR = CEILING(real(DT/DTSOLID)) 61 ELSE 62 FACTOR = 1 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 IF(ENERGY_EQ) CALL CONV_GS_DES1 92 IF(ANY_SPECIES_EQ) CALL DES_REACTION_MODEL 93 ELSE 94 IF(ANY_SPECIES_EQ) CALL ZERO_RRATE_DES 95 IF(ENERGY_EQ) CALL ZERO_ENERGY_SOURCE 96 ENDIF 97 CALL CALC_PG_GRAD 98 ENDIF 99 100 IF(any(CALC_RADT_DES)) CALL CALC_avgTs 101 102 103 ! Main DEM time loop 104 !----------------------------------------------------------------->>> 105 DO NN = 1, FACTOR 106 107 IF(DES_CONTINUUM_COUPLED) THEN 108 ! If the current time in the discrete loop exceeds the current time in 109 ! the continuum simulation, exit the discrete loop 110 IF(S_TIME.GE.(TIME+DT)) EXIT 111 ! If next time step in the discrete loop will exceed the current time 112 ! in the continuum simulation, modify the discrete time step so final 113 ! time will match 114 IF((S_TIME+DTSOLID).GT.(TIME+DT)) & 115 DTSOLID = TIME + DT - S_TIME 116 ENDIF 117 118 ! Calculate inter particle forces acting (collisional, cohesion) 119 CALL CALC_FORCE_DEM 120 ! Calculate or distribute fluid-particle drag force. 121 CALL CALC_DRAG_DES 122 ! Calculate heat conduction to/from wall 123 IF(ENERGY_EQ)CALL CALC_DEM_THERMO_WITH_WALL_STL 124 125 ! Update the old values of particle position and velocity with the new 126 ! values computed 127 IF (DO_OLD) CALL CFUPDATEOLD 128 ! Calculate thermochemical sources (energy and rates of formation). 129 CALL CALC_THERMO_DES 130 ! Call user functions. 131 132 IF(CALL_USR) CALL USR1_DES 133 134 ! sap = multisap%saps(0) 135 ! ! CHECK SORT 136 ! do ii=2, sap%x_endpoints_len 137 ! if (sap%x_endpoints(ii)%value < sap%x_endpoints(ii-1)%value) then 138 ! print *,"****************************************************************************************" 139 ! print *,"ii:",ii," endpoints(ii):",sap%x_endpoints(ii)%box_id,sap%x_endpoints(ii)%value 140 ! print *,"****************************************************************************************" 141 ! stop __LINE__ 142 ! endif 143 ! enddo 144 145 ! do nnn=0, size(multisap%saps)-1 146 ! !print *,"nnn = ",nnn 147 ! if (.not.check_boxes(multisap%saps(nnn))) stop __LINE__ 148 ! if (.not.check_sort(multisap%saps(nnn))) stop __LINE__ 149 ! enddo 150 151 ! Update position and velocities 152 CALL CFNEWVALUES 153 154 ! do nnn=0, size(multisap%saps)-1 155 ! !print *,"nnn = ",nnn 156 ! if (.not.check_boxes(multisap%saps(nnn))) stop __LINE__ 157 ! if (.not.check_sort(multisap%saps(nnn))) stop __LINE__ 158 ! enddo 159 160 ! sap = multisap%saps(0) 161 ! ! CHECK SORT 162 ! do ii=2, sap%x_endpoints_len 163 ! if (sap%x_endpoints(ii)%value < sap%x_endpoints(ii-1)%value) then 164 ! print *,"****************************************************************************************" 165 ! print *,"ii:",ii," endpoints(ii):",sap%x_endpoints(ii)%box_id,sap%x_endpoints(ii)%value 166 ! print *,"****************************************************************************************" 167 ! stop __LINE__ 168 ! endif 169 ! enddo 170 171 172 ! Update particle temperatures 173 CALL DES_THERMO_NEWVALUES 174 175 ! Set DO_NSEARCH before calling DES_PAR_EXCHANGE. 176 DO_NSEARCH = (NN == 1 .OR. MOD(NN,NEIGHBOR_SEARCH_N) == 0) 177 178 ! Add/Remove particles to the system via flow BCs. 179 IF(DEM_BCMI > 0) CALL MASS_INFLOW_DEM 180 IF(DEM_BCMO > 0) CALL MASS_OUTFLOW_DEM(DO_NSEARCH) 181 182 ! Call exchange particles - this will exchange particle crossing 183 ! boundaries as well as updates ghost particles information 184 IF (DO_NSEARCH .OR. (numPEs>1) .OR. DES_PERIODIC_WALLS) THEN 185 CALL DESGRID_PIC(.TRUE.) 186 CALL DES_PAR_EXCHANGE 187 ENDIF 188 189 IF(DO_NSEARCH) CALL NEIGHBOUR 190 191 ! Explicitly coupled simulations do not need to rebin particles to 192 ! the fluid grid every time step. However, this implies that the 193 ! fluid cell information and interpolation weights become stale. 194 IF(DES_CONTINUUM_COUPLED .AND. & 195 .NOT.DES_EXPLICITLY_COUPLED) THEN 196 ! Bin particles to fluid grid. 197 CALL PARTICLES_IN_CELL 198 ! Calculate interpolation weights 199 CALL CALC_INTERP_WEIGHTS 200 ! Calculate mean fields (EPg). 201 CALL COMP_MEAN_FIELDS 202 ENDIF 203 204 ! Update time to reflect changes 205 S_TIME = S_TIME + DTSOLID 206 207 ! The following section targets data writes for DEM only cases: 208 IF(.NOT.DES_CONTINUUM_COUPLED) THEN 209 ! Keep track of TIME and number of steps for DEM simulations 210 TIME = S_TIME 211 NSTEP = NSTEP + 1 212 ! Call the output manager to write RES and SPx data. 213 CALL OUTPUT_MANAGER(.FALSE., .FALSE.) 214 ENDIF ! end if (.not.des_continuum_coupled) 215 216 IF(CALL_USR) CALL USR2_DES 217 218 ENDDO ! end do NN = 1, FACTOR 219 220 ! END DEM time loop 221 !-----------------------------------------------------------------<<< 222 223 IF(CALL_USR) CALL USR3_DES 224 225 ! Reset the discrete time step to original value. 226 DTSOLID = DTSOLID_TMP 227 228 IF(DES_CONTINUUM_COUPLED) CALL DESMPI_SEND_RECV_FIELD_VARS 229 230 TMP_WALL = WALL_TIME() - TMP_WALL 231 IF(TMP_WALL > 1.0d-10) THEN 232 WRITE(ERR_MSG, 9000) trim(iVal(dble(FACTOR)/TMP_WALL)) 233 ELSE 234 WRITE(ERR_MSG, 9000) '+Inf' 235 ENDIF 236 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE., LOG=.FALSE.) 237 238 9000 FORMAT(' NITs/SEC = ',A) 239 240 RETURN 241 END SUBROUTINE DES_TIME_MARCH 242