des_time_march.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!     Subroutine: DES_TIME_MARCH                                       !
!     Author: Jay Boyalakuntla                        Date: 21-Jun-04  !
!                                                                      !
!     Purpose: Main DEM driver routine                                 !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE des_time_march

! Modules
!---------------------------------------------------------------------//
      use calc_collision_wall, only: calc_dem_thermo_with_wall_stl
      use des_bc, only: dem_bcmi, dem_bcmo
      use des_thermo, only: calc_radt_des
      use desgrid, only: desgrid_pic
      use discretelement
      use error_manager
      use functions
      use machine
      use mpi_funs_des, only: desmpi_send_recv_field_vars
      use mpi_funs_des, only: des_par_exchange
      use mpi_utility
      use output_man, only: output_manager
      use run, only: any_species_eq
      use run, only: call_usr
      use run, only: energy_eq
      use run, only: nstep
      use run, only: time, tstop, dt
      use sendrecv
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! Total number of particles
      INTEGER, SAVE :: NP=0

      !type(sap_t) :: sap

! time step loop counter index
      INTEGER :: NN,ii,nnn
! loop counter index for any initial particle settling incoupled cases
      INTEGER :: FACTOR
! Temporary variables when des_continuum_coupled is T to track
! changes in solid time step
      DOUBLE PRECISION :: DTSOLID_TMP
! Numbers to calculate wall time spent in DEM calculations.
      DOUBLE PRECISION :: TMP_WALL

!......................................................................!

! In case of restarts assign S_TIME from MFIX TIME
      s_time = time
      dtsolid_tmp = dtsolid
      tmp_wall = wall_time()

! Initialize time stepping variables for coupled gas/solids simulations.
      IF(des_continuum_coupled) THEN
         IF(dt.GE.dtsolid) THEN
            factor = ceiling(real(dt/dtsolid))
         ELSE
            factor = 1
            dtsolid = dt
         ENDIF

! Initialize time stepping variable for pure granular simulations.
      ELSE
         factor = ceiling(real((tstop-time)/dtsolid))
         dt = dtsolid
         CALL output_manager(.false., .false.)
      ENDIF   ! end if/else (des_continuum_coupled)

      np = pip - ighost_cnt
      CALL global_all_sum(np)

      IF(des_continuum_coupled) THEN
         WRITE(err_msg, 1000) trim(ival(factor)), trim(ival(np))
         CALL flush_err_msg(header=.false., footer=.false., log=.false.)
      ELSE
         WRITE(err_msg, 1100) time, dtsolid, trim(ival(factor))
         CALL flush_err_msg(header=.false., footer=.false., log=.false.)
      ENDIF
 1000 FORMAT(/'DEM NITs: ',a,3x,'Total PIP: ', a)
 1100 FORMAT(/'Time: ',g12.5,3x,'DT: ',g12.5,3x,'DEM NITs: ',a)

      IF(call_usr) CALL usr0_des

      IF(des_continuum_coupled) THEN
         IF(des_explicitly_coupled) THEN
            CALL drag_gs_des1
            IF(energy_eq) CALL conv_gs_des1
            IF(any_species_eq) CALL des_reaction_model
         ELSE
            IF(any_species_eq) CALL zero_rrate_des
            IF(energy_eq) CALL zero_energy_source
         ENDIF
         CALL calc_pg_grad
      ENDIF

      IF(any(calc_radt_des)) CALL calc_avgts


! Main DEM time loop
!----------------------------------------------------------------->>>
      DO nn = 1, factor

         IF(des_continuum_coupled) THEN
! If the current time in the discrete loop exceeds the current time in
! the continuum simulation, exit the discrete loop
            IF(s_time.GE.(time+dt)) EXIT
! If next time step in the discrete loop will exceed the current time
! in the continuum simulation, modify the discrete time step so final
! time will match
            IF((s_time+dtsolid).GT.(time+dt)) &
               dtsolid = time + dt - s_time
         ENDIF

! Calculate inter particle forces acting (collisional, cohesion)
         CALL calc_force_dem
! Calculate or distribute fluid-particle drag force.
         CALL calc_drag_des
! Calculate heat conduction to/from wall
         IF(energy_eq)CALL calc_dem_thermo_with_wall_stl

! Update the old values of particle position and velocity with the new
! values computed
         IF (do_old) CALL cfupdateold
! Calculate thermochemical sources (energy and  rates of formation).
         CALL calc_thermo_des
! Call user functions.

         IF(call_usr) CALL usr1_des

         ! sap = multisap%saps(0)
         ! ! CHECK SORT
         ! do ii=2, sap%x_endpoints_len
         !    if (sap%x_endpoints(ii)%value < sap%x_endpoints(ii-1)%value) then
         !       print *,"****************************************************************************************"
         !       print *,"ii:",ii,"  endpoints(ii):",sap%x_endpoints(ii)%box_id,sap%x_endpoints(ii)%value
         !       print *,"****************************************************************************************"
         !       stop __LINE__
         !    endif
         ! enddo

         ! do nnn=0, size(multisap%saps)-1
         !    !print *,"nnn = ",nnn
         !    if (.not.check_boxes(multisap%saps(nnn))) stop __LINE__
         !    if (.not.check_sort(multisap%saps(nnn))) stop __LINE__
         ! enddo

! Update position and velocities
         CALL cfnewvalues

         ! do nnn=0, size(multisap%saps)-1
         !    !print *,"nnn = ",nnn
         !    if (.not.check_boxes(multisap%saps(nnn))) stop __LINE__
         !    if (.not.check_sort(multisap%saps(nnn))) stop __LINE__
         ! enddo

         ! sap = multisap%saps(0)
         ! ! CHECK SORT
         ! do ii=2, sap%x_endpoints_len
         !    if (sap%x_endpoints(ii)%value < sap%x_endpoints(ii-1)%value) then
         !       print *,"****************************************************************************************"
         !       print *,"ii:",ii,"  endpoints(ii):",sap%x_endpoints(ii)%box_id,sap%x_endpoints(ii)%value
         !       print *,"****************************************************************************************"
         !       stop __LINE__
         !    endif
         ! enddo


! Update particle temperatures
         CALL des_thermo_newvalues

! Set DO_NSEARCH before calling DES_PAR_EXCHANGE.
         do_nsearch = (nn == 1 .OR. mod(nn,neighbor_search_n) == 0)

! Add/Remove particles to the system via flow BCs.
         IF(dem_bcmi > 0) CALL mass_inflow_dem
         IF(dem_bcmo > 0) CALL mass_outflow_dem(do_nsearch)

! Call exchange particles - this will exchange particle crossing
! boundaries as well as updates ghost particles information
         IF (do_nsearch .OR. (numpes>1) .OR. des_periodic_walls) THEN
            CALL desgrid_pic(.true.)
            CALL des_par_exchange
         ENDIF

         IF(do_nsearch) CALL neighbour

! Explicitly coupled simulations do not need to rebin particles to
! the fluid grid every time step. However, this implies that the
! fluid cell information and interpolation weights become stale.
         IF(des_continuum_coupled .AND. &
            .NOT.des_explicitly_coupled) THEN
! Bin particles to fluid grid.
            CALL particles_in_cell
! Calculate interpolation weights
            CALL calc_interp_weights
! Calculate mean fields (EPg).
            CALL comp_mean_fields
         ENDIF

! Update time to reflect changes
         s_time = s_time + dtsolid

! The following section targets data writes for DEM only cases:
         IF(.NOT.des_continuum_coupled) THEN
! Keep track of TIME and number of steps for DEM simulations
            time = s_time
            nstep = nstep + 1
! Call the output manager to write RES and SPx data.
            CALL output_manager(.false., .false.)
         ENDIF  ! end if (.not.des_continuum_coupled)

         IF(call_usr) CALL usr2_des

      ENDDO ! end do NN = 1, FACTOR

! END DEM time loop
!-----------------------------------------------------------------<<<

      IF(call_usr) CALL usr3_des

! Reset the discrete time step to original value.
      dtsolid = dtsolid_tmp

      IF(des_continuum_coupled) CALL desmpi_send_recv_field_vars

      tmp_wall = wall_time() - tmp_wall
      IF(tmp_wall > 1.0d-10) THEN
         WRITE(err_msg, 9000) trim(ival(dble(factor)/tmp_wall))
      ELSE
         WRITE(err_msg, 9000) '+Inf'
      ENDIF
      CALL flush_err_msg(header=.false., footer=.false., log=.false.)

 9000 FORMAT('    NITs/SEC = ',a)

      RETURN
      END SUBROUTINE des_time_march