File: N:\mfix\model\des\init_settling_dem.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Module name: MAKE_ARRAYS_DES                                        !
3     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
4     !                                                                      !
5     !  Purpose: DES - allocating DES arrays
6     !                                                                      !
7     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
8     
9     #include "version.inc"
10     
11           SUBROUTINE INIT_SETTLING_DEM
12     
13           USE desgrid, ONLY: desgrid_pic
14           USE derived_types, only: multisap, boxhandle
15           USE discretelement
16           USE error_manager
17           USE mpi_funs_des, ONLY: DES_PAR_EXCHANGE
18           USE run
19           use functions, only: is_nonexistent
20           use multi_sweep_and_prune, only: aabb_t, init_multisap, multisap_add, multisap_quicksort, multisap_sweep
21           use geometry
22     
23           IMPLICIT NONE
24     !-----------------------------------------------
25     ! Local variables
26     !-----------------------------------------------
27           INTEGER :: FACTOR, nn
28     
29           type(aabb_t) :: aabb
30     
31           real :: mins(3), maxs(3), rad
32     
33     !-----------------------------------------------
34     ! Include statement functions
35     !-----------------------------------------------
36     
37     
38     ! Skip this routine if there are no particles.
39           IF(PARTICLES == 0) RETURN
40     ! Skip this routine if not a new run.
41           IF(RUN_TYPE /= 'NEW') RETURN
42     
43     ! Skip if not coupled.
44           IF(.NOT.DES_CONTINUUM_COUPLED) RETURN
45     
46     ! Write the initial configuration before settling
47           IF(PRINT_DES_DATA .AND. NFACTOR>0) CALL WRITE_DES_DATA
48     
49           WRITE(ERR_MSG, 1100) trim(iVal(NFACTOR))
50           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
51      1100 FORMAT('Beginning DEM settling period: ',A,' steps.')
52     
53     ! Disable the coupling flag.
54           DES_CONTINUUM_COUPLED = .FALSE.
55     
56           mins(1) = 0
57           mins(2) = 0
58           mins(3) = 0
59           maxs(1) = XLENGTH
60           maxs(2) = YLENGTH
61           maxs(3) = ZLENGTH
62     
63     #ifdef do_sap
64           rad = 100*maxval(des_radius)
65           print *,"rad = ",rad
66           print *,"XLENGTH = ",XLENGTH
67           print *,"YLENGTH = ",YLENGTH
68           print *,"ZLENGTH = ",ZLENGTH
69              call init_multisap(multisap,floor(XLENGTH/rad),floor(YLENGTH/rad),floor(ZLENGTH/rad),mins,maxs)
70              ! initialize SAP
71              do nn=1, MAX_PIP
72                 if(is_nonexistent(nn)) cycle
73                 aabb%minendpoint(:) = DES_POS_NEW(nn,:)-DES_RADIUS(nn)
74                 aabb%maxendpoint(:) = DES_POS_NEW(nn,:)+DES_RADIUS(nn)
75     
76                 if ( any(DES_RADIUS(nn)*multisap%one_over_cell_length(1:merge(2,3,NO_K)) > 0.5 ) ) then
77                    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(:)
78                    ERROR_STOP __LINE__
79                 endif
80     
81                 call multisap_add(multisap,aabb,nn,boxhandle(nn))
82              enddo
83     
84              call multisap_quicksort(multisap)
85              call multisap_sweep(multisap)
86     #endif
87     
88           DO FACTOR = 1, NFACTOR
89     ! calculate forces
90     
91              CALL CALC_FORCE_DEM
92     ! update particle position/velocity
93     
94              CALL CFNEWVALUES
95     ! set the flag do_nsearch before calling particle in cell (for mpi)
96              DO_NSEARCH = (MOD(FACTOR,NEIGHBOR_SEARCH_N)==0)
97     
98     ! Bin the particles to the DES grid.
99              CALL DESGRID_PIC(.TRUE.)
100     ! exchange particle crossing boundaries and updates ghost particles
101              CALL DES_PAR_EXCHANGE
102     ! find particles on grid
103              CALL PARTICLES_IN_CELL
104     ! perform neighbor search
105              IF(DO_NSEARCH) CALL NEIGHBOUR
106           ENDDO
107     
108     ! Reset the comoupling flag.
109           DES_CONTINUUM_COUPLED = .TRUE.
110     
111           WRITE(ERR_MSG, 1200)
112           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
113      1200 FORMAT('DEM settling period complete.')
114     
115     ! this write_des_data is needed to properly show the initial state of
116     ! the simulation (granular or coupled). In the coupled case, the
117     ! particles may have 'settled' according to above.  In the granular
118     ! case, the initial state won't be written until after the particles
119     ! have moved without this call.
120     !      IF(PRINT_DES_DATA) CALL WRITE_DES_DATA
121     
122           RETURN
123           END SUBROUTINE INIT_SETTLING_DEM
124