cfnewvalues.f

File: RELATIVE:/../../../mfix.git/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 discretelement
27           USE fldvar
28           USE mfix_pic
29           USE mpi_utility
30           USE parallel
31           USE param
32           USE param1
33           USE physprop
34           use geometry, only: DO_K, NO_K
35           IMPLICIT NONE
36     !-----------------------------------------------
37     ! Local Variables
38     !-----------------------------------------------
39           INTEGER :: L
40           DOUBLE PRECISION :: DD(3), NEIGHBOR_SEARCH_DIST
41           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
42           DOUBLE PRECISION :: OMEGA_MAG,OMEGA_UNIT(3),ROT_ANGLE
43     !-----------------------------------------------
44     
45           IF(MPPIC) THEN
46              IF(MPPIC_SOLID_STRESS_SNIDER) THEN
47                 CALL CFNEWVALUES_MPPIC_SNIDER
48              ELSE
49     ! call the coloring function like approach
50                 CALL CFNEWVALUES_MPPIC
51              ENDIF
52              RETURN
53           ENDIF
54     
55     ! Adams-Bashforth defaults to Euler for the first time step.
56           IF(FIRST_PASS .AND. INTG_ADAMS_BASHFORTH) THEN
57              DO L =1, MAX_PIP
58                 IF(IS_NONEXISTENT(L)) CYCLE                       ! Only real particles
59                 IF(IS_ENTERING(L).or.IS_ENTERING_GHOST(L)) CYCLE  ! Only non-entering
60                 IF(IS_GHOST(L)) CYCLE                             ! Skip ghost particles
61                 DES_ACC_OLD(:,L) = FC(:,L)/PMASS(L) + GRAV(:)
62                 ROT_ACC_OLD(:,L) = TOW(:,L)
63              ENDDO
64           ENDIF
65     
66     !$omp parallel do if(max_pip .ge. 10000) default(none)                    &
67     !$omp shared(MAX_PIP,INTG_EULER,INTG_ADAMS_BASHFORTH,fc,tow,              &
68     !$omp       omega_new,omega_old,pmass,grav,des_vel_new,des_pos_new,       &
69     !$omp       des_vel_old,des_pos_old,dtsolid,omoi,des_acc_old,rot_acc_old, &
70     !$omp       ppos,neighbor_search_rad_ratio,des_radius,DO_OLD, iGlobal_ID, &
71     !$omp       particle_orientation, orientation) &
72     !$omp private(l,dd,neighbor_search_dist,rot_angle,omega_mag,omega_unit)   &
73     !$omp reduction(.or.:do_nsearch) schedule (auto)
74     
75           DO L = 1, MAX_PIP
76     ! only process particles that exist
77              IF(IS_NONEXISTENT(L)) CYCLE
78     ! skip ghost particles
79              IF(IS_GHOST(L).or.IS_ENTERING_GHOST(L).or.IS_EXITING_GHOST(L)) CYCLE
80     
81     ! If a particle is classified as new, then forces are ignored.
82     ! Classification from new to existing is performed in routine
83     ! des_check_new_particle.f
84              IF(.NOT.IS_ENTERING(L) .AND. .NOT.IS_ENTERING_GHOST(L))THEN
85                 FC(:,L) = FC(:,L)/PMASS(L) + GRAV(:)
86              ELSE
87                 FC(:,L) = ZERO
88                 TOW(:,L) = ZERO
89              ENDIF
90     
91     ! Advance particle position, velocity
92             IF (INTG_EULER) THEN
93     ! first-order method
94                 DES_VEL_NEW(:,L) = DES_VEL_NEW(:,L) + FC(:,L)*DTSOLID
95                 DD(:) = DES_VEL_NEW(:,L)*DTSOLID
96                 DES_POS_NEW(:,L) = DES_POS_NEW(:,L) + DD(:)
97                 OMEGA_NEW(:,L)   = OMEGA_NEW(:,L) + TOW(:,L)*OMOI(L)*DTSOLID
98              ELSEIF (INTG_ADAMS_BASHFORTH) THEN
99     
100     ! Second-order Adams-Bashforth/Trapezoidal scheme
101                 DES_VEL_NEW(:,L) = DES_VEL_OLD(:,L) + 0.5d0*&
102                    ( 3.d0*FC(:,L)-DES_ACC_OLD(:,L) )*DTSOLID
103                 OMEGA_NEW(:,L)   =  OMEGA_OLD(:,L) + 0.5d0*&
104                    ( 3.d0*TOW(:,L)*OMOI(L)-ROT_ACC_OLD(:,L) )*DTSOLID
105                 DD(:) = 0.5d0*( DES_VEL_OLD(:,L)+DES_VEL_NEW(:,L) )*DTSOLID
106                 DES_POS_NEW(:,L) = DES_POS_OLD(:,L) + DD(:)
107                 DES_ACC_OLD(:,L) = FC(:,L)
108                 ROT_ACC_OLD(:,L) = TOW(:,L)*OMOI(L)
109              ENDIF
110     
111     ! Update particle orientation - Always first order
112     ! When omega is non-zero, compute the rotation angle, and apply the
113     ! Rodrigues' rotation formula
114     
115              IF(PARTICLE_ORIENTATION) THEN
116                 OMEGA_MAG = OMEGA_NEW(1,L)**2 +OMEGA_NEW(2,L)**2 + OMEGA_NEW(3,L)**2
117     
118                 IF(OMEGA_MAG>ZERO) THEN
119                    OMEGA_MAG=DSQRT(OMEGA_MAG)
120                    OMEGA_UNIT(:) = OMEGA_NEW(:,L)/OMEGA_MAG
121                    ROT_ANGLE = OMEGA_MAG * DTSOLID
122     
123                    ORIENTATION(:,L) = ORIENTATION(:,L)*DCOS(ROT_ANGLE) &
124                                      + DES_CROSSPRDCT(OMEGA_UNIT,ORIENTATION(:,L))*DSIN(ROT_ANGLE) &
125                                      + OMEGA_UNIT(:)*DOT_PRODUCT(OMEGA_UNIT,ORIENTATION(:,L))*(ONE-DCOS(ROT_ANGLE))
126                 ENDIF
127              ENDIF
128     
129     ! Check if the particle has moved a distance greater than or equal to
130     ! its radius during one solids time step. if so, call stop
131              IF(dot_product(DD,DD).GE.DES_RADIUS(L)**2) THEN
132                 WRITE(*,1002) iGlobal_ID(L), sqrt(dot_product(DD,DD)), DES_RADIUS(L)
133                 IF (DO_OLD) WRITE(*,'(5X,A,3(ES17.9))') &
134                    'old particle pos = ', DES_POS_OLD(:,L)
135                 WRITE(*,'(5X,A,3(ES17.9))') &
136                    'new particle pos = ', DES_POS_NEW(:,L)
137                 WRITE(*,'(5X,A,3(ES17.9))')&
138                    'new particle vel = ', DES_VEL_NEW(:,L)
139                 WRITE(*,1003)
140                 STOP 1
141              ENDIF
142     
143     ! Check if the particle has moved a distance greater than or equal to
144     ! its radius since the last time a neighbor search was called. if so,
145     ! make sure that neighbor is called in des_time_march
146              IF(.NOT.DO_NSEARCH) THEN
147                 DD(:) = DES_POS_NEW(:,L) - PPOS(:,L)
148                 NEIGHBOR_SEARCH_DIST = NEIGHBOR_SEARCH_RAD_RATIO*&
149                    DES_RADIUS(L)
150                 IF(dot_product(DD,DD).GE.NEIGHBOR_SEARCH_DIST**2) DO_NSEARCH = .TRUE.
151              ENDIF
152     
153     ! Reset total contact force and torque
154              FC(:,L) = ZERO
155              TOW(:,L) = ZERO
156     
157           ENDDO
158     !$omp end parallel do
159     
160           FIRST_PASS = .FALSE.
161     
162     
163      1002 FORMAT(/1X,70('*')//&
164              ' From: CFNEWVALUES -',/&
165              ' Message: Particle ',I10, ' moved a distance ', ES17.9, &
166              ' during a',/10X, 'single solids time step, which is ',&
167              ' greater than',/10X,'its radius: ', ES17.9)
168      1003 FORMAT(1X,70('*')/)
169     
170           RETURN
171     
172           contains
173     
174             include 'functions.inc'
175     
176           END SUBROUTINE CFNEWVALUES
177     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
178     !                                                                      C
179     !  Module name: CFNEWVALUES_MPPIC_SNIDER                               C
180     !
181     !                                                                      C
182     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
183           SUBROUTINE CFNEWVALUES_MPPIC_SNIDER
184     
185           USE param
186           USE param1
187           USE parallel
188           USE physprop
189           USE constant
190           USE fldvar
191           USE discretelement
192           USE des_bc
193           USE mpi_utility
194           USE fldvar
195           USE cutcell
196           USE mfix_pic
197           USE randomno
198           USE geometry, only: DO_K, NO_K
199           USE fun_avg
200           USE functions
201     
202           IMPLICIT NONE
203     !-----------------------------------------------
204     ! Local Variables
205     !-----------------------------------------------
206           INTEGER L, M, IDIM
207           INTEGER I, J, K, IJK, IJK_OLD
208     
209           DOUBLE PRECISION DD(3), DIST, &
210                            DP_BAR, COEFF_EN, MEANVEL(3), D_GRIDUNITS(3)
211     
212           DOUBLE PRECISION DELUP(3), UPRIMETAU(3), UPRIMETAU_INT(3), MEAN_FREE_PATH, PS_FORCE(3)
213     ! index to track accounted for particles
214           INTEGER PC
215     
216     ! Logical for local debug warnings
217           LOGICAL DES_LOC_DEBUG
218     
219     ! maximum distance particles can move in MPPIC
220           DOUBLE PRECISION UPRIMEMOD
221     
222     ! dt's in each direction  based on cfl_pic for the mppic case
223     
224           DOUBLE PRECISION DTPIC_TMPX, DTPIC_TMPY , DTPIC_TMPZ, THREEINTOSQRT2, RAD_EFF
225     
226           LOGICAL :: DELETE_PART, INSIDE_DOMAIN
227           INTEGER :: PIP_DEL_COUNT
228     
229           INTEGER :: LPIP_DEL_COUNT_ALL(0:numPEs-1), PIJK_OLD(5)
230     
231           double precision  sig_u, mean_u
232           double precision, allocatable, dimension(:,:) ::  rand_vel
233     !-----------------------------------------------
234     
235           PC = 1
236           DTPIC_CFL = LARGE_NUMBER
237     
238           if(NO_K) THREEINTOSQRT2 = 2.D0*SQRT(2.D0)
239           if(DO_K) THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
240           THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
241           DES_VEL_MAX(:) = ZERO
242           PIP_DEL_COUNT = 0
243     
244           !EPG_MIN2 = MINVAL(EP_G(:))
245           !epg_min_loc = MINLOC(EP_G)
246           !IJK = epg_min_loc(1)
247     
248           allocate(rand_vel(3, MAX_PIP))
249           do idim = 1, merge(2,3,NO_K)
250              mean_u = zero
251              sig_u = 1.d0
252              CALL NOR_RNO(RAND_VEL(IDIM, 1:MAX_PIP), MEAN_U, SIG_U)
253           enddo
254     
255           DO L = 1, MAX_PIP
256              DELETE_PART = .false.
257     ! pradeep skip ghost particles
258              if(pc.gt.pip) exit
259              if(is_nonexistent(l)) cycle
260              pc = pc+1
261              if(is_ghost(l) .or. is_entering_ghost(l) .or. is_exiting_ghost(l)) cycle
262     
263              DES_LOC_DEBUG = .FALSE.
264     
265              !if(mppic) then
266              !   IJK = PIJK(L, 4)
267     
268              !   COUNT_BC = DES_CELLWISE_BCDATA(IJK)%COUNT_DES_BC
269              !   IF(COUNT_BC.GE.1) CALL MPPIC_ADD_FRIC_FORCE(L)
270              !ENDIF
271     ! If a particle is classified as new, then forces are ignored.
272     ! Classification from new to existing is performed in routine
273     ! des_check_new_particle.f
274              IF(.NOT.IS_ENTERING(L) .AND. .NOT.IS_ENTERING_GHOST(L))THEN
275                 FC(:,L) = FC(:,L)/PMASS(L) + GRAV(:)
276              ELSE
277                 FC(:,L) = ZERO
278              ENDIF
279     
280              !DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
281              !By comparing MFIX equations and equations in those papers,
282              !D_p = Beta/(EP_S*RHOP)
283              !F_gp in drag_fgs.f  = Beta*PVOL/EP_S
284              !Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
285              DP_BAR = F_gp(L)/(PMASS(L))
286              IF(.NOT.DES_CONTINUUM_COUPLED) DP_BAR = ZERO
287     
288              if(.not.MPPIC_PDRAG_IMPLICIT) DP_BAR = ZERO
289     
290              M = PIJK(L,5)
291              IJK = PIJK(L,4)
292     
293              IJK_OLD = IJK
294     
295              PIJK_OLD(:) = PIJK(L,:)
296     
297              I = I_OF(IJK)
298              J = J_OF(IJK)
299              K = K_OF(IJK)
300              COEFF_EN =  MPPIC_COEFF_EN1
301              UPRIMETAU(:) = ZERO
302     
303              DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L) + &
304              & FC(:,L)*DTSOLID)/(1.d0+DP_BAR*DTSOLID)
305     
306              do idim = 1, merge(2,3,NO_K)
307                 SIG_U = 0.05D0
308                 rand_vel(idim, L)  = sig_u*DES_VEL_NEW(IDIM,L)*rand_vel(idim, L)
309              enddo
310     
311              !MEANVEL(1) = DES_U_S(IJK_OLD,M)
312              !MEANVEL(2) = DES_V_S(IJK_OLD,M)
313              !IF(DO_K) MEANVEL(3) = DES_W_S(IJK_OLD,M)
314     
315              PS_FORCE(:) = PS_GRAD(L, :)
316              DELUP(:) = -( DTSOLID*PS_FORCE(:))/((1.d0+DP_BAR*DTSOLID))
317              DELUP(:) = DELUP(:)/ROP_S(IJK_OLD,M)
318     
319              MEANVEL(:) = AVGSOLVEL_P(:,L)
320     
321              DO IDIM = 1, merge(2,3,NO_K)
322                 IF(PS_FORCE(IDIM).LE.ZERO) THEN
323                    UPRIMETAU(IDIM) = MIN(DELUP(IDIM), (1+COEFF_EN)*(MEANVEL(IDIM)-DES_VEL_NEW(IDIM,L)))
324                    UPRIMETAU_INT(IDIM) = UPRIMETAU(IDIM)
325                    UPRIMETAU(IDIM) = MAX(UPRIMETAU(IDIM), ZERO)
326                 ELSE
327                    UPRIMETAU(IDIM) = MAX(DELUP(IDIM), (1+COEFF_EN)*(MEANVEL(IDIM)-DES_VEL_NEW(IDIM,L)))
328                    UPRIMETAU_INT(IDIM) = UPRIMETAU(IDIM)
329                    UPRIMETAU(IDIM) = MIN(UPRIMETAU(IDIM), ZERO)
330                 END IF
331     
332              ENDDO
333     
334              DES_VEL_NEW(:,L) = DES_VEL_NEW(:,L) + UPRIMETAU(:)
335              DD(:) = DES_VEL_NEW(:,L)*DTSOLID
336     
337              UPRIMEMOD = SQRT(DOT_PRODUCT(DES_VEL_NEW(1:,L), DES_VEL_NEW(1:,L)))
338     
339              RAD_EFF = DES_RADIUS(L)
340                    !RAD_EFF = (DES_STAT_WT(L)**(1.d0/3.d0))*DES_RADIUS(L)
341              MEAN_FREE_PATH  = MAX(1.d0/(1.d0-EP_STAR), 1.d0/(1.D0-EP_G(IJK_OLD)))
342              MEAN_FREE_PATH = MEAN_FREE_PATH*RAD_EFF/THREEINTOSQRT2
343     
344              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) then
345              !   DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L)/UPRIMEMOD)*MEAN_FREE_PATH/DTSOLID
346              !ENDIF
347     
348               DIST = SQRT(DOT_PRODUCT(DD,DD))
349     
350              CALL PIC_FIND_NEW_CELL(L)
351     
352              IJK = PIJK(L,4)
353     
354              INSIDE_DOMAIN = .true.
355     
356              INSIDE_DOMAIN = FLUID_AT(IJK)
357     
358              IF(EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN.and.IJK.NE.IJK_OLD) then
359                 DD(:) = rand_vel(:, L)*dtsolid
360                 DES_VEL_NEW(:,L) = 0.8d0*DES_VEL_NEW(:,L)
361              ENDIF
362     
363     
364     
365              PIJK(L,:) = PIJK_OLD(:)
366              DIST = SQRT(DOT_PRODUCT(DD,DD))
367     
368               IF(DIST.GT.MEAN_FREE_PATH) THEN
369                  DD(:) =  DES_VEL_NEW(:,L)*DTSOLID*MEAN_FREE_PATH/DIST
370               ENDIF
371     
372              DES_POS_NEW(:,L) = DES_POS_NEW(:,L) + DD(:)
373     
374              D_GRIDUNITS(1) = ABS(DD(1))/DX(PIJK(L,1))
375              D_GRIDUNITS(2) = ABS(DD(2))/DY(PIJK(L,2))
376              D_GRIDUNITS(3) = 1.d0
377              IF(DO_K) D_GRIDUNITS(3) = ABS(DD(3))/DZ(PIJK(L,3))
378     
379              DIST = SQRT(DOT_PRODUCT(DD,DD))
380              DTPIC_TMPX = (CFL_PIC*DX(PIJK(L,1)))/(ABS(DES_VEL_NEW(1,L))+SMALL_NUMBER)
381              DTPIC_TMPY = (CFL_PIC*DY(PIJK(L,2)))/(ABS(DES_VEL_NEW(2,L))+SMALL_NUMBER)
382              DTPIC_TMPZ = LARGE_NUMBER
383              IF(DO_K) DTPIC_TMPZ = (CFL_PIC*DZ(PIJK(L,3)))/(ABS(DES_VEL_NEW(3,L))+SMALL_NUMBER)
384     
385     
386     ! Check if the particle has moved a distance greater than or equal to grid spacing
387     ! if so, then exit
388     
389              DO IDIM = 1, merge(2,3,NO_K)
390                 IF(D_GRIDUNITS(IDIM).GT.ONE) THEN
391                    IF(DMP_LOG.OR.myPE.eq.pe_IO) THEN
392     
393                       WRITE(UNIT_LOG, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
394                       WRITE(*, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
395                       DELETE_PART = .true.
396     
397                    ENDIF
398                    !CALL mfix_exit(myPE)
399                 END IF
400     
401              END DO
402              IF(.not.DELETE_PART) DTPIC_CFL = MIN(DTPIC_CFL, DTPIC_TMPX, DTPIC_TMPY, DTPIC_TMPZ)
403              FC(:,L) = ZERO
404     
405              IF(DELETE_PART) THEN
406                 CALL SET_NONEXISTENT(l)
407                 PIP_DEL_COUNT = PIP_DEL_COUNT + 1
408              ENDIF
409              IF (DES_LOC_DEBUG) WRITE(*,1001)
410           ENDDO
411     
412     
413           IF(MPPIC) THEN
414              CALL global_all_max(DTPIC_CFL)
415              PIP = PIP - PIP_DEL_COUNT
416     
417              LPIP_DEL_COUNT_ALL(:) = 0
418              LPIP_DEL_COUNT_ALL(mype) = PIP_DEL_COUNT
419              CALL GLOBAL_ALL_SUM(LPIP_DEL_COUNT_ALL)
420              IF((DMP_LOG).AND.SUM(LPIP_DEL_COUNT_ALL(:)).GT.0) THEN
421                 IF(PRINT_DES_SCREEN) THEN
422                    WRITE(*,'(/,2x,A,2x,i10,/,A)') &
423                         'TOTAL NUMBER OF PARTS STEPPING MORE THAN ONE GRID SPACE = ', SUM(LPIP_DEL_COUNT_ALL(:)), &
424                         'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
425                 ENDIF
426     
427                 WRITE(UNIT_LOG,'(/,2x,A,2x,i10,/,A)') &
428                      'TOTAL NUMBER OF PARTS  STEPPING MORE THAN ONE GRID SPACEC= ', SUM(LPIP_DEL_COUNT_ALL(:)), &
429                      'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
430                 !DO IPROC = 0, NUMPES-1
431                 !   WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
432                 !   'PARTICLES OUTSIDE DOMAIN (PIC)  ON PROC:', IPROC,' EQUAL TO', LPIP_DEL_COUNT_ALL(IPROC)
433                 !ENDDO
434     
435              ENDIF
436     
437              DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
438     
439              IF(DTSOLID.GT.DTPIC_MAX) THEN
440                 !IF(DMP_LOG) WRITE(UNIT_LOG, 2001) DTSOLID
441                 IF(myPE.eq.pe_IO) WRITE(*, 2004) DTSOLID
442              ELSEIF(DTSOLID.LT.DTPIC_MAX) THEN
443     
444                 !IF(DMP_LOG) WRITE(UNIT_LOG, 2002) DTSOLID
445                 IF(myPE.eq.pe_IO) WRITE(*, 2002) DTPIC_MAX
446              ELSE
447                 !WRITE(*,'(A40,2x,g17.8)') 'DT
448                 !IF(mype.eq.pe_IO) WRITE(*,2003) DTSOLID
449              END IF
450     
451     
452           ENDIF
453      1001 FORMAT(3X,'<---------- END CFNEWVALUES ----------')
454     
455     2001  FORMAT(/1X,70('*'),//,10X,  &
456                & 'MOVEMENT UNDESIRED IN CFNEWVALUES: PARTICLE', i5, /,10X, &
457                & ' MOVED MORE THAN A GRID SPACING IN ONE TIME STEP', /,10X, &
458                & 'MOVEMENT IN GRID UNITS = ', 3(2x, g17.8),/,10X,  &
459                & 'TERMINAL ERROR: NOT STOPPING, BUT DELETING THE PARTICLE', &
460                & /1X,70('*'), /10X, &
461                & 'DES_VEL_NEW = ',  3(2x, g17.8))
462     
463      2004 FORMAT(/10x, &
464           & 'DTSOLID SHUD BE REDUCED TO', g17.8)
465     
466      2002 FORMAT(/10x, &
467           & 'DTSOLID CAN BE INCREASED TO', g17.8)
468     
469           RETURN
470           END SUBROUTINE CFNEWVALUES_MPPIC_SNIDER
471     
472     
473           SUBROUTINE CFNEWVALUES_MPPIC
474     
475           USE param
476           USE param1
477           USE parallel
478           USE physprop
479           USE constant
480           USE fldvar
481           USE discretelement
482           USE des_bc
483           USE mpi_utility
484           USE mfix_pic
485           USE randomno
486           USE cutcell
487           USE fun_avg
488           USE functions
489           IMPLICIT NONE
490     !-----------------------------------------------
491     ! Local Variables
492     !-----------------------------------------------
493           INTEGER L, M, IDIM
494           INTEGER I, J, K, IJK, IJK_OLD
495     
496           DOUBLE PRECISION DD(3), DIST, &
497                            DP_BAR, MEANVEL(3), D_GRIDUNITS(3)
498     
499           DOUBLE PRECISION MEAN_FREE_PATH
500     ! index to track accounted for particles
501           INTEGER PC
502     
503     ! Logical for local debug warnings
504           LOGICAL DES_LOC_DEBUG
505     
506     ! maximum distance particles can move in MPPIC
507           DOUBLE PRECISION UPRIMEMOD
508     
509     ! dt's in each direction  based on cfl_pic for the mppic case
510     
511           DOUBLE PRECISION DTPIC_TMPX, DTPIC_TMPY , DTPIC_TMPZ, THREEINTOSQRT2, RAD_EFF, POS_Z
512           DOUBLE PRECISION :: VELP_INT(3)
513     
514           LOGICAL :: DELETE_PART, INSIDE_DOMAIN
515           INTEGER :: PIP_DEL_COUNT
516     
517           INTEGER :: LPIP_DEL_COUNT_ALL(0:numPEs-1), PIJK_OLD(5)
518     
519           double precision  sig_u, mean_u, DTPIC_MIN_X,  DTPIC_MIN_Y,  DTPIC_MIN_Z
520           double precision, allocatable, dimension(:,:) ::  rand_vel
521     
522           double precision :: norm1, norm2, norm3
523           Logical :: OUTER_STABILITY_COND, DES_FIXED_BED
524     !-----------------------------------------------
525     
526           OUTER_STABILITY_COND = .false.
527           DES_FIXED_BED = .false.
528           PC = 1
529           DTPIC_CFL = LARGE_NUMBER
530           DTPIC_MIN_X =  LARGE_NUMBER
531           DTPIC_MIN_Y =  LARGE_NUMBER
532           DTPIC_MIN_Z =  LARGE_NUMBER
533     
534           if(NO_K) THREEINTOSQRT2 = 2.D0*SQRT(2.D0)
535           if(DO_K) THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
536           THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
537           DES_VEL_MAX(:) = ZERO
538           PIP_DEL_COUNT = 0
539     
540           !EPG_MIN2 = MINVAL(EP_G(:))
541           !epg_min_loc = MINLOC(EP_G)
542           !IJK = epg_min_loc(1)
543           allocate(rand_vel(3, MAX_PIP))
544           do idim = 1, merge(2,3,NO_K)
545              mean_u = zero
546              sig_u = 1.d0
547              CALL NOR_RNO(RAND_VEL(IDIM, 1:MAX_PIP), MEAN_U, SIG_U)
548           enddo
549           !WRITE(*, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(:,3)), MAXVAL(FC(:,3))
550           !WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(:,3)), MAXVAL(FC(:,3))
551     
552           DO L = 1, MAX_PIP
553              DELETE_PART = .false.
554     ! pradeep skip ghost particles
555              if(pc.gt.pip) exit
556              if(is_nonexistent(l)) cycle
557              pc = pc+1
558              if(is_ghost(l) .or. is_entering_ghost(l) .or. is_exiting_ghost(l)) cycle
559     
560              DES_LOC_DEBUG = .FALSE.
561     
562              IF(.NOT.IS_ENTERING(L))THEN
563                 FC(:,L) = FC(:,L)/PMASS(L) + GRAV(:)
564              ELSE
565                 FC(:,L) = ZERO
566              ENDIF
567     
568              IF(DES_FIXED_BED) FC(:,L) = ZERO
569     
570              !DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
571              !By comparing the MFIX and equations in those papers,
572              !D_p = Beta/(EP_S*RHOP)
573              !F_gp in drag_fgs.f  = Beta*PVOL/EP_S
574              !Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
575              DP_BAR = F_gp(L)/(PMASS(L))
576              IF(.NOT.DES_CONTINUUM_COUPLED) DP_BAR = ZERO
577     
578              if(.not.MPPIC_PDRAG_IMPLICIT) DP_BAR = ZERO
579              M = PIJK(L,5)
580              IJK = PIJK(L,4)
581              IJK_OLD = IJK
582              PIJK_OLD(:) = PIJK(L,:)
583     
584              I = I_OF(IJK)
585              J = J_OF(IJK)
586              K = K_OF(IJK)
587     
588              DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L) + &
589              & FC(:,L)*DTSOLID)/(1.d0+DP_BAR*DTSOLID)
590     
591              IF(DES_FIXED_BED) DES_VEL_NEW(:,L) = ZERO
592     
593              !MPPIC_VPTAU(L,:) = DES_VEL_NEW(:,L)
594     
595              VELP_INT(:) = DES_VEL_NEW(:,L)
596     
597              MEANVEL(1) = DES_U_S(IJK_OLD,M)
598              MEANVEL(2) = DES_V_S(IJK_OLD,M)
599              IF(DO_K) MEANVEL(3) = DES_W_S(IJK_OLD,M)
600     
601              RAD_EFF = DES_RADIUS(L)
602              !RAD_EFF = (DES_STAT_WT(L)**(1.d0/3.d0))*DES_RADIUS(L)
603              IF(.not.DES_ONEWAY_COUPLED) then
604                 MEAN_FREE_PATH  = MAX(1.d0/(1.d0-EP_STAR), 1.d0/(1.D0-EP_G(IJK_OLD)))
605                 MEAN_FREE_PATH  = MEAN_FREE_PATH*RAD_EFF/THREEINTOSQRT2
606              endif
607     
608              DO IDIM = 1, merge(2,3,NO_K)
609                 !SIG_U = 0.05D0*MEANVEL(IDIM)
610                 !DES_VEL_NEW(IDIM,L) = DES_VEL_NEW(IDIM,L) + SIG_U*RAND_VEL(IDIM, L )
611                 !PART_TAUP = RO_Sol(L)*((2.d0*DES_RADIUS(L))**2.d0)/(18.d0* MU_G(IJK))
612                 SIG_U = 0.005D0      !*MEANVEL(IDIM)
613                 RAND_VEL(IDIM, L)  = SIG_U*RAND_VEL(IDIM, L)*DES_VEL_NEW(IDIM,L)
614                 IF(DES_FIXED_BED) RAND_VEL(IDIM,L) = ZERO
615                 !rand_vel(idim, L)  = sig_u*rand_vel(idim, L)/part_taup
616                 !rand_vel(idim, L)  = sig_u* mean_free_path*rand_vel(idim, L)/part_taup
617                 DES_VEL_NEW(idim,L) = DES_VEL_NEW(idim,L) + rand_vel(idim, L)
618              enddo
619     
620              IF(.not.DES_ONEWAY_COUPLED.and.(.not.des_fixed_bed)) CALL MPPIC_APPLY_PS_GRAD_PART(L)
621     
622              UPRIMEMOD = SQRT(DOT_PRODUCT(DES_VEL_NEW(1:,L), DES_VEL_NEW(1:,L)))
623     
624              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) then
625              !   DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L)/UPRIMEMOD)*MEAN_FREE_PATH/DTSOLID
626              !ENDIF
627     
628              IF(DES_FIXED_BED) THEN
629                 DD(:) = ZERO
630              ELSE
631                 DD(:) = DES_VEL_NEW(:,L)*DTSOLID !+ rand_vel(:, L)*dtsolid
632              ENDIF
633     
634              CALL PIC_FIND_NEW_CELL(L)
635     
636              IJK = PIJK(L,4)
637     
638              !IF((EP_G(IJK).LT.0.35.and.fluid_at(ijk)).or.(ep_g(ijk_old).lt.0.35)) then !.and.(ijk.ne.ijk_old)) then
639              !IF((EP_G(IJK).LT.EP_STAR.and.fluid_at(ijk)).and.(ijk.ne.ijk_old)) then
640              INSIDE_DOMAIN = .true.
641              INSIDE_DOMAIN = FLUID_AT(IJK)
642     
643              IF(CUT_CELL_AT(IJK)) THEN
644                 POS_Z = zero
645                 IF(DO_K) POS_Z = DES_POS_NEW(3,L)
646                 CALL GET_DEL_H_DES(IJK,'SCALAR', &
647                 & DES_POS_NEW(1,L),  DES_POS_NEW(2,L), &
648                 & POS_Z, &
649                 & DIST, NORM1, NORM2, NORM3, .true.)
650     
651                 IF(DIST.LE.ZERO) INSIDE_DOMAIN = .false.
652              ENDIF
653     
654              !IF((EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN)) then
655              !IF(EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN) then
656              IF(1.d0 - EP_G(IJK).GT. 1.3d0*(1.d0 - EP_STAR).and.INSIDE_DOMAIN.and.IJK.NE.IJK_OLD.and.OUTER_STABILITY_COND) then
657                 DD(:) = ZERO
658                 DES_VEL_NEW(:,L) = 0.8d0*DES_VEL_NEW(:,L)
659              ENDIF
660     
661              PIJK(L,:) = PIJK_OLD(:)
662     
663              DES_POS_NEW(:,L) = DES_POS_NEW(:,L) + DD(:)
664     
665              DIST = SQRT(DOT_PRODUCT(DD,DD))
666     
667              D_GRIDUNITS(1) = ABS(DD(1))/DX(PIJK(L,1))
668              D_GRIDUNITS(2) = ABS(DD(2))/DY(PIJK(L,2))
669              D_GRIDUNITS(3) = 1.d0
670              IF(DO_K) D_GRIDUNITS(3) = ABS(DD(3))/DZ(PIJK(L,3))
671     
672              DIST = SQRT(DOT_PRODUCT(DD,DD))
673              DTPIC_TMPX = (CFL_PIC*DX(PIJK(L,1)))/(ABS(DES_VEL_NEW(1,L))+SMALL_NUMBER)
674              DTPIC_TMPY = (CFL_PIC*DY(PIJK(L,2)))/(ABS(DES_VEL_NEW(2,L))+SMALL_NUMBER)
675              DTPIC_TMPZ = LARGE_NUMBER
676              IF(DO_K) DTPIC_TMPZ = (CFL_PIC*DZ(PIJK(L,3)))/(ABS(DES_VEL_NEW(3,L))+SMALL_NUMBER)
677     
678     
679     ! Check if the particle has moved a distance greater than or equal to grid spacing
680     ! if so, then exit
681     
682              DO IDIM = 1, merge(2,3,NO_K)
683                 IF(D_GRIDUNITS(IDIM).GT.ONE) THEN
684                    IF(DMP_LOG.OR.myPE.eq.pe_IO) THEN
685     
686                       WRITE(UNIT_LOG, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
687                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,L)
688                       IF (DO_OLD) WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(:,l)
689                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(:,L)
690                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC          = ', FC(:,L)
691     
692                       WRITE(*, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
693     
694                       WRITE(*, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,L)
695                       IF (DO_OLD) WRITE(*, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(:,l)
696                       WRITE(*, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(:,L)
697                       WRITE(*, '(A,2x,3(g17.8))') 'FC          = ', FC(:,L)
698                       read(*,*)
699                       DELETE_PART = .true.
700     
701                    ENDIF
702                    !CALL mfix_exit(myPE)
703                 END IF
704     
705              END DO
706     
707              IF(.not.DELETE_PART) then
708                 DTPIC_CFL = MIN(DTPIC_CFL, DTPIC_TMPX, DTPIC_TMPY, DTPIC_TMPZ)
709                 DTPIC_MIN_X = MIN(DTPIC_MIN_X, DTPIC_TMPX)
710                 DTPIC_MIN_Y = MIN(DTPIC_MIN_Y, DTPIC_TMPY)
711                 DTPIC_MIN_Z = MIN(DTPIC_MIN_Z, DTPIC_TMPZ)
712              ENDIF
713              FC(:,L) = ZERO
714     
715              IF(DELETE_PART) THEN
716                 CALL SET_NONEXISTENT(L)
717                 PIP_DEL_COUNT = PIP_DEL_COUNT + 1
718              ENDIF
719              IF (DES_LOC_DEBUG) WRITE(*,1001)
720           ENDDO
721     
722     
723           DEALLOCATE(RAND_VEL)
724     
725           CALL global_all_max(DTPIC_CFL)
726           PIP = PIP - PIP_DEL_COUNT
727     
728           LPIP_DEL_COUNT_ALL(:) = 0
729           LPIP_DEL_COUNT_ALL(mype) = PIP_DEL_COUNT
730           CALL GLOBAL_ALL_SUM(LPIP_DEL_COUNT_ALL)
731           IF((DMP_LOG).AND.SUM(LPIP_DEL_COUNT_ALL(:)).GT.0) THEN
732              IF(PRINT_DES_SCREEN) THEN
733                 WRITE(*,'(/,2x,A,2x,i10,/,A)') &
734                      'TOTAL NUMBER OF PARTS STEPPING MORE THAN ONE GRID SPACE = ', SUM(LPIP_DEL_COUNT_ALL(:)), &
735                      'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
736              ENDIF
737     
738              WRITE(UNIT_LOG,'(/,2x,A,2x,i10,/,A)') &
739                   'TOTAL NUMBER OF PARTS  STEPPING MORE THAN ONE GRID SPACEC= ', SUM(LPIP_DEL_COUNT_ALL(:)), &
740                   'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
741              !DO IPROC = 0, NUMPES-1
742              !   WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
743              !     'PARTICLES OUTSIDE DOMAIN (PIC)  ON PROC:', IPROC,' EQUAL TO', LPIP_DEL_COUNT_ALL(IPROC)
744              !ENDDO
745     
746           ENDIF
747     
748           DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
749     
750           RETURN
751           IF(DTSOLID.GT.DTPIC_MAX) THEN
752              !IF(DMP_LOG) WRITE(UNIT_LOG, 2001) DTSOLID
753              IF(myPE.eq.pe_IO) WRITE(*, 2004) DTSOLID
754           ELSEIF(DTSOLID.LT.DTPIC_MAX) THEN
755     
756              !IF(DMP_LOG) WRITE(UNIT_LOG, 2002) DTSOLID
757              IF(myPE.eq.pe_IO) WRITE(*, 2002) DTPIC_MAX
758           ELSE
759             !WRITE(*,'(A40,2x,g17.8)') 'DT
760             !IF(mype.eq.pe_IO) WRITE(*,2003) DTSOLID
761           END IF
762     
763           WRITE(UNIT_LOG, '(10x,A,2x,3(g17.8))') 'DTPIC MINS IN EACH DIRECTION = ', DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
764           WRITE(*, '(10x,A,2x,3(g17.8))') 'DTPIC MINS IN EACH DIRECTION = ', DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
765     
766      1001 FORMAT(3X,'<---------- END CFNEWVALUES ----------')
767     
768     2001  FORMAT(/1X,70('*'),//,10X,  &
769                & 'MOVEMENT UNDESIRED IN CFNEWVALUES: PARTICLE', i5, /,10X, &
770                & ' MOVED MORE THAN A GRID SPACING IN ONE TIME STEP', /,10X, &
771                & 'MOVEMENT IN GRID UNITS = ', 3(2x, g17.8),/,10X,  &
772                & 'TERMINAL ERROR: NOT STOPPING, BUT DELETING THE PARTICLE', &
773                & /1X,70('*'), /10X, &
774                & 'DES_VEL_NEW = ',  3(2x, g17.8))
775     
776      2004 FORMAT(/10x, &
777           & 'DTSOLID SHUD BE REDUCED TO', g17.8)
778     
779      2002 FORMAT(/10x, &
780           & 'DTSOLID CAN BE INCREASED TO', g17.8)
781     
782           RETURN
783           END SUBROUTINE CFNEWVALUES_MPPIC
784     
785     
786     !------------------------------------------------------------------------
787     ! subroutine       : des_dbgpic
788     ! Author           : Pradeep G.
789     ! Purpose          : For debugging the pic values
790     ! Parameters       : pstart - start indices of the particle
791     !                    pend - end indices of the particle
792     !                    pfreq - optional frequency (when the local count matches the
793     !                    frequency the filw will be written)
794     !                    if not send then it prints the file
795     !------------------------------------------------------------------------
796     
797           subroutine des_dbgpic (pstart,pend,pfreq)
798     
799     !-----------------------------------------------
800     ! Modules
801     !-----------------------------------------------
802           use discretelement
803           use fldvar
804           use functions
805           implicit none
806     !-----------------------------------------------
807     ! Dummy arguments
808     !-----------------------------------------------
809           INTEGER, INTENT(IN) :: pstart,pend
810           INTEGER, INTENT(IN), OPTIONAL :: pfreq
811     !-----------------------------------------------
812     ! Local variables
813     !-----------------------------------------------
814           integer lp,lijk
815           integer, save :: lfcount = 0 ,lfreq =0
816           character(255) :: filename
817     !-----------------------------------------------
818           if (present(pfreq)) then
819              lfreq = lfreq+1
820              if (lfreq .ne. pfreq) return
821              lfreq =0
822           end if
823           lfcount = lfcount + 1
824           write(filename,'("debug",I3.3)') lfcount
825           open (unit=100,file=filename)
826           do lp = pstart,pend
827              if (is_normal(lp) .or. is_entering(lp) .or. is_exiting(lp)) then
828                 lijk = pijk(lp,4)
829                 write(100,*)"positon =",lijk,pijk(lp,1),pijk(lp,2), &
830                    pijk(lp,3),ep_g(lijk),DES_U_s(lijk,1)
831                 write(100,*)"forces =", FC(2,lp),tow(1,lp)
832              end if
833           end do
834           close (100)
835     
836           RETURN
837           END SUBROUTINE des_dbgpic
838     
839     !------------------------------------------------------------------------
840     ! subroutine       : des_dbgtecplot
841     ! Author           : Pradeep G.
842     ! Purpose          : prints the tecplot file for particle location
843     ! Parameters       : pstart - start indices of the particle
844     !                    pend - end indices of the particle
845     !                    pfreq - optional frequency (when the local count matches the
846     !                    frequency the filw will be written)
847     !                    if not send then it prints the file
848     !------------------------------------------------------------------------
849     
850           subroutine des_dbgtecplot (pstart,pend,pfreq)
851     
852     !-----------------------------------------------
853     ! Modules
854     !-----------------------------------------------
855           USE discretelement
856           USE fldvar
857           USE functions
858           implicit none
859     !-----------------------------------------------
860     ! Dummy arguments
861     !-----------------------------------------------
862           INTEGER, INTENT(IN) :: pstart,pend
863           INTEGER, INTENT(IN), OPTIONAL :: pfreq
864     !-----------------------------------------------
865     ! Local variables
866     !-----------------------------------------------
867           integer lp,lijk
868           integer, save :: lfcount = 0 ,lfreq =0
869           character(255) :: filename
870     !-----------------------------------------------
871     
872           if (present(pfreq)) then
873              lfreq = lfreq+1
874              if (lfreq .ne. pfreq) return
875              lfreq =0
876           end if
877           lfcount = lfcount + 1
878           write(filename,'("new_tec",I3.3,".dat")') lfcount
879           open (unit=100,file=filename)
880           write(100,'(9(A,3X),A)') &
881              'VARIABLES = ', '"ijk"', '"x"', '"y"', '"vx"', '"vy"', &
882              '"ep_g"', '"FCX"' ,'"FCY"', '"TOW"'
883           write(100,'(A,F14.7,A)') 'zone T = "' , s_time , '"'
884           do lp = pstart,pend
885              if (.not.is_nonexistent(lp)) then
886                 lijk = pijk(lp,4)
887                 write(100,*)lijk,des_pos_new(1,lp),des_pos_new(2,lp), &
888                    des_vel_new(1,lp),des_vel_new(2,lp),ep_g(lijk),&
889                    fc(1,lp),fc(2,lp),tow(lp,1)
890              endif
891           enddo
892           close (100)
893           RETURN
894           END SUBROUTINE DES_DBGTECPLOT
895