File: /nfs/home/0/users/jenkins/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 param
25           USE param1
26           USE parallel
27           USE physprop
28           USE constant
29           USE fldvar
30           USE discretelement
31           USE des_bc
32           USE mpi_utility
33           USE mfix_pic
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(.NOT.PEA(L,1)) CYCLE  ! Only real particles
59                 IF(PEA(L,2)) CYCLE       ! Only non-entering
60                 IF(PEA(L,4)) 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,pea,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(.NOT.PEA(L,1)) CYCLE
78     ! skip ghost particles
79              IF(PEA(L,4)) 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.PEA(L,2))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           END SUBROUTINE CFNEWVALUES
172     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
173     !                                                                      C
174     !  Module name: CFNEWVALUES_MPPIC_SNIDER                               C
175     !
176     !                                                                      C
177     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
178           SUBROUTINE CFNEWVALUES_MPPIC_SNIDER
179     
180           USE param
181           USE param1
182           USE parallel
183           USE physprop
184           USE constant
185           USE fldvar
186           USE discretelement
187           USE des_bc
188           USE mpi_utility
189           USE fldvar
190           USE cutcell
191           USE mfix_pic
192           USE randomno
193           USE geometry, only: DO_K, NO_K
194           USE fun_avg
195           USE functions
196     
197           IMPLICIT NONE
198     !-----------------------------------------------
199     ! Local Variables
200     !-----------------------------------------------
201           INTEGER L, M, IDIM
202           INTEGER I, J, K, IJK, IJK_OLD, IJK2, IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
203     
204           DOUBLE PRECISION DD(3), DIST, &
205                            NEIGHBOR_SEARCH_DIST, DP_BAR, COEFF_EN, MEANVEL(3), D_GRIDUNITS(3)
206     
207           DOUBLE PRECISION DELUP(3), UPRIMETAU(3), UPRIMETAU_INT(3), MEAN_FREE_PATH, PS_FORCE(3)
208     ! index to track accounted for particles
209           INTEGER PC
210     
211     ! Logical for local debug warnings
212           LOGICAL DES_LOC_DEBUG
213     
214     ! maximum distance particles can move in MPPIC
215           DOUBLE PRECISION MAXDIST_PIC, UPRIMEMOD, UPRIMEMODNEW, signvel
216     
217     ! dt's in each direction  based on cfl_pic for the mppic case
218     
219           DOUBLE PRECISION DTPIC_TMPX, DTPIC_TMPY , DTPIC_TMPZ, THREEINTOSQRT2, RAD_EFF, MEANUS(3, MMAX)
220           DOUBLE PRECISION :: DPS_DXE, DPS_DXW, DPS_DYN, DPS_DYS, DPS_DZT, DPS_DZB
221     
222           LOGICAL :: DELETE_PART, INSIDE_DOMAIN
223           INTEGER :: PIP_DEL_COUNT
224     
225           DOUBLE PRECISION MEANUS_e(3, MMAX), MEANUS_w(3, MMAX),MEANUS_n(3, MMAX),MEANUS_s(3, MMAX),MEANUS_t(3, MMAX), MEANUS_b(3, MMAX)
226           INTEGER :: LPIP_DEL_COUNT_ALL(0:numPEs-1), PIJK_OLD(5)
227     
228     
229           double precision  sig_u, mean_u
230           double precision, allocatable, dimension(:,:) ::  rand_vel
231     !-----------------------------------------------
232     
233           PC = 1
234           DTPIC_CFL = LARGE_NUMBER
235     
236           if(NO_K) THREEINTOSQRT2 = 2.D0*SQRT(2.D0)
237           if(DO_K) THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
238           THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
239           DES_VEL_MAX(:) = ZERO
240           PIP_DEL_COUNT = 0
241     
242           !EPG_MIN2 = MINVAL(EP_G(:))
243           !epg_min_loc = MINLOC(EP_G)
244           !IJK = epg_min_loc(1)
245     
246           allocate(rand_vel(3, MAX_PIP))
247           do idim = 1, merge(2,3,NO_K)
248              mean_u = zero
249              sig_u = 1.d0
250              CALL NOR_RNO(RAND_VEL(IDIM, 1:MAX_PIP), MEAN_U, SIG_U)
251           enddo
252     
253           DO L = 1, MAX_PIP
254              DELETE_PART = .false.
255     ! pradeep skip ghost particles
256              if(pc.gt.pip) exit
257              if(.not.pea(l,1)) cycle
258              pc = pc+1
259              if(pea(l,4)) cycle
260     
261              DES_LOC_DEBUG = .FALSE.
262     
263              !if(mppic) then
264              !   IJK = PIJK(L, 4)
265     
266              !   COUNT_BC = DES_CELLWISE_BCDATA(IJK)%COUNT_DES_BC
267              !   IF(COUNT_BC.GE.1) CALL MPPIC_ADD_FRIC_FORCE(L)
268              !ENDIF
269     ! If a particle is classified as new, then forces are ignored.
270     ! Classification from new to existing is performed in routine
271     ! des_check_new_particle.f
272              IF(.NOT.PEA(L,2))THEN
273                 FC(:,L) = FC(:,L)/PMASS(L) + GRAV(:)
274              ELSE
275                 FC(:,L) = ZERO
276              ENDIF
277     
278              !DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
279              !By comparing MFIX equations and equations in those papers,
280              !D_p = Beta/(EP_S*RHOP)
281              !F_gp in drag_fgs.f  = Beta*PVOL/EP_S
282              !Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
283              DP_BAR = F_gp(L)/(PMASS(L))
284              IF(.NOT.DES_CONTINUUM_COUPLED) DP_BAR = ZERO
285     
286              if(.not.MPPIC_PDRAG_IMPLICIT) DP_BAR = ZERO
287     
288              M = PIJK(L,5)
289              IJK = PIJK(L,4)
290     
291              IJK_OLD = IJK
292     
293              PIJK_OLD(:) = PIJK(L,:)
294     
295              I = I_OF(IJK)
296              J = J_OF(IJK)
297              K = K_OF(IJK)
298              COEFF_EN =  MPPIC_COEFF_EN1
299              UPRIMETAU(:) = ZERO
300     
301              DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L) + &
302              & FC(:,L)*DTSOLID)/(1.d0+DP_BAR*DTSOLID)
303     
304              do idim = 1, merge(2,3,NO_K)
305                 SIG_U = 0.05D0
306                 rand_vel(idim, L)  = sig_u*DES_VEL_NEW(IDIM,L)*rand_vel(idim, L)
307              enddo
308     
309              !MEANVEL(1) = DES_U_S(IJK_OLD,M)
310              !MEANVEL(2) = DES_V_S(IJK_OLD,M)
311              !IF(DO_K) MEANVEL(3) = DES_W_S(IJK_OLD,M)
312     
313              PS_FORCE(:) = PS_GRAD(L, :)
314              DELUP(:) = -( DTSOLID*PS_FORCE(:))/((1.d0+DP_BAR*DTSOLID))
315              DELUP(:) = DELUP(:)/ROP_S(IJK_OLD,M)
316     
317              MEANVEL(:) = AVGSOLVEL_P(:,L)
318     
319              DO IDIM = 1, merge(2,3,NO_K)
320                 IF(PS_FORCE(IDIM).LE.ZERO) THEN
321                    UPRIMETAU(IDIM) = MIN(DELUP(IDIM), (1+COEFF_EN)*(MEANVEL(IDIM)-DES_VEL_NEW(IDIM,L)))
322                    UPRIMETAU_INT(IDIM) = UPRIMETAU(IDIM)
323                    UPRIMETAU(IDIM) = MAX(UPRIMETAU(IDIM), ZERO)
324                 ELSE
325                    UPRIMETAU(IDIM) = MAX(DELUP(IDIM), (1+COEFF_EN)*(MEANVEL(IDIM)-DES_VEL_NEW(IDIM,L)))
326                    UPRIMETAU_INT(IDIM) = UPRIMETAU(IDIM)
327                    UPRIMETAU(IDIM) = MIN(UPRIMETAU(IDIM), ZERO)
328                 END IF
329     
330              ENDDO
331     
332              DES_VEL_NEW(:,L) = DES_VEL_NEW(:,L) + UPRIMETAU(:)
333              DD(:) = DES_VEL_NEW(:,L)*DTSOLID
334     
335              UPRIMEMOD = SQRT(DOT_PRODUCT(DES_VEL_NEW(1:,L), DES_VEL_NEW(1:,L)))
336     
337              RAD_EFF = DES_RADIUS(L)
338                    !RAD_EFF = (DES_STAT_WT(L)**(1.d0/3.d0))*DES_RADIUS(L)
339              MEAN_FREE_PATH  = MAX(1.d0/(1.d0-EP_STAR), 1.d0/(1.D0-EP_G(IJK_OLD)))
340              MEAN_FREE_PATH = MEAN_FREE_PATH*RAD_EFF/THREEINTOSQRT2
341     
342              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) then
343              !   DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L)/UPRIMEMOD)*MEAN_FREE_PATH/DTSOLID
344              !ENDIF
345     
346               DIST = SQRT(DOT_PRODUCT(DD,DD))
347     
348              CALL PIC_FIND_NEW_CELL(L)
349     
350              IJK = PIJK(L,4)
351     
352              INSIDE_DOMAIN = .true.
353     
354              INSIDE_DOMAIN = FLUID_AT(IJK)
355     
356              IF(EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN.and.IJK.NE.IJK_OLD) then
357                 DD(:) = rand_vel(:, L)*dtsolid
358                 DES_VEL_NEW(:,L) = 0.8d0*DES_VEL_NEW(:,L)
359              ENDIF
360     
361     
362     
363              PIJK(L,:) = PIJK_OLD(:)
364              DIST = SQRT(DOT_PRODUCT(DD,DD))
365     
366               IF(DIST.GT.MEAN_FREE_PATH) THEN
367                  DD(:) =  DES_VEL_NEW(:,L)*DTSOLID*MEAN_FREE_PATH/DIST
368               ENDIF
369     
370              DES_POS_NEW(:,L) = DES_POS_NEW(:,L) + DD(:)
371     
372              D_GRIDUNITS(1) = ABS(DD(1))/DX(PIJK(L,1))
373              D_GRIDUNITS(2) = ABS(DD(2))/DY(PIJK(L,2))
374              D_GRIDUNITS(3) = 1.d0
375              IF(DO_K) D_GRIDUNITS(3) = ABS(DD(3))/DZ(PIJK(L,3))
376     
377              DIST = SQRT(DOT_PRODUCT(DD,DD))
378              DTPIC_TMPX = (CFL_PIC*DX(PIJK(L,1)))/(ABS(DES_VEL_NEW(1,L))+SMALL_NUMBER)
379              DTPIC_TMPY = (CFL_PIC*DY(PIJK(L,2)))/(ABS(DES_VEL_NEW(2,L))+SMALL_NUMBER)
380              DTPIC_TMPZ = LARGE_NUMBER
381              IF(DO_K) DTPIC_TMPZ = (CFL_PIC*DZ(PIJK(L,3)))/(ABS(DES_VEL_NEW(3,L))+SMALL_NUMBER)
382     
383     
384     ! Check if the particle has moved a distance greater than or equal to grid spacing
385     ! if so, then exit
386     
387              DO IDIM = 1, merge(2,3,NO_K)
388                 IF(D_GRIDUNITS(IDIM).GT.ONE) THEN
389                    IF(DMP_LOG.OR.myPE.eq.pe_IO) THEN
390     
391                       WRITE(UNIT_LOG, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
392                       WRITE(*, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
393                       DELETE_PART = .true.
394     
395                    ENDIF
396                    !CALL mfix_exit(myPE)
397                 END IF
398     
399              END DO
400              IF(.not.DELETE_PART) DTPIC_CFL = MIN(DTPIC_CFL, DTPIC_TMPX, DTPIC_TMPY, DTPIC_TMPZ)
401              FC(:,L) = ZERO
402     
403              IF(DELETE_PART) THEN
404                 PEA(L,1) = .false.
405                 PIP_DEL_COUNT = PIP_DEL_COUNT + 1
406              ENDIF
407              IF (DES_LOC_DEBUG) WRITE(*,1001)
408           ENDDO
409     
410     
411           IF(MPPIC) THEN
412              CALL global_all_max(DTPIC_CFL)
413              PIP = PIP - PIP_DEL_COUNT
414     
415              LPIP_DEL_COUNT_ALL(:) = 0
416              LPIP_DEL_COUNT_ALL(mype) = PIP_DEL_COUNT
417              CALL GLOBAL_ALL_SUM(LPIP_DEL_COUNT_ALL)
418              IF((DMP_LOG).AND.SUM(LPIP_DEL_COUNT_ALL(:)).GT.0) THEN
419                 IF(PRINT_DES_SCREEN) THEN
420                    WRITE(*,'(/,2x,A,2x,i10,/,A)') &
421                         'TOTAL NUMBER OF PARTS STEPPING MORE THAN ONE GRID SPACE = ', SUM(LPIP_DEL_COUNT_ALL(:)), &
422                         'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
423                 ENDIF
424     
425                 WRITE(UNIT_LOG,'(/,2x,A,2x,i10,/,A)') &
426                      'TOTAL NUMBER OF PARTS  STEPPING MORE THAN ONE GRID SPACEC= ', SUM(LPIP_DEL_COUNT_ALL(:)), &
427                      'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
428                 !DO IPROC = 0, NUMPES-1
429                 !   WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
430                 !   'PARTICLES OUTSIDE DOMAIN (PIC)  ON PROC:', IPROC,' EQUAL TO', LPIP_DEL_COUNT_ALL(IPROC)
431                 !ENDDO
432     
433              ENDIF
434     
435              DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
436     
437              IF(DTSOLID.GT.DTPIC_MAX) THEN
438                 !IF(DMP_LOG) WRITE(UNIT_LOG, 2001) DTSOLID
439                 IF(myPE.eq.pe_IO) WRITE(*, 2004) DTSOLID
440              ELSEIF(DTSOLID.LT.DTPIC_MAX) THEN
441     
442                 !IF(DMP_LOG) WRITE(UNIT_LOG, 2002) DTSOLID
443                 IF(myPE.eq.pe_IO) WRITE(*, 2002) DTPIC_MAX
444              ELSE
445                 !WRITE(*,'(A40,2x,g17.8)') 'DT
446                 !IF(mype.eq.pe_IO) WRITE(*,2003) DTSOLID
447              END IF
448     
449     
450           ENDIF
451      1000 FORMAT(3X,'---------- FROM CFNEWVALUES ---------->')
452      1001 FORMAT(3X,'<---------- END CFNEWVALUES ----------')
453     
454      1002 FORMAT(/1X,70('*')//&
455              ' From: CFNEWVALUES -',/&
456              ' Message: Particle ',I10, ' moved a distance ', ES17.9, &
457              ' during a',/10X, 'single solids time step, which is ',&
458              ' greater than',/10X,'its radius: ', ES17.9)
459      1003 FORMAT(1X,70('*')/)
460     
461     2001  FORMAT(/1X,70('*'),//,10X,  &
462                & 'MOVEMENT UNDESIRED IN CFNEWVALUES: PARTICLE', i5, /,10X, &
463                & ' MOVED MORE THAN A GRID SPACING IN ONE TIME STEP', /,10X, &
464                & 'MOVEMENT IN GRID UNITS = ', 3(2x, g17.8),/,10X,  &
465                & 'TERMINAL ERROR: NOT STOPPING, BUT DELETING THE PARTICLE', &
466                & /1X,70('*'), /10X, &
467                & 'DES_VEL_NEW = ',  3(2x, g17.8))
468     
469      2004 FORMAT(/10x, &
470           & 'DTSOLID SHUD BE REDUCED TO', g17.8)
471     
472      2002 FORMAT(/10x, &
473           & 'DTSOLID CAN BE INCREASED TO', g17.8)
474     
475      2003 FORMAT(/10x, &
476           & 'DTSOLID REMAINS UNCHANGED AT = ', g17.8)
477     
478     
479           RETURN
480           END SUBROUTINE CFNEWVALUES_MPPIC_SNIDER
481     
482     
483           SUBROUTINE CFNEWVALUES_MPPIC
484     
485           USE param
486           USE param1
487           USE parallel
488           USE physprop
489           USE constant
490           USE fldvar
491           USE discretelement
492           USE des_bc
493           USE mpi_utility
494           USE mfix_pic
495           USE randomno
496           USE cutcell
497           USE fun_avg
498           USE functions
499           IMPLICIT NONE
500     !-----------------------------------------------
501     ! Local Variables
502     !-----------------------------------------------
503           INTEGER L, M, IDIM
504           INTEGER I, J, K, IJK, IJK_OLD
505     
506           DOUBLE PRECISION DD(3), DIST, &
507                            DP_BAR, MEANVEL(3), D_GRIDUNITS(3)
508     
509           DOUBLE PRECISION MEAN_FREE_PATH
510     ! index to track accounted for particles
511           INTEGER PC
512     
513     ! Logical for local debug warnings
514           LOGICAL DES_LOC_DEBUG
515     
516     ! maximum distance particles can move in MPPIC
517           DOUBLE PRECISION UPRIMEMOD
518     
519     ! dt's in each direction  based on cfl_pic for the mppic case
520     
521           DOUBLE PRECISION DTPIC_TMPX, DTPIC_TMPY , DTPIC_TMPZ, THREEINTOSQRT2, RAD_EFF, POS_Z
522           DOUBLE PRECISION :: VELP_INT(3)
523     
524           LOGICAL :: DELETE_PART, INSIDE_DOMAIN
525           INTEGER :: PIP_DEL_COUNT
526     
527           INTEGER :: LPIP_DEL_COUNT_ALL(0:numPEs-1), PIJK_OLD(5)
528     
529           double precision  sig_u, mean_u, DTPIC_MIN_X,  DTPIC_MIN_Y,  DTPIC_MIN_Z
530           double precision, allocatable, dimension(:,:) ::  rand_vel
531     
532           double precision :: norm1, norm2, norm3
533           Logical :: OUTER_STABILITY_COND, DES_FIXED_BED
534     !-----------------------------------------------
535     
536           OUTER_STABILITY_COND = .false.
537           DES_FIXED_BED = .false.
538           PC = 1
539           DTPIC_CFL = LARGE_NUMBER
540           DTPIC_MIN_X =  LARGE_NUMBER
541           DTPIC_MIN_Y =  LARGE_NUMBER
542           DTPIC_MIN_Z =  LARGE_NUMBER
543     
544           if(NO_K) THREEINTOSQRT2 = 2.D0*SQRT(2.D0)
545           if(DO_K) THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
546           THREEINTOSQRT2 = 3.D0*SQRT(2.D0)
547           DES_VEL_MAX(:) = ZERO
548           PIP_DEL_COUNT = 0
549     
550           !EPG_MIN2 = MINVAL(EP_G(:))
551           !epg_min_loc = MINLOC(EP_G)
552           !IJK = epg_min_loc(1)
553           allocate(rand_vel(3, MAX_PIP))
554           do idim = 1, merge(2,3,NO_K)
555              mean_u = zero
556              sig_u = 1.d0
557              CALL NOR_RNO(RAND_VEL(IDIM, 1:MAX_PIP), MEAN_U, SIG_U)
558           enddo
559           !WRITE(*, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(:,3)), MAXVAL(FC(:,3))
560           !WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC MIN AND MAX IN Z = ', MINVAL(FC(:,3)), MAXVAL(FC(:,3))
561     
562           DO L = 1, MAX_PIP
563              DELETE_PART = .false.
564     ! pradeep skip ghost particles
565              if(pc.gt.pip) exit
566              if(.not.pea(l,1)) cycle
567              pc = pc+1
568              if(pea(l,4)) cycle
569     
570              DES_LOC_DEBUG = .FALSE.
571     
572              IF(.NOT.PEA(L,2))THEN
573                 FC(:,L) = FC(:,L)/PMASS(L) + GRAV(:)
574              ELSE
575                 FC(:,L) = ZERO
576              ENDIF
577     
578              IF(DES_FIXED_BED) FC(:,L) = ZERO
579     
580              !DP_BAR is the D_p in the Snider paper (JCP, vol 170, 523-549, 2001)
581              !By comparing the MFIX and equations in those papers,
582              !D_p = Beta/(EP_S*RHOP)
583              !F_gp in drag_fgs.f  = Beta*PVOL/EP_S
584              !Therefore, D_p = F_gp/(PVOL*RHOP) = F_gp/PMASS
585              DP_BAR = F_gp(L)/(PMASS(L))
586              IF(.NOT.DES_CONTINUUM_COUPLED) DP_BAR = ZERO
587     
588              if(.not.MPPIC_PDRAG_IMPLICIT) DP_BAR = ZERO
589              M = PIJK(L,5)
590              IJK = PIJK(L,4)
591              IJK_OLD = IJK
592              PIJK_OLD(:) = PIJK(L,:)
593     
594              I = I_OF(IJK)
595              J = J_OF(IJK)
596              K = K_OF(IJK)
597     
598              DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L) + &
599              & FC(:,L)*DTSOLID)/(1.d0+DP_BAR*DTSOLID)
600     
601              IF(DES_FIXED_BED) DES_VEL_NEW(:,L) = ZERO
602     
603              !MPPIC_VPTAU(L,:) = DES_VEL_NEW(:,L)
604     
605              VELP_INT(:) = DES_VEL_NEW(:,L)
606     
607              MEANVEL(1) = DES_U_S(IJK_OLD,M)
608              MEANVEL(2) = DES_V_S(IJK_OLD,M)
609              IF(DO_K) MEANVEL(3) = DES_W_S(IJK_OLD,M)
610     
611              RAD_EFF = DES_RADIUS(L)
612              !RAD_EFF = (DES_STAT_WT(L)**(1.d0/3.d0))*DES_RADIUS(L)
613              IF(.not.DES_ONEWAY_COUPLED) then
614                 MEAN_FREE_PATH  = MAX(1.d0/(1.d0-EP_STAR), 1.d0/(1.D0-EP_G(IJK_OLD)))
615                 MEAN_FREE_PATH  = MEAN_FREE_PATH*RAD_EFF/THREEINTOSQRT2
616              endif
617     
618              DO IDIM = 1, merge(2,3,NO_K)
619                 !SIG_U = 0.05D0*MEANVEL(IDIM)
620                 !DES_VEL_NEW(IDIM,L) = DES_VEL_NEW(IDIM,L) + SIG_U*RAND_VEL(IDIM, L )
621                 !PART_TAUP = RO_Sol(L)*((2.d0*DES_RADIUS(L))**2.d0)/(18.d0* MU_G(IJK))
622                 SIG_U = 0.005D0      !*MEANVEL(IDIM)
623                 RAND_VEL(IDIM, L)  = SIG_U*RAND_VEL(IDIM, L)*DES_VEL_NEW(IDIM,L)
624                 IF(DES_FIXED_BED) RAND_VEL(IDIM,L) = ZERO
625                 !rand_vel(idim, L)  = sig_u*rand_vel(idim, L)/part_taup
626                 !rand_vel(idim, L)  = sig_u* mean_free_path*rand_vel(idim, L)/part_taup
627                 DES_VEL_NEW(idim,L) = DES_VEL_NEW(idim,L) + rand_vel(idim, L)
628              enddo
629     
630              IF(.not.DES_ONEWAY_COUPLED.and.(.not.des_fixed_bed)) CALL MPPIC_APPLY_PS_GRAD_PART(L)
631     
632              UPRIMEMOD = SQRT(DOT_PRODUCT(DES_VEL_NEW(1:,L), DES_VEL_NEW(1:,L)))
633     
634              !IF(UPRIMEMOD*DTSOLID.GT.MEAN_FREE_PATH) then
635              !   DES_VEL_NEW(:,L) = (DES_VEL_NEW(:,L)/UPRIMEMOD)*MEAN_FREE_PATH/DTSOLID
636              !ENDIF
637     
638              IF(DES_FIXED_BED) THEN
639                 DD(:) = ZERO
640              ELSE
641                 DD(:) = DES_VEL_NEW(:,L)*DTSOLID !+ rand_vel(:, L)*dtsolid
642              ENDIF
643     
644              CALL PIC_FIND_NEW_CELL(L)
645     
646              IJK = PIJK(L,4)
647     
648              !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
649              !IF((EP_G(IJK).LT.EP_STAR.and.fluid_at(ijk)).and.(ijk.ne.ijk_old)) then
650              INSIDE_DOMAIN = .true.
651              INSIDE_DOMAIN = FLUID_AT(IJK)
652     
653              IF(CUT_CELL_AT(IJK)) THEN
654                 POS_Z = zero
655                 IF(DO_K) POS_Z = DES_POS_NEW(3,L)
656                 CALL GET_DEL_H_DES(IJK,'SCALAR', &
657                 & DES_POS_NEW(1,L),  DES_POS_NEW(2,L), &
658                 & POS_Z, &
659                 & DIST, NORM1, NORM2, NORM3, .true.)
660     
661                 IF(DIST.LE.ZERO) INSIDE_DOMAIN = .false.
662              ENDIF
663     
664              !IF((EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN)) then
665              !IF(EP_G(IJK).LT.EP_STAR.and.INSIDE_DOMAIN) then
666              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
667                 DD(:) = ZERO
668                 DES_VEL_NEW(:,L) = 0.8d0*DES_VEL_NEW(:,L)
669              ENDIF
670     
671              PIJK(L,:) = PIJK_OLD(:)
672     
673              DES_POS_NEW(:,L) = DES_POS_NEW(:,L) + DD(:)
674     
675              DIST = SQRT(DOT_PRODUCT(DD,DD))
676     
677              D_GRIDUNITS(1) = ABS(DD(1))/DX(PIJK(L,1))
678              D_GRIDUNITS(2) = ABS(DD(2))/DY(PIJK(L,2))
679              D_GRIDUNITS(3) = 1.d0
680              IF(DO_K) D_GRIDUNITS(3) = ABS(DD(3))/DZ(PIJK(L,3))
681     
682              DIST = SQRT(DOT_PRODUCT(DD,DD))
683              DTPIC_TMPX = (CFL_PIC*DX(PIJK(L,1)))/(ABS(DES_VEL_NEW(1,L))+SMALL_NUMBER)
684              DTPIC_TMPY = (CFL_PIC*DY(PIJK(L,2)))/(ABS(DES_VEL_NEW(2,L))+SMALL_NUMBER)
685              DTPIC_TMPZ = LARGE_NUMBER
686              IF(DO_K) DTPIC_TMPZ = (CFL_PIC*DZ(PIJK(L,3)))/(ABS(DES_VEL_NEW(3,L))+SMALL_NUMBER)
687     
688     
689     ! Check if the particle has moved a distance greater than or equal to grid spacing
690     ! if so, then exit
691     
692              DO IDIM = 1, merge(2,3,NO_K)
693                 IF(D_GRIDUNITS(IDIM).GT.ONE) THEN
694                    IF(DMP_LOG.OR.myPE.eq.pe_IO) THEN
695     
696                       WRITE(UNIT_LOG, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
697                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,L)
698                       IF (DO_OLD) WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(:,l)
699                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(:,L)
700                       WRITE(UNIT_LOG, '(A,2x,3(g17.8))') 'FC          = ', FC(:,L)
701     
702                       WRITE(*, 2001) L, D_GRIDUNITS(:), DES_VEL_NEW(:,L)
703     
704                       WRITE(*, '(A,2x,3(g17.8))') 'rand_vel = ', rand_vel(:,L)
705                       IF (DO_OLD) WRITE(*, '(A,2x,3(g17.8))') 'des_pos_old = ', des_pos_old(:,l)
706                       WRITE(*, '(A,2x,3(g17.8))') 'des_pos_new = ', des_pos_new(:,L)
707                       WRITE(*, '(A,2x,3(g17.8))') 'FC          = ', FC(:,L)
708                       read(*,*)
709                       DELETE_PART = .true.
710     
711                    ENDIF
712                    !CALL mfix_exit(myPE)
713                 END IF
714     
715              END DO
716     
717              IF(.not.DELETE_PART) then
718                 DTPIC_CFL = MIN(DTPIC_CFL, DTPIC_TMPX, DTPIC_TMPY, DTPIC_TMPZ)
719                 DTPIC_MIN_X = MIN(DTPIC_MIN_X, DTPIC_TMPX)
720                 DTPIC_MIN_Y = MIN(DTPIC_MIN_Y, DTPIC_TMPY)
721                 DTPIC_MIN_Z = MIN(DTPIC_MIN_Z, DTPIC_TMPZ)
722              ENDIF
723              FC(:,L) = ZERO
724     
725              IF(DELETE_PART) THEN
726                 PEA(L,1) = .false.
727                 PIP_DEL_COUNT = PIP_DEL_COUNT + 1
728              ENDIF
729              IF (DES_LOC_DEBUG) WRITE(*,1001)
730           ENDDO
731     
732     
733           DEALLOCATE(RAND_VEL)
734     
735           CALL global_all_max(DTPIC_CFL)
736           PIP = PIP - PIP_DEL_COUNT
737     
738           LPIP_DEL_COUNT_ALL(:) = 0
739           LPIP_DEL_COUNT_ALL(mype) = PIP_DEL_COUNT
740           CALL GLOBAL_ALL_SUM(LPIP_DEL_COUNT_ALL)
741           IF((DMP_LOG).AND.SUM(LPIP_DEL_COUNT_ALL(:)).GT.0) THEN
742              IF(PRINT_DES_SCREEN) THEN
743                 WRITE(*,'(/,2x,A,2x,i10,/,A)') &
744                      'TOTAL NUMBER OF PARTS STEPPING MORE THAN ONE GRID SPACE = ', SUM(LPIP_DEL_COUNT_ALL(:)), &
745                      'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
746              ENDIF
747     
748              WRITE(UNIT_LOG,'(/,2x,A,2x,i10,/,A)') &
749                   'TOTAL NUMBER OF PARTS  STEPPING MORE THAN ONE GRID SPACEC= ', SUM(LPIP_DEL_COUNT_ALL(:)), &
750                   'THIS SHOULD NOT HAPPEN FREQUENTLY: MONITOR THIS MESSAGE'
751              !DO IPROC = 0, NUMPES-1
752              !   WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
753              !     'PARTICLES OUTSIDE DOMAIN (PIC)  ON PROC:', IPROC,' EQUAL TO', LPIP_DEL_COUNT_ALL(IPROC)
754              !ENDDO
755     
756           ENDIF
757     
758           DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
759     
760           RETURN
761           IF(DTSOLID.GT.DTPIC_MAX) THEN
762              !IF(DMP_LOG) WRITE(UNIT_LOG, 2001) DTSOLID
763              IF(myPE.eq.pe_IO) WRITE(*, 2004) DTSOLID
764           ELSEIF(DTSOLID.LT.DTPIC_MAX) THEN
765     
766              !IF(DMP_LOG) WRITE(UNIT_LOG, 2002) DTSOLID
767              IF(myPE.eq.pe_IO) WRITE(*, 2002) DTPIC_MAX
768           ELSE
769             !WRITE(*,'(A40,2x,g17.8)') 'DT
770             !IF(mype.eq.pe_IO) WRITE(*,2003) DTSOLID
771           END IF
772     
773           WRITE(UNIT_LOG, '(10x,A,2x,3(g17.8))') 'DTPIC MINS IN EACH DIRECTION = ', DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
774           WRITE(*, '(10x,A,2x,3(g17.8))') 'DTPIC MINS IN EACH DIRECTION = ', DTPIC_MIN_X, DTPIC_MIN_Y, DTPIC_MIN_Z
775     
776     
777      1000 FORMAT(3X,'---------- FROM CFNEWVALUES ---------->')
778      1001 FORMAT(3X,'<---------- END CFNEWVALUES ----------')
779     
780      1002 FORMAT(/1X,70('*')//&
781              ' From: CFNEWVALUES -',/&
782              ' Message: Particle ',I10, ' moved a distance ', ES17.9, &
783              ' during a',/10X, 'single solids time step, which is ',&
784              ' greater than',/10X,'its radius: ', ES17.9)
785      1003 FORMAT(1X,70('*')/)
786     
787     2001  FORMAT(/1X,70('*'),//,10X,  &
788                & 'MOVEMENT UNDESIRED IN CFNEWVALUES: PARTICLE', i5, /,10X, &
789                & ' MOVED MORE THAN A GRID SPACING IN ONE TIME STEP', /,10X, &
790                & 'MOVEMENT IN GRID UNITS = ', 3(2x, g17.8),/,10X,  &
791                & 'TERMINAL ERROR: NOT STOPPING, BUT DELETING THE PARTICLE', &
792                & /1X,70('*'), /10X, &
793                & 'DES_VEL_NEW = ',  3(2x, g17.8))
794     
795      2004 FORMAT(/10x, &
796           & 'DTSOLID SHUD BE REDUCED TO', g17.8)
797     
798      2002 FORMAT(/10x, &
799           & 'DTSOLID CAN BE INCREASED TO', g17.8)
800     
801      2003 FORMAT(/10x, &
802           & 'DTSOLID REMAINS UNCHANGED AT = ', g17.8)
803     
804     
805           RETURN
806           END SUBROUTINE CFNEWVALUES_MPPIC
807     
808     
809     !------------------------------------------------------------------------
810     ! subroutine       : des_dbgpic
811     ! Author           : Pradeep G.
812     ! Purpose          : For debugging the pic values
813     ! Parameters       : pstart - start indices of the particle
814     !                    pend - end indices of the particle
815     !                    pfreq - optional frequency (when the local count matches the
816     !                    frequency the filw will be written)
817     !                    if not send then it prints the file
818     !------------------------------------------------------------------------
819     
820           subroutine des_dbgpic (pstart,pend,pfreq)
821     
822     !-----------------------------------------------
823     ! Modules
824     !-----------------------------------------------
825           use discretelement
826           USE fldvar
827           implicit none
828     !-----------------------------------------------
829     ! Dummy arguments
830     !-----------------------------------------------
831           INTEGER, INTENT(IN) :: pstart,pend
832           INTEGER, INTENT(IN), OPTIONAL :: pfreq
833     !-----------------------------------------------
834     ! Local variables
835     !-----------------------------------------------
836           integer lp,lijk
837           integer, save :: lfcount = 0 ,lfreq =0
838           character(30) :: filename
839     !-----------------------------------------------
840           if (present(pfreq)) then
841              lfreq = lfreq+1
842              if (lfreq .ne. pfreq) return
843              lfreq =0
844           end if
845           lfcount = lfcount + 1
846           write(filename,'("debug",I3.3)') lfcount
847           open (unit=100,file=filename)
848           do lp = pstart,pend
849              if (pea(lp,1) .and. .not.pea(lp,4)) then
850                 lijk = pijk(lp,4)
851                 write(100,*)"positon =",lijk,pijk(lp,1),pijk(lp,2), &
852                    pijk(lp,3),ep_g(lijk),DES_U_s(lijk,1)
853                 write(100,*)"forces =", FC(2,lp),tow(1,lp)
854              end if
855           end do
856           close (100)
857     
858           RETURN
859           END SUBROUTINE des_dbgpic
860     
861     
862     !------------------------------------------------------------------------
863     ! subroutine       : des_dbgtecplot
864     ! Author           : Pradeep G.
865     ! Purpose          : prints the tecplot file for particle location
866     ! Parameters       : pstart - start indices of the particle
867     !                    pend - end indices of the particle
868     !                    pfreq - optional frequency (when the local count matches the
869     !                    frequency the filw will be written)
870     !                    if not send then it prints the file
871     !------------------------------------------------------------------------
872     
873           subroutine des_dbgtecplot (pstart,pend,pfreq)
874     
875     !-----------------------------------------------
876     ! Modules
877     !-----------------------------------------------
878           use discretelement
879           USE fldvar
880           implicit none
881     !-----------------------------------------------
882     ! Dummy arguments
883     !-----------------------------------------------
884           INTEGER, INTENT(IN) :: pstart,pend
885           INTEGER, INTENT(IN), OPTIONAL :: pfreq
886     !-----------------------------------------------
887     ! Local variables
888     !-----------------------------------------------
889           integer lp,lijk
890           integer, save :: lfcount = 0 ,lfreq =0
891           character(30) :: filename
892     !-----------------------------------------------
893     
894           if (present(pfreq)) then
895              lfreq = lfreq+1
896              if (lfreq .ne. pfreq) return
897              lfreq =0
898           end if
899           lfcount = lfcount + 1
900           write(filename,'("new_tec",I3.3,".dat")') lfcount
901           open (unit=100,file=filename)
902           write(100,'(9(A,3X),A)') &
903              'VARIABLES = ', '"ijk"', '"x"', '"y"', '"vx"', '"vy"', &
904              '"ep_g"', '"FCX"' ,'"FCY"', '"TOW"'
905           write(100,'(A,F14.7,A)') 'zone T = "' , s_time , '"'
906           do lp = pstart,pend
907              if (pea(lp,1)) then
908                 lijk = pijk(lp,4)
909                 write(100,*)lijk,des_pos_new(1,lp),des_pos_new(2,lp), &
910                    des_vel_new(1,lp),des_vel_new(2,lp),ep_g(lijk),&
911                    fc(1,lp),fc(2,lp),tow(lp,1)
912              endif
913           enddo
914           close (100)
915           RETURN
916           END SUBROUTINE DES_DBGTECPLOT
917     
918     
919     
920