init_settling_dem.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Module name: MAKE_ARRAYS_DES                                        !
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
!                                                                      !
!  Purpose: DES - allocating DES arrays
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

#include "version.inc"

      SUBROUTINE init_settling_dem

      USE desgrid, ONLY: desgrid_pic
      USE derived_types, only: multisap, boxhandle
      USE discretelement
      USE error_manager
      USE mpi_funs_des, ONLY: des_par_exchange
      USE run
      use functions, only: is_nonexistent
      use multi_sweep_and_prune, only: aabb_t, init_multisap, multisap_add, multisap_quicksort, multisap_sweep
      use geometry

      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
      INTEGER :: FACTOR, nn

      type(aabb_t) :: aabb

      real :: mins(3), maxs(3), rad

!-----------------------------------------------
! Include statement functions
!-----------------------------------------------


! Skip this routine if there are no particles.
      IF(particles == 0) RETURN
! Skip this routine if not a new run.
      IF(run_type /= 'NEW') RETURN

! Skip if not coupled.
      IF(.NOT.des_continuum_coupled) RETURN

! Write the initial configuration before settling
      IF(print_des_data .AND. nfactor>0) CALL write_des_data

      WRITE(err_msg, 1100) trim(ival(nfactor))
      CALL flush_err_msg(header=.false., footer=.false.)
 1100 FORMAT('Beginning DEM settling period: ',a,' steps.')

! Disable the coupling flag.
      des_continuum_coupled = .false.

      mins(1) = 0
      mins(2) = 0
      mins(3) = 0
      maxs(1) = xlength
      maxs(2) = ylength
      maxs(3) = zlength

#ifdef do_sap
      rad = 100*maxval(des_radius)
      print *,"rad = ",rad
      print *,"XLENGTH = ",xlength
      print *,"YLENGTH = ",ylength
      print *,"ZLENGTH = ",zlength
         call init_multisap(multisap,floor(xlength/rad),floor(ylength/rad),floor(zlength/rad),mins,maxs)
         ! initialize SAP
         do nn=1, max_pip
            if(is_nonexistent(nn)) cycle
            aabb%minendpoint(:) = des_pos_new(nn,:)-des_radius(nn)
            aabb%maxendpoint(:) = des_pos_new(nn,:)+des_radius(nn)

            if ( any(des_radius(nn)*multisap%one_over_cell_length(1:merge(2,3,no_k)) > 0.5 ) ) then
               print *,"BAD RADIUS..grid too fine, need to have radius=",des_radius(nn),"  less than half cell length= ",0.5/multisap%one_over_cell_length(:)
               error_stop __line__
            endif

            call multisap_add(multisap,aabb,nn,boxhandle(nn))
         enddo

         call multisap_quicksort(multisap)
         call multisap_sweep(multisap)
#endif

      DO factor = 1, nfactor
! calculate forces

         CALL calc_force_dem
! update particle position/velocity

         CALL cfnewvalues
! set the flag do_nsearch before calling particle in cell (for mpi)
         do_nsearch = (mod(factor,neighbor_search_n)==0)

! Bin the particles to the DES grid.
         CALL desgrid_pic(.true.)
! exchange particle crossing boundaries and updates ghost particles
         CALL des_par_exchange
! find particles on grid
         CALL particles_in_cell
! perform neighbor search
         IF(do_nsearch) CALL neighbour
      ENDDO

! Reset the comoupling flag.
      des_continuum_coupled = .true.

      WRITE(err_msg, 1200)
      CALL flush_err_msg(header=.false., footer=.false.)
 1200 FORMAT('DEM settling period complete.')

! this write_des_data is needed to properly show the initial state of
! the simulation (granular or coupled). In the coupled case, the
! particles may have 'settled' according to above.  In the granular
! case, the initial state won't be written until after the particles
! have moved without this call.
!      IF(PRINT_DES_DATA) CALL WRITE_DES_DATA

      RETURN
      END SUBROUTINE init_settling_dem