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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CFNEWVALUES                                            C
4     !
5     !  Purpose: DES - Calculate the new values of particle velocity,
6     !           position, angular velocity etc
7     !
8     !                                                                      C
9     !  Comments: Implements Eqns 1, 2, 3, 4 & 5  from the following paper:
10     !    Tsuji Y., Kawaguchi T., and Tanak T., "Lagrangian numerical
11     !    simulation of plug glow of cohesionless particles in a
12     !    horizontal pipe", Powder technology, 71, 239-250, 1992
13     !
14     !  pradeep : changes for parallel processing
15     !          1. periodic boundaries might lie in different proc. so adjust
16     !             particle position for periodic removed
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19           SUBROUTINE CFNEWVALUES
20     
21     !-----------------------------------------------
22     ! Modules
23     !-----------------------------------------------
24           USE constant
25           USE des_bc
26           USE derived_types, only: multisap, boxhandle
27           USE discretelement
28           USE fldvar
29           USE mfix_pic
30           USE mpi_utility
31           USE parallel
32           USE run
33           USE param
34           USE param1
35           USE physprop
36           use geometry, only: DO_K, NO_K
37           use multi_sweep_and_prune, only: aabb_t, multisap_sort, multisap_update
38     
39           IMPLICIT NONE
40     !-----------------------------------------------
41     ! Local Variables
42     !-----------------------------------------------
43           INTEGER :: L
44           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
45           DOUBLE PRECISION :: OMEGA_MAG,OMEGA_UNIT(3),ROT_ANGLE
46     
47           type(aabb_t) aabb
48     
49     !-----------------------------------------------
50     
51     ! Adams-Bashforth defaults to Euler for the first time step.
52           IF(FIRST_PASS .AND. INTG_ADAMS_BASHFORTH) THEN
53              DO L =1, MAX_PIP
54                 IF(IS_NONEXISTENT(L)) CYCLE                       ! Only real particles
55                 IF(IS_ENTERING(L).or.IS_ENTERING_GHOST(L)) CYCLE  ! Only non-entering
56                 IF(IS_GHOST(L)) CYCLE                             ! Skip ghost particles
57                 DES_ACC_OLD(L,:) = FC(L,:)/PMASS(L) + GRAV(:)
58                 ROT_ACC_OLD(L,:) = TOW(L,:)
59              ENDDO
60           ENDIF
61     
62     
63     !!$omp parallel default(none)                    &
64     !!$omp shared(MAX_PIP,INTG_EULER,INTG_ADAMS_BASHFORTH,fc,tow,do_nsearch,   &
65     !!$omp       omega_new,omega_old,pmass,grav,des_vel_new,des_pos_new,       &
66     !!$omp       des_vel_old,des_pos_old,dtsolid,omoi,des_acc_old,rot_acc_old, &
67     !!$omp       ppos,neighbor_search_rad_ratio,des_radius,DO_OLD, iGlobal_ID, &
68     !!$omp       particle_orientation,orientation,boxhandle,multisap,particle_state) &
69     !!$omp private(l,rot_angle,omega_mag,omega_unit,aabb)
70     
71     ! If a particle is classified as new, then forces are ignored.
72     ! Classification from new to existing is performed in routine
73     ! des_check_new_particle.f
74     
75     ! Advance particle position, velocity
76     ! first-order method
77           IF (INTG_EULER) THEN
78     !!$omp sections
79     !!$omp section
80              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) DES_VEL_NEW(:,1) =   &
81                 DES_VEL_NEW(:,1) + DTSOLID*(FC(:,1)/PMASS(:) + GRAV(1))
82              WHERE(PARTICLE_STATE < NORMAL_GHOST) DES_POS_NEW(:,1) =       &
83                 DES_POS_NEW(:,1) + DES_VEL_NEW(:,1)*DTSOLID
84              FC(:,1) = ZERO
85     
86     !!$omp section
87              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) DES_VEL_NEW(:,2) =   &
88                 DES_VEL_NEW(:,2) + DTSOLID*(FC(:,2)/PMASS(:) + GRAV(2))
89              WHERE(PARTICLE_STATE < NORMAL_GHOST) DES_POS_NEW(:,2) =       &
90                 DES_POS_NEW(:,2) + DES_VEL_NEW(:,2)*DTSOLID
91              FC(:,2) = ZERO
92     
93     !!$omp section
94              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) DES_VEL_NEW(:,3) =   &
95                 DES_VEL_NEW(:,3) + DTSOLID*(FC(:,3)/PMASS(:) + GRAV(3))
96              WHERE(PARTICLE_STATE < NORMAL_GHOST) DES_POS_NEW(:,3) =       &
97                 DES_POS_NEW(:,3) + DES_VEL_NEW(:,3)*DTSOLID
98              FC(:,3) = ZERO
99     
100     !!$omp section
101              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) OMEGA_NEW(:,1) =     &
102                 OMEGA_NEW(:,1) + TOW(:,1)*OMOI(:)*DTSOLID
103              TOW(:,1) = ZERO
104     
105     !!$omp section
106              WHERE(PARTICLE_STATE == NORMAL_PARTICLE)OMEGA_NEW(:,2) =      &
107                 OMEGA_NEW(:,2) + TOW(:,2)*OMOI(:)*DTSOLID
108              TOW(:,2) = ZERO
109     
110     !!$omp section
111              WHERE(PARTICLE_STATE == NORMAL_PARTICLE)OMEGA_NEW(:,3) =      &
112                 OMEGA_NEW(:,3) + TOW(:,3)*OMOI(:)*DTSOLID
113              TOW(:,3) = ZERO
114     !!$omp end sections
115     
116     ! Second-order Adams-Bashforth/Trapezoidal scheme
117           ELSEIF (INTG_ADAMS_BASHFORTH) THEN
118     
119     !!$omp sections
120     !!$omp section
121              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
122                 FC(:MAX_PIP,1) = FC(:MAX_PIP,1)/PMASS(:MAX_PIP) + GRAV(1)
123              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
124                 DES_VEL_NEW(:MAX_PIP,1) = DES_VEL_OLD(:MAX_PIP,1) + 0.5d0* &
125                    (3.d0*FC(:MAX_PIP,1)-DES_ACC_OLD(:MAX_PIP,1) )*DTSOLID
126                 DES_ACC_OLD(:MAX_PIP,1) = FC(:MAX_PIP,1)
127     
128              WHERE(PARTICLE_STATE < NORMAL_GHOST)                          &
129                 DES_POS_NEW(:MAX_PIP,1) = DES_POS_OLD(:MAX_PIP,1) + 0.5d0* &
130                    (DES_VEL_OLD(:MAX_PIP,1)+DES_VEL_NEW(:MAX_PIP,1))*DTSOLID
131              FC(:MAX_PIP,1) = ZERO
132     
133     !!$omp section
134              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
135                 FC(:MAX_PIP,2) = FC(:MAX_PIP,2)/PMASS(:MAX_PIP) + GRAV(2)
136              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
137                 DES_VEL_NEW(:MAX_PIP,2) = DES_VEL_OLD(:MAX_PIP,2) + 0.5d0* &
138                    (3.d0*FC(:MAX_PIP,2)-DES_ACC_OLD(:MAX_PIP,2) )*DTSOLID
139                 DES_ACC_OLD(:MAX_PIP,2) = FC(:MAX_PIP,2)
140     
141              WHERE(PARTICLE_STATE < NORMAL_GHOST)                          &
142                 DES_POS_NEW(:MAX_PIP,2) = DES_POS_OLD(:MAX_PIP,2) + 0.5d0* &
143                    (DES_VEL_OLD(:MAX_PIP,2)+DES_VEL_NEW(:MAX_PIP,2))*DTSOLID
144              FC(:MAX_PIP,2) = ZERO
145     
146     !!$omp section
147              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
148                 FC(:MAX_PIP,3) = FC(:MAX_PIP,3)/PMASS(:MAX_PIP) + GRAV(3)
149              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
150                 DES_VEL_NEW(:MAX_PIP,3) = DES_VEL_OLD(:MAX_PIP,3) + 0.5d0* &
151                      (3.d0*FC(:MAX_PIP,3)-DES_ACC_OLD(:MAX_PIP,3) )*DTSOLID
152                 DES_ACC_OLD(:MAX_PIP,3) = FC(:MAX_PIP,3)
153     
154              WHERE(PARTICLE_STATE < NORMAL_GHOST)                          &
155                 DES_POS_NEW(:MAX_PIP,3) = DES_POS_OLD(:MAX_PIP,3) + 0.5d0* &
156                    (DES_VEL_OLD(:MAX_PIP,3)+DES_VEL_NEW(:MAX_PIP,3))*DTSOLID
157              FC(:MAX_PIP,3) = ZERO
158     
159     !!$omp section
160             WHERE(PARTICLE_STATE(:MAX_PIP) == NORMAL_PARTICLE) &
161                OMEGA_NEW(:MAX_PIP,1) = OMEGA_OLD(:MAX_PIP,1) + 0.5d0*     &
162                    (3.d0*TOW(:MAX_PIP,1)*OMOI(:MAX_PIP) -                  &
163                    ROT_ACC_OLD(:MAX_PIP,1))*DTSOLID
164             WHERE(PARTICLE_STATE(:MAX_PIP) == NORMAL_PARTICLE) &
165                 ROT_ACC_OLD(:MAX_PIP,1) = TOW(:MAX_PIP,1)*OMOI(:MAX_PIP)
166              TOW(:MAX_PIP,1) = ZERO
167     
168     !!$omp section
169              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
170                 OMEGA_NEW(:MAX_PIP,2) = OMEGA_OLD(:MAX_PIP,2) + 0.5d0*     &
171                      (3.d0*TOW(:MAX_PIP,2)*OMOI(:MAX_PIP)-                 &
172                      ROT_ACC_OLD(:MAX_PIP,2) )*DTSOLID
173              WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
174                 ROT_ACC_OLD(:MAX_PIP,2) = TOW(:MAX_PIP,2)*OMOI(:MAX_PIP)
175              TOW(:MAX_PIP,2) = ZERO
176     
177     !!$omp section
178             WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
179                 OMEGA_NEW(:MAX_PIP,3) = OMEGA_OLD(:MAX_PIP,3) + 0.5d0*     &
180                    (3.d0*TOW(:MAX_PIP,3)*OMOI(:MAX_PIP)-                   &
181                    ROT_ACC_OLD(:MAX_PIP,3) )*DTSOLID
182             WHERE(PARTICLE_STATE == NORMAL_PARTICLE) &
183                 ROT_ACC_OLD(:MAX_PIP,3) = TOW(:MAX_PIP,3)*OMOI(:MAX_PIP)
184              TOW(:MAX_PIP,3) = ZERO
185     
186     !!$omp end sections
187           ENDIF
188     
189     #ifdef do_sap
190     !!$omp single
191              DO L = 1, MAX_PIP
192                 aabb%minendpoint(:) = DES_POS_NEW(L,:)-DES_RADIUS(L)
193                 aabb%maxendpoint(:) = DES_POS_NEW(L,:)+DES_RADIUS(L)
194                 call multisap_update(multisap,aabb,boxhandle(L))
195              ENDDO
196     !!$omp end single
197     #endif
198     
199     ! Update particle orientation - Always first order
200     ! When omega is non-zero, compute the rotation angle, and apply the
201     ! Rodrigues' rotation formula
202     
203           IF(PARTICLE_ORIENTATION) THEN
204              DO L = 1, MAX_PIP
205                 OMEGA_MAG = OMEGA_NEW(L,1)**2 +OMEGA_NEW(L,2)**2 + OMEGA_NEW(L,3)**2
206     
207                 IF(OMEGA_MAG>ZERO) THEN
208                    OMEGA_MAG=DSQRT(OMEGA_MAG)
209                    OMEGA_UNIT(:) = OMEGA_NEW(L,:)/OMEGA_MAG
210                    ROT_ANGLE = OMEGA_MAG * DTSOLID
211     
212                    ORIENTATION(:,L) = ORIENTATION(:,L)*DCOS(ROT_ANGLE) + &
213                       DES_CROSSPRDCT(OMEGA_UNIT,ORIENTATION(:,L))*DSIN(ROT_ANGLE) + &
214                       OMEGA_UNIT(:)*DOT_PRODUCT(OMEGA_UNIT,ORIENTATION(:,L))*&
215                       (ONE-DCOS(ROT_ANGLE))
216                 ENDIF
217              ENDDO
218           ENDIF
219     
220     ! Check if the particle has moved a distance greater than or equal to
221     ! its radius since the last time a neighbor search was called. if so,
222     ! make sure that neighbor is called in des_time_march
223           IF(.NOT.DO_NSEARCH) THEN
224     !!$omp do reduction (.or.:do_nsearch)
225              DO L = 1, MAX_PIP
226                 DO_NSEARCH = DO_NSEARCH .or. &
227                    (DES_POS_NEW(L,1) - PPOS(L,1))**2+              &
228                    (DES_POS_NEW(L,2) - PPOS(L,2))**2+              &
229                    (DES_POS_NEW(L,3) - PPOS(L,3))**2  >=           &
230                    (NEIGHBOR_SEARCH_RAD_RATIO*DES_RADIUS(L))**2
231              ENDDO
232           ENDIF
233     
234     !!$omp end parallel
235     
236     #ifdef do_sap
237                 call multisap_sort(multisap)
238     #endif
239     
240           FIRST_PASS = .FALSE.
241     
242      1002 FORMAT(/1X,70('*')//&
243              ' From: CFNEWVALUES -',/&
244              ' Message: Particle ',I10, ' moved a distance ', ES17.9, &
245              ' during a',/10X, 'single solids time step, which is ',&
246              ' greater than',/10X,'its radius: ', ES17.9)
247      1003 FORMAT(1X,70('*')/)
248     
249           RETURN
250     
251           contains
252     
253           include 'functions.inc'
254     
255           END SUBROUTINE CFNEWVALUES
256