File: /nfs/home/0/users/jenkins/mfix.git/model/des/pic_routines.f
1
2
3
4
5
6
7 SUBROUTINE PIC_TIME_MARCH
8
9 USE param
10 USE param1
11 USE run
12 USE output
13 USE funits
14 USE discretelement
15 USE constant
16 USE sendrecv
17 USE des_bc
18 USE cutcell
19 USE mfix_pic
20 USE pic_bc
21 USE error_manager
22 USE fldvar, only: P_g
23
24 use mpi_utility, only: GLOBAL_ALL_SUM
25 use mpi_funs_des, only: DES_PAR_EXCHANGE
26
27 USE fun_avg
28 USE functions
29
30 IMPLICIT NONE
31
32
33
34 INTEGER IJK, I, J, K
35 INTEGER NN, FACTOR, NP, BCV_I, LL, PC, TIME_LOOP_COUNT
36
37
38
39 DOUBLE PRECISION DES_RES_TIME, DES_SPX_TIME
40
41
42
43 DOUBLE PRECISION TMP_DTS, DTSOLID_TMP
44 CHARACTER(LEN=5) :: FILENAME
45
46
47
48 LOGICAL,SAVE:: FIRST_PASS = .TRUE.
49
50
51 DOUBLE PRECISION xpos,ypos, zpos, NORM_CF(3), DIST
52
53
54 DOUBLE PRECISION :: DES_KE_VEC(DIMN)
55
56
57 double precision :: TEND_PIC_LOOP
58
59 Integer :: PIC_ITERS
60
61 LOGICAL FOCUS
62
63 INTEGER :: gPIP
64
65
66
67
68 CALL INIT_ERR_MSG("PIC_TIME_MARCH")
69
70 S_TIME = TIME
71 TIME_LOOP_COUNT = 0
72
73
74
75
76 IF(DES_CONTINUUM_COUPLED) CALL CALC_PG_GRAD
77
78 TEND_PIC_LOOP = MERGE(TIME+DT, TSTOP, DES_CONTINUUM_COUPLED)
79 PIC_ITERS = 0
80
81
82
83 DO WHILE(S_TIME.LT.TEND_PIC_LOOP)
84
85
86 = MERGE(MIN(DTPIC_MAX, DT), DTPIC_MAX, DES_CONTINUUM_COUPLED)
87 DTSOLID_ORIG = DTSOLID
88 IF(MOD(PIC_ITERS, 10).eq.0) then
89 IF(DES_CONTINUUM_COUPLED) then
90 WRITE(ERR_MSG, 2000) DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
91 ELSE
92 WRITE(ERR_MSG, 2001) S_TIME, DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
93 ENDIF
94 CALL FLUSH_ERR_MSG(HEADER = .FALSE., FOOTER = .FALSE.)
95 ENDIF
96
97 2000 FORMAT(/5x,'DTSOLID CURRENT = ',g17.8,/5x,'DTPIC_CFL',8x,'= ', &
98 g17.8, /5x,'DTPIC TAUP',7x,'= ',g17.8,/5x,'DT FLOW',10x,'= ', &
99 g17.8)
100
101 2001 FORMAT(/5x,'TIME',13X,'= ',g17.8,/5x,'DTSOLID CURRENT = ',g17.8,&
102 /5x,'DTPIC_CFL',8X,'= ', g17.8,/5x,'DTPIC TAUP',7x,'= ',g17.8,&
103 /5x,'DT FLOW',10X,'= ', g17.8)
104
105
106 PIC_ITERS = PIC_ITERS + 1
107
108
109
110
111
112 IF(S_TIME + DTSOLID > TEND_PIC_LOOP) THEN
113 WRITE(ERR_MSG, 2010) DTSOLID, TIME + DT - S_TIME
114 CALL FLUSH_ERR_MSG(HEADER = .FALSE., FOOTER = .FALSE.)
115 DTSOLID = TEND_PIC_LOOP - S_TIME
116 ENDIF
117
118 2010 FORMAT(/5X,'REDUCING DTSOLID TO ENSURE STIME + DTSOLID LE ', &
119 'TIME + DT', /5x,'DTSOLID ORIG = ', g17.8, /5x, &
120 'DTSOLID ACTUAL = ', g17.8)
121
122
123 CALL DES_PAR_EXCHANGE
124
125 CALL PARTICLES_IN_CELL
126
127 CALL COMP_MEAN_FIELDS
128
129
130
131 CALL REPORT_PIC_STATS
132
133 CALL MPPIC_COMPUTE_PS_GRAD
134 IF(DES_CONTINUUM_COUPLED) CALL CALC_DRAG_DES
135 IF (DO_OLD) CALL CFUPDATEOLD
136
137 CALL CFNEWVALUES
138
139
140
141 CALL PIC_APPLY_WALLBC_STL
142
143
144
145 = S_TIME + DTSOLID
146
147
148 DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
149
150
151
152
153 IF(.NOT.DES_CONTINUUM_COUPLED) THEN
154
155 = S_TIME
156
157
158
159 IF(PRINT_DES_DATA) THEN
160 IF ( (S_TIME+0.1d0*DTSOLID >= DES_SPX_TIME) .OR. &
161 (S_TIME+0.1d0*DTSOLID >= TSTOP)) then
162 DES_SPX_TIME = &
163 ( INT((S_TIME+0.1d0*DTSOLID)/DES_SPX_DT) &
164 + 1 )*DES_SPX_DT
165 CALL WRITE_DES_DATA
166 IF(DMP_LOG) WRITE(UNIT_LOG,'(3X,A,1X,ES15.5)') &
167 'DES data file written at time =', S_TIME
168 ENDIF
169 ENDIF
170
171 IF ( (S_TIME+0.1d0*DTSOLID >= DES_RES_TIME) .OR. &
172 (S_TIME+0.1d0*DTSOLID >= TSTOP) .OR. &
173 (NN == FACTOR) ) THEN
174 DES_RES_TIME = &
175 ( INT((S_TIME+0.1d0*DTSOLID)/DES_RES_DT) &
176 + 1 )*DES_RES_DT
177 CALL WRITE_RES0_DES
178
179
180 CALL WRITE_RES1
181 IF(DMP_LOG) WRITE(UNIT_LOG,'(3X,A,1X,ES15.5)') &
182 'DES.RES and .RES files written at time =', S_TIME
183 ENDIF
184 ENDIF
185
186
187
188
189 ENDDO
190
191 CALL GLOBAL_ALL_SUM(PIP, gPIP)
192 WRITE(ERR_MSG, 3000) trim(iVal(PIC_ITERS)), trim(iVal(gPIP))
193 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
194
195 3000 FORMAT(/'PIC NITs: ',A,3x,'Total PIP: ', A)
196
197
198
199
200
201
202
203
204 IF(.NOT.DES_CONTINUUM_COUPLED)then
205 if(dmp_log)write(unit_log,'(1X,A)')&
206 '<---------- END MPPIC_TIME_MARCH ----------'
207
208 else
209
210
211
212
213
214
215
216
217 END IF
218
219
220 CALL FINL_ERR_MSG
221 end SUBROUTINE PIC_TIME_MARCH
222
223
224
225
226
227
228
229 SUBROUTINE MPPIC_COMP_EULERIAN_VELS_NON_CG
230 USE compar
231 USE constant
232 USE discretelement
233 USE functions
234 USE geometry
235 USE indices
236 USE mfix_pic
237 USE parallel
238 USE param
239 USE param1
240 use desmpi
241 IMPLICIT NONE
242
243
244
245
246 LOGICAL :: FOCUS
247
248 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,&
249 IPJPK, IPJKP, IJPKP, IPJPKP, II, JJ, KK
250
251
252 INTEGER :: M
253
254 DO IJK = ijkstart3, ijkend3
255 I = I_OF(IJK)
256 J = J_OF(IJK)
257 K = K_OF(IJK)
258
259
260
261 (IJK, :) = ZERO
262 PIC_V_S(IJK, :) = ZERO
263 PIC_W_S(IJK, :) = ZERO
264 IF (.NOT.IS_ON_myPE_plus1layer(I,J,K)) CYCLE
265 IF(.NOT.FLUID_AT(IJK)) CYCLE
266
267 IF(I.GE.IMIN1.AND.I.LT.IMAX1) then
268
269 = IP_OF(IJK)
270 DO M = 1, DES_MMAX
271 PIC_U_S(IJK,M) = 0.5d0*(DES_U_S(IJK, M) + DES_U_S(IPJK,M))
272 ENDDO
273 ENDIF
274
275 if(J.GE.JMIN1.AND.J.LT.JMAX1) then
276
277 = JP_OF(IJK)
278 DO M = 1, DES_MMAX
279 PIC_V_S(IJK,M) = 0.5d0*(DES_V_S(IJK, M) + DES_V_S(IJPK,M))
280 ENDDO
281 ENDIF
282
283
284 if(K.GE.KMIN1.AND.K.LT.KMAX1.AND.DO_K) then
285
286 = KP_OF(IJK)
287 DO M = 1, DES_MMAX
288 PIC_W_S(IJK,M) = 0.5d0*(DES_W_S(IJK, M) + DES_W_S(IJKP,M))
289 ENDDO
290 ENDIF
291 ENDDO
292
293
294
295
296
297 CALL MPPIC_BC_U_S
298 CALL MPPIC_BC_V_S
299 IF(DO_K) CALL MPPIC_BC_W_S
300
301 END SUBROUTINE MPPIC_COMP_EULERIAN_VELS_NON_CG
302
303
304
305
306
307
308
309
310
311 SUBROUTINE MPPIC_COMP_EULERIAN_VELS_CG
312 USE compar
313 USE constant
314 USE cutcell
315 USE desmpi
316 USE discretelement
317 USE functions
318 USE geometry
319 USE indices
320 USE mfix_pic
321 USE parallel
322 USE param
323 USE param1
324 IMPLICIT NONE
325
326
327
328
329 LOGICAL :: FOCUS
330
331 INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,&
332 IPJPK, IPJKP, IJPKP, IPJPKP, II, JJ, KK
333
334
335 INTEGER :: M
336
337 DO IJK = ijkstart3, ijkend3
338 I = I_OF(IJK)
339 J = J_OF(IJK)
340 K = K_OF(IJK)
341
342
343
344 (IJK, :) = ZERO
345 PIC_V_S(IJK, :) = ZERO
346 PIC_W_S(IJK, :) = ZERO
347
348 IF (.NOT.IS_ON_myPE_plus1layer(I,J,K)) CYCLE
349
350 IF(WALL_U_AT(IJK)) THEN
351 PIC_U_S(IJK, :) = ZERO
352
353
354
355
356 ELSE
357 if(.not.FLUID_AT(IJK)) cycle
358
359 IPJK = IP_OF(IJK)
360 IF(FLUID_AT(IPJK)) THEN
361 DO M = 1, DES_MMAX
362 PIC_U_S(IJK,M) = 0.5d0*(DES_U_S(IJK, M) + DES_U_S(IPJK,M))
363 ENDDO
364 ELSE
365 PIC_U_S(IJK,:) = DES_U_S(IJK, :)
366 ENDIF
367 ENDIF
368
369 IF(WALL_V_AT(IJK)) THEN
370 PIC_V_S(IJK, :) = ZERO
371 ELSE
372 if(.not.FLUID_AT(IJK)) cycle
373 IJPK = JP_OF(IJK)
374 IF(FLUID_AT(IJPK)) THEN
375 DO M = 1, DES_MMAX
376 PIC_V_S(IJK,M) = 0.5d0*(DES_V_S(IJK, M) + DES_V_S(IJPK,M))
377 ENDDO
378 ELSE
379 PIC_V_S(IJK,:) = DES_V_S(IJK, :)
380 ENDIF
381 ENDIF
382
383 IF(DO_K) THEN
384 IF(WALL_W_AT(IJK)) THEN
385 PIC_W_S(IJK, :) = ZERO
386 ELSE
387 if(.not.FLUID_AT(IJK)) cycle
388 IJKP = KP_OF(IJK)
389 IF(FLUID_AT(IJKP)) THEN
390 DO M = 1, DES_MMAX
391 PIC_W_S(IJK,M) = 0.5d0*(DES_W_S(IJK, M) + DES_W_S(IJKP,M))
392 ENDDO
393 ELSE
394 PIC_W_S(IJK,:) = DES_W_S(IJK, :)
395 ENDIF
396 ENDIF
397 ENDIF
398 ENDDO
399
400 END SUBROUTINE MPPIC_COMP_EULERIAN_VELS_CG
401
402
403
404
405
406
407
408 SUBROUTINE MPPIC_APPLY_PS_GRAD_PART(L)
409
410 USE param
411 USE param1
412 USE parallel
413 USE constant
414 USE discretelement
415 USE mpi_utility
416 USE mfix_pic
417 USE cutcell
418 USE fldvar, only: ep_g
419 USE fun_avg
420 USE functions
421 IMPLICIT NONE
422 INTEGER, INTENT(IN) :: L
423
424
425
426 INTEGER M, IDIM
427 INTEGER I, J, K, IJK, IJK_C
428
429 DOUBLE PRECISION D(DIMN), DIST, DP_BAR, COEFF_EN, COEFF_EN2
430
431 DOUBLE PRECISION DELUP(DIMN), UPRIMETAU(DIMN), UPRIMETAU_INT(DIMN), PS_FORCE(DIMN), VEL_ORIG(DIMN)
432
433 INTEGER PC , epg_min_loc(1)
434
435
436 LOGICAL DES_LOC_DEBUG
437
438
439 DOUBLE PRECISION MAXDIST_PIC, UPRIMEMOD, UPRIMEMODNEW, signvel
440
441
442
443 DOUBLE PRECISION MEANUS(DIMN, DES_MMAX), RELVEL(DIMN), MEANVEL(DIMN), VEL_NEW(DIMN)
444
445
446 LOGICAL :: INSIDE_DOMAIN
447
448
449 = PIJK(L,5)
450 IJK = PIJK(L,4)
451 COEFF_EN = MPPIC_COEFF_EN1
452 COEFF_EN2 = MPPIC_COEFF_EN2
453
454 VEL_ORIG(1:DIMN) = DES_VEL_NEW(:,L)
455 VEL_NEW (1:DIMN) = DES_VEL_NEW(:,L)
456
457 IF(L.EQ.FOCUS_PARTICLE) THEN
458
459 WRITE(*,'(A20,2x,3(2x,i4))') 'CELL ID = ', PIJK(L,1:3)
460 WRITE(*,'(A20,2x,3(2x,g17.8))') 'EPS = ', 1.d0 - EP_g(PIJK(L,4))
461 WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', DES_VEL_NEW(:,L)
462
463 WRITE(*,'(A20,2x,3(2x,g17.8))') 'FC = ', FC(:,L)
464 ENDIF
465
466 MEANVEL(1) = DES_U_S(IJK,M)
467 MEANVEL(2) = DES_V_S(IJK,M)
468 IF(DO_K) MEANVEL(3) = DES_W_S(IJK,M)
469
470 PS_FORCE(:) = PS_GRAD(L, :)
471
472 (:) = -PS_FORCE(:)
473
474 MEANUS(:,M) = AVGSOLVEL_P (:,L)
475
476 (:) = DES_VEL_NEW(:,L) - MEANUS(:,M)
477
478
479
480 DO IDIM = 1, DIMN
481
482 IF(ABS(PS_FORCE(IDIM)).eq.zero) cycle
483
484 IF(VEL_ORIG(IDIM)*MEANUS(IDIM,M).GT.ZERO) THEN
485
486 IF(VEL_ORIG(IDIM)*DELUP(IDIM).GT.ZERO) THEN
487
488 IF(ABS(MEANUS(IDIM,M)) .GT. ABS(VEL_ORIG(IDIM))) THEN
489 IJK_C = IJK
490 IF(IDIM.eq.1) then
491 if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
492 if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
493 ELSEIF(IDIM.eq.2) then
494 if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
495 if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
496 ELSEIF(IDIM.eq.3) then
497 if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
498 if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
499 ENDIF
500 INSIDE_DOMAIN = .false.
501 INSIDE_DOMAIN = fluid_at(IJK_C)
502
503 if(INSIDE_DOMAIN) then
504 VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
505 endif
506
507 ELSE
508
509 ENDIF
510 ELSE
511 IF(ABS(VEL_ORIG(IDIM)).GT.ABS(MEANUS(IDIM,M))) then
512
513 = IJK
514 IF(IDIM.eq.1) then
515 if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
516 if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
517 ELSEIF(IDIM.eq.2) then
518 if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
519 if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
520 ELSEIF(IDIM.eq.3) then
521 if(VEL_ORIG(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
522 if(VEL_ORIG(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
523 ENDIF
524
525
526 = .false.
527 INSIDE_DOMAIN = fluid_at(IJK_C)
528
529 if(INSIDE_DOMAIN) then
530 VEL_NEW(IDIM) = (MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM))
531 else
532 VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
533 endif
534
535 IJK_C = IJK
536 IF(IDIM.eq.1) then
537 if(VEL_NEW(IDIM).LT.ZERO) IJK_C = IM_OF(IJK)
538 if(VEL_NEW(IDIM).GT.ZERO) IJK_C = IP_OF(IJK)
539 ELSEIF(IDIM.eq.2) then
540 if(VEL_NEW(IDIM).LT.ZERO) IJK_C = JM_OF(IJK)
541 if(VEL_NEW(IDIM).GT.ZERO) IJK_C = JP_OF(IJK)
542 ELSEIF(IDIM.eq.3) then
543 if(VEL_NEW(IDIM).LT.ZERO) IJK_C = KM_OF(IJK)
544 if(VEL_NEW(IDIM).GT.ZERO) IJK_C = KP_OF(IJK)
545 ENDIF
546 INSIDE_DOMAIN = .false.
547 INSIDE_DOMAIN = fluid_at(IJK_C)
548
549
550
551
552
553
554 ELSE
555
556 (IDIM) = COEFF_EN2 * VEL_ORIG(IDIM)
557
558
559
560 ENDIF
561 ENDIF
562 ELSE
563 IF(MEANUS(IDIM,M)*DELUP(IDIM).GT.ZERO) THEN
564 VEL_NEW(IDIM) = MEANUS(IDIM,M) - COEFF_EN*RELVEL(IDIM)
565
566 ELSE
567
568
569 ENDIF
570 ENDIF
571
572 IF(MPPIC_GRAV_TREATMENT) THEN
573 IF(DELUP(IDIM)*GRAV(IDIM).LT.ZERO.AND.VEL_ORIG(IDIM)*GRAV(IDIM).GT.ZERO) THEN
574 VEL_NEW(IDIM) = -COEFF_EN*VEL_ORIG(IDIM)
575
576 ENDIF
577 ENDIF
578
579 (IDIM, L) = VEL_NEW(IDIM)
580 ENDDO
581
582
583 if(L.eq.FOCUS_PARTICLE) THEN
584
585
586
587 WRITE(*,'(A20,2x,3(2x,i5))') 'PIJK I, J, K =', I_OF(IJK),J_OF(IJK),K_OF(IJK)
588 WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL ORIG = ', VEL_ORIG(:)
589 WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_VEL NEW = ', DES_VEL_NEW(:,L)
590 WRITE(*,'(A20,2x,3(2x,g17.8))') 'MEANUS = ', MEANUS(:,1)
591 WRITE(*,'(A20,2x,3(2x,g17.8))') 'DES_POS_NEW = ', DES_POS_NEW(:,L)
592 WRITE(*,'(A20,2x,3(2x,g17.8))') 'GRAD PS = ', PS_FORCE(:)
593 WRITE(*,'(A20,2x,3(2x,g17.8))') 'DELUP = ', DELUP(:)
594
595
596
597
598 read(*,*)
599 ENDIF
600
601
602
603
604
605
606
607 RETURN
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625 END SUBROUTINE MPPIC_APPLY_PS_GRAD_PART
626
627
628
629
630
631
632
633
634 SUBROUTINE MPPIC_COMPUTE_PS_GRAD
635 USE param
636 USE param1
637 USE parallel
638 USE physprop
639 USE fldvar, only: ep_g, u_g, v_g, w_g
640 USE run
641 USE geometry
642 USE indices
643 USE bc
644 USE compar
645 USE sendrecv
646 USE discretelement
647 USE constant
648 USE cutcell
649 USE interpolation
650 USE mfix_pic
651 USE fun_avg
652 USE functions
653 implicit none
654
655
656 INTEGER I, J, K, IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM, IDIM, M
657
658
659 DOUBLE PRECISION TEMP1, TEMP2, avg_factor, VOL_TOT_VEC(3), VOL_TOT_SCAL
660
661 integer :: korder, iw,ie,js,jn,kb,ktp, onew, pcell(3), cur_ijk, NP, nindx
662
663 integer :: ii,jj,kk, ipjpk, ijpkp, ipjkp, ipjpkp, I1, J1, K1
664
665 double precision :: vol_ijk, vol_ipjk, vol_ijpk, vol_ipjpk
666 double precision :: vol_ijkp, vol_ipjkp, vol_ijpkp, vol_ipjpkp
667
668 if(MPPIC_SOLID_STRESS_SNIDER) then
669
670 DO IJK = IJKSTART3, IJKEND3
671 IF(FLUID_AT(IJK)) THEN
672 PIC_P_S(IJK,1) = PSFAC_FRIC_PIC * ((1.d0 - EP_G(IJK))**FRIC_EXP_PIC)/ MAX(EP_G(IJK) &
673 - EP_STAR, FRIC_NON_SING_FAC*EP_G(IJK))
674
675 ELSE
676
677 (IJK,1) = PSFAC_FRIC_PIC * ((1.d0 - 0.1d0)**FRIC_EXP_PIC)/ MAX(0.1d0 - EP_STAR, FRIC_NON_SING_FAC*0.1d0)
678
679 ENDIF
680 ENDDO
681
682
683 ELSE
684 DO IJK = IJKSTART3, IJKEND3
685 PS_FORCE_PIC(IJK,:) = ZERO
686
687 IF(FLUID_AT(IJK)) THEN
688 if(EP_G(IJK).lt.ep_star) then
689 PIC_P_S(IJK,1) = one*(one-ep_g(ijk))
690 else
691
692 PIC_P_S(IJK,1) = ZERO
693 endif
694
695 ELSE
696 PIC_P_S(IJK,1) = 1.
697 ENDIF
698 ENDDO
699
700 ENDIF
701
702 CALL SEND_RECV(PIC_P_S,1)
703
704
705
706
707
708
709
710 DO IJK = IJKSTART3, IJKEND3
711 PS_FORCE_PIC(IJK,:) = ZERO
712 IF(.NOT.FLUID_AT(IJK)) CYCLE
713
714 I = I_OF(IJK)
715 J = J_OF(IJK)
716 K = K_OF(IJK)
717 IPJK = IP_OF(IJK)
718 IJPK = JP_OF(IJK)
719 IJKP = KP_OF(IJK)
720
721 if(fluid_at(ipjk)) then
722 PS_FORCE_PIC(IJK,1) = 2.d0*(PIC_P_S(IPJK,1) - PIC_P_S(IJK,1))/(DX(I)+DX(I_of(ipjk)))
723 else
724 if(PIC_P_S(IJK,1).GT.ZERO) then
725
726 (IJK,1) = 2.d0*(PIC_P_S(IPJK,1) - PIC_P_S(IJK,1))/(DX(I)+DX(I_of(ipjk)))
727 ELSE
728 PS_FORCE_PIC(IJK,1) = zero
729 endif
730 ENDIF
731
732 IF(FLUID_AT(IJPK)) THEN
733 PS_FORCE_PIC(IJK,2) = 2.d0*(PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/(DY(j)+Dy(j_of(ijpk)))
734 ELSE
735 IF(PIC_P_S(IJK,1).GT.ZERO) then
736 PS_FORCE_PIC(IJK,2) = 2.d0*(PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/(DY(j)+Dy(j_of(ijpk)))
737 ELSE
738 PS_FORCE_PIC(IJK,2) = zero
739 ENDIF
740 ENDIF
741
742 IF(DO_K) THEN
743 IF(FLUID_AT(IJKP)) then
744 PS_FORCE_PIC(IJK,3) = 2.d0*(PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/(Dz(k)+Dz(k_of(ijkp)))
745 ELSE
746 IF(PIC_P_S(IJK,1).GT.ZERO) then
747 PS_FORCE_PIC(IJK,3) = 2.d0*(PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/(Dz(k)+Dz(k_of(ijkp)))
748 ELSE
749 PS_FORCE_PIC(IJK,3) = ZERO
750 ENDIF
751 ENDIF
752 ENDIF
753 ENDDO
754
755
756
757 = 1
758 DO K1 = KSTART3, KEND3
759 DO J1 = JSTART3, JEND3
760
761 IF(I1.NE.ISTART2) EXIT
762 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
763 IJK = FUNIJK(I1,J1,K1)
764 I = I_OF(IJK)
765 J = J_OF(IJK)
766 K = K_OF(IJK)
767
768 IPJK = IP_OF(IJK)
769 IF(PIC_P_S(IPJK,1).GT.ZERO) then
770 PS_FORCE_PIC(IJK,1) = 2.d0*(PIC_P_S(IPJK,1) - PIC_P_S(IJK,1))/(DX(I)+DX(I_of(ipjk)))
771 else
772 PS_FORCE_PIC(IJK,1) = zero
773 endif
774
775 ENDDO
776 ENDDO
777 J1 = 1
778 DO K1 = KSTART3, KEND3
779 DO I1 = ISTART3, IEND3
780
781 IF(J1.NE.JSTART2) EXIT
782 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
783 IJK = FUNIJK(I1,J1,K1)
784 I = I_OF(IJK)
785 J = J_OF(IJK)
786 K = K_OF(IJK)
787
788
789 IJPK = JP_OF(IJK)
790
791 IF(PIC_P_S(IJPK,1).GT.ZERO) then
792
793 PS_FORCE_PIC(IJK,2) = 2.d0*(PIC_P_S(IJPK,1) - PIC_P_S(IJK,1))/(DY(j)+Dy(j_of(ijpk)))
794 ELSE
795 PS_FORCE_PIC(IJK,2) = ZERO
796 ENDIF
797
798 END DO
799 END DO
800
801 IF(DO_K) then
802 K1 = 1
803 DO J1 = JSTART3, JEND3
804 DO I1 = ISTART3, IEND3
805 IF(K1.NE.KSTART2) EXIT
806 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
807 IJK = FUNIJK(I1,J1,K1)
808 I = I_OF(IJK)
809 J = J_OF(IJK)
810 K = K_OF(IJK)
811
812
813 IJKP = KP_OF(IJK)
814
815 IF(PIC_P_S(IJKP,1).GT.ZERO) then
816
817 PS_FORCE_PIC(IJK,3) = 2.d0*(PIC_P_S(IJKP,1) - PIC_P_S(IJK,1))/(Dz(k)+Dz(k_of(ijkp)))
818 ELSE
819 PS_FORCE_PIC(IJK, 3) = ZERO
820 ENDIF
821
822 END DO
823 END DO
824 ENDIF
825
826 DO IDIM = 1, merge(2,3,NO_K)
827 CALL SEND_RECV(PS_FORCE_PIC(:,IDIM),1)
828 ENDDO
829
830 CALL SET_INTERPOLATION_SCHEME(2)
831
832 KORDER = merge ( 1, 2, no_k)
833
834
835
836 = merge(0.5d0, 0.25D0, NO_K)
837
838 do ijk = ijkstart3,ijkend3
839
840 if(.not.fluid_at(ijk) .or. pinc(ijk).eq.0) cycle
841 i = i_of(ijk)
842 j = j_of(ijk)
843 k = k_of(ijk)
844 pcell(1) = i-1
845 pcell(2) = j-1
846 pcell(3) = merge(1, k-1, no_k)
847 call set_interpolation_stencil(pcell,iw,ie,js,jn,kb,&
848 ktp,interp_scheme,merge(2,3,no_k),ordernew = onew)
849
850
851 do k = 1, merge(1, onew, no_K)
852 do j = 1,onew
853 do i = 1,onew
854 ii = iw + i-1
855 jj = js + j-1
856 kk = kb + k-1
857 cur_ijk = funijk_map_c(ii,jj,kk)
858 ipjk = funijk_map_c(ii+1,jj,kk)
859 ijpk = funijk_map_c(ii,jj+1,kk)
860 ipjpk = funijk_map_c(ii+1,jj+1,kk)
861
862 vol_ijk = zero
863 vol_ipjk = zero
864 vol_ijpk = zero
865 vol_ipjpk = zero
866
867 vol_ijkp = zero
868 vol_ipjkp = zero
869 vol_ijpkp = zero
870 vol_ipjpkp = zero
871
872 if(fluid_at(cur_ijk)) vol_ijk = vol(cur_ijk)
873 if(fluid_at(ipjk)) vol_ipjk = vol(ipjk)
874 if(fluid_at(ijpk)) vol_ijpk = vol(ijpk)
875 if(fluid_at(ipjpk)) vol_ipjpk = vol(ipjpk)
876
877
878 if(DO_K) then
879 ijkp = funijk_map_c(ii,jj,kk+1)
880 ijpkp = funijk_map_c(ii,jj+1,kk+1)
881 ipjkp = funijk_map_c(ii+1,jj,kk+1)
882 ipjpkp = funijk_map_c(ii+1,jj+1,kk+1)
883
884
885 if(fluid_at(ijkp)) vol_ijkp = vol(ijkp)
886 if(fluid_at(ipjkp)) vol_ipjkp = vol(ipjkp)
887 if(fluid_at(ijpkp)) vol_ijpkp = vol(ijpkp)
888 if(fluid_at(ipjpkp)) vol_ipjpkp = vol(ipjpkp)
889
890 endif
891 gstencil(i,j,k,1) = xe(ii)
892 gstencil(i,j,k,2) = yn(jj)
893 gstencil(i,j,k,3) = merge(DZ(1), ZT(KK), NO_K)
894
895 VOL_TOT_SCAL = ZERO
896
897 VOL_TOT_SCAL = vol_ijk + vol_ipjk + vol_ijpk + vol_ipjpk + &
898 & vol_ijkp + vol_ipjkp + vol_ijpkp + vol_ipjpkp
899
900 VOL_TOT_VEC = ZERO
901
902 VOL_TOT_VEC(1) = VOL(CUR_IJK) + VOL(IJPK)
903 VOL_TOT_VEC(2) = VOL(CUR_IJK) + VOL(IPJK)
904
905 DO M = 1, DES_MMAX
906 VEL_SOL_STENCIL(i,j,k,1,M) = pic_u_s(cur_ijk,M)*vol(cur_ijk) + pic_u_s(ijpk,M)*vol(ijpk)
907
908 VEL_SOL_STENCIL(i,j,k,2,M) = pic_v_s(cur_ijk,M)*vol(cur_ijk) + pic_v_s(ipjk,M)*vol(ipjk)
909 ENDDO
910
911
912 (i,j,k) = ep_g(cur_ijk)*vol_ijk+ ep_g(ipjk)*vol_ipjk+ &
913 & ep_g(ijpk)*vol_ijpk + ep_g(ipjpk)*vol_ipjpk
914
915 psgradstencil(i,j,k,1) = PS_FORCE_PIC(cur_ijk,1)*VOL(CUR_IJK) + PS_FORCE_PIC(ijpk,1)*VOL(IJPK)
916
917 psgradstencil(i,j,k,2) = PS_FORCE_PIC(cur_ijk,2)*VOL(CUR_IJK) + PS_FORCE_PIC(ipjk,2)*VOL(IPJK)
918
919 vstencil(i,j,k,1) = u_g(cur_ijk)*vol(cur_ijk) + u_g(ijpk)*vol(ijpk)
920 vstencil(i,j,k,2) = v_g(cur_ijk)*vol(cur_ijk) + v_g(ipjk)*vol(ipjk)
921 if(DO_K) then
922 VOL_TOT_VEC(1) = VOL_TOT_VEC(1) + VOL(IJKP) + VOL(IJPKP)
923 VOL_TOT_VEC(2) = VOL_TOT_VEC(2) + VOL(IJKP) + VOL(IPJKP)
924 VOL_TOT_VEC(3) = VOL(CUR_IJK) + VOL(IPJK) + VOL(IJPK) + VOL(IPJPK)
925 sstencil(i,j,k) = sstencil(i,j,k) + ep_g(ijkp)*vol_ijkp + ep_g(ipjkp)*vol_ipjkp &
926 & + ep_g(ijpkp)*vol_ijpkp + ep_g(ipjpkp)*vol_ipjpkp
927
928 psgradstencil(i,j,k,1) = psgradstencil(i,j,k,1) &
929 + PS_FORCE_PIC(ijkp,1)*VOL(ijkp) + PS_FORCE_PIC(ijpkp,1)*vol(ijpkp)
930 psgradstencil(i,j,k,2) = psgradstencil(i,j,k,2) &
931 + PS_FORCE_PIC(ijkp,2)*vol(ijkp) + PS_FORCE_PIC(ipjkp,2)*vol(ipjkp)
932 psgradstencil(i,j,k,3) = PS_FORCE_PIC(cur_ijk,3)*vol(cur_ijk)+&
933 & PS_FORCE_PIC(ijpk,3)*vol(ijpk)+PS_FORCE_PIC(ipjk,3)*vol(ipjk)+&
934 & PS_FORCE_PIC(ipjpk,3)*vol(ipjpk)
935
936 vstencil(i,j,k,1) = vstencil(i,j,k,1) + u_g(ijkp)*vol(ijkp) + u_g(ijpkp)*vol(ijpkp)
937
938 vstencil(i,j,k,2) = vstencil(i,j,k,2) + v_g(ijkp)*vol(ijkp) + v_g(ipjkp)*vol(ipjkp)
939 vstencil(i,j,k,3) = w_g(cur_ijk)*vol(cur_ijk)+&
940 & w_g(ijpk)*vol(ijpk)+w_g(ipjk)*vol(ipjk)+w_g(ipjpk)*vol(ipjpk)
941
942 DO M = 1, DES_MMAX
943 VEL_SOL_STENCIL(i,j,k,1, M) = VEL_SOL_STENCIL(i,j,k,1,M) &
944 & + PIC_u_s(ijkp,M)*vol(ijkp) + PIC_u_s(ijpkp,M)*vol(ijpkp)
945 VEL_SOL_STENCIL(i,j,k,2, M) = VEL_SOL_STENCIL(i,j,k,2,M) &
946 & + PIC_v_s(ijkp,M)*vol(ijkp) + PIC_v_s(ipjkp,M)*vol(ipjkp)
947 VEL_SOL_STENCIL(i,j,k,3, M) = PIC_w_s(cur_ijk,M)*vol(cur_ijk) +&
948 PIC_w_s(ijpk,M)*vol(ijpk)+PIC_w_s(ipjk,M)*vol(ipjk)+PIC_w_s(ipjpk,M)*vol(ipjpk)
949 ENDDO
950 else
951 psgradstencil(i,j,k,3) = 0.d0
952 VEL_SOL_STENCIL(i,j,k,3, 1:DES_MMAX) = 0.d0
953 vstencil(i,j,k,3) = 0.d0
954
955 endif
956 DO IDIM = 1, merge(2,3,NO_K)
957 IF(VOL_TOT_VEC(IDIM).GT.ZERO) THEN
958 psgradstencil(i,j,k,idim) = psgradstencil(i,j,k,idim)/VOL_TOT_VEC(idim)
959
960 VEL_SOL_STENCIL(i,j,k,idim, 1:DES_MMAX) = VEL_SOL_STENCIL(i,j,k,idim, 1:DES_MMAX)/VOL_TOT_VEC(idim)
961
962 vstencil(i,j,k,idim) = vstencil(i,j,k,idim)/VOL_TOT_VEC(idim)
963
964
965
966 ENDIF
967
968 ENDDO
969
970
971
972
973
974
975
976 if(VOL_TOT_SCAL.gt.zero) sstencil(i,j,k) = sstencil(i,j,k)/VOL_TOT_SCAL
977
978
979 enddo
980 enddo
981 enddo
982
983
984 do nindx = 1,pinc(ijk)
985 np = pic(ijk)%p(nindx)
986 m = pijk(np,5)
987
988 if(NO_K) then
989 call interpolator(gstencil(1:onew,1:onew,1,1:dimn), &
990 psgradstencil(1:onew,1:onew,1,1:dimn), &
991 des_pos_new(1:dimn,np),PS_GRAD(np,1:dimn), &
992 onew,interp_scheme,weightp)
993 else
994 call interpolator(gstencil(1:onew,1:onew,1:onew,1:dimn), &
995 psgradstencil(1:onew,1:onew,1:onew,1:dimn), &
996 des_pos_new(1:dimn,np),PS_GRAD(np,1:dimn), &
997 onew,interp_scheme,weightp)
998 endif
999
1000 do idim = 1, merge(2,3,NO_K)
1001 AVGSOLVEL_P(IDIM,NP) = ARRAY_DOT_PRODUCT( &
1002 VEL_SOL_STENCIL(:,:,:,IDIM,M),WEIGHTP(:,:,:))
1003 ENDDO
1004
1005 EPG_P(NP) = ARRAY_DOT_PRODUCT(SSTENCIL(:,:,:),WEIGHTP(:,:,:))
1006
1007 END DO
1008 END DO
1009
1010
1011 END SUBROUTINE MPPIC_COMPUTE_PS_GRAD
1012
1013
1014
1015
1016
1017
1018
1019
1020 SUBROUTINE MPPIC_BC_U_S
1021
1022
1023
1024
1025 USE parallel
1026 USE constant
1027 USE geometry
1028 USE indices
1029 USE compar
1030 USE mfix_pic
1031 USE functions
1032 IMPLICIT NONE
1033
1034
1035
1036
1037 INTEGER :: IER
1038
1039 INTEGER :: L
1040
1041 INTEGER :: I, J, K, IM, I1, I2, J1, J2, K1, K2, IJK,&
1042 JM, KM, IJKW, IMJK, IPJK, IP, IJK_WALL
1043
1044
1045
1046 IF (DO_K) THEN
1047 K1 = 1
1048 DO J1 = jmin3,jmax3
1049 DO I1 = imin3, imax3
1050 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1051 IJK_WALL = FUNIJK(I1,J1,K1)
1052 IJK = FUNIJK(I1,J1,K1+1)
1053 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1054 END DO
1055 END DO
1056
1057 K1 = KMAX2
1058 DO J1 = jmin3,jmax3
1059 DO I1 = imin3, imax3
1060 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1061 IJK_WALL = FUNIJK(I1,J1,K1)
1062 IJK = FUNIJK(I1,J1,K1-1)
1063 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1064 END DO
1065 END DO
1066 ENDIF
1067
1068 J1 = 1
1069 DO K1 = kmin3, kmax3
1070 DO I1 = imin3, imax3
1071 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1072 IJK_WALL = FUNIJK(I1,J1,K1)
1073 IJK = FUNIJK(I1,J1+1,K1)
1074 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1075 END DO
1076 END DO
1077 J1 = JMAX2
1078 DO K1 = kmin3, kmax3
1079 DO I1 = imin3, imax3
1080 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1081 IJK_WALL = FUNIJK(I1,J1,K1)
1082 IJK = FUNIJK(I1,J1-1,K1)
1083 PIC_U_S(IJK_WALL, :) = PIC_U_S(IJK,:)
1084 END DO
1085 END DO
1086
1087 RETURN
1088 END SUBROUTINE MPPIC_BC_U_S
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098 SUBROUTINE MPPIC_BC_V_S
1099
1100
1101
1102
1103
1104 USE parallel
1105 USE constant
1106 USE geometry
1107 USE indices
1108 USE compar
1109 USE mfix_pic
1110 USE functions
1111 IMPLICIT NONE
1112
1113
1114
1115
1116 INTEGER IER
1117
1118 INTEGER L
1119
1120 INTEGER I, J, K, JM, I1, I2, J1, J2, K1, K2, IJK,&
1121 IM, KM, IJKS, IJMK, IJPK, IJK_WALL
1122
1123
1124
1125 IF (DO_K) THEN
1126 K1 = 1
1127 DO J1 = jmin3,jmax3
1128 DO I1 = imin3, imax3
1129 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1130 IJK_WALL = FUNIJK(I1,J1,K1)
1131 IJK = FUNIJK(I1,J1,K1+1)
1132 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1133 END DO
1134 END DO
1135 K1 = KMAX2
1136 DO J1 = jmin3,jmax3
1137 DO I1 = imin3, imax3
1138 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1139 IJK_WALL = FUNIJK(I1,J1,K1)
1140 IJK = FUNIJK(I1,J1,K1-1)
1141 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1142 END DO
1143 END DO
1144 ENDIF
1145
1146 I1 = 1
1147 DO K1 = kmin3, kmax3
1148 DO J1 = jmin3, jmax3
1149 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1150 IJK_WALL = FUNIJK(I1,J1,K1)
1151 IJK = FUNIJK(I1+1,J1,K1)
1152 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1153
1154 END DO
1155 END DO
1156 I1 = IMAX2
1157 DO K1 = kmin3, kmax3
1158 DO J1 = jmin3, jmax3
1159 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1160 IJK_WALL = FUNIJK(I1,J1,K1)
1161 IJK = FUNIJK(I1-1,J1,K1)
1162 PIC_V_S(IJK_WALL, :) = PIC_V_S(IJK,:)
1163
1164 END DO
1165 END DO
1166 END SUBROUTINE MPPIC_BC_V_S
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176 SUBROUTINE MPPIC_BC_W_S
1177
1178
1179
1180
1181 USE parallel
1182 USE constant
1183 USE geometry
1184 USE indices
1185 USE compar
1186 USE mfix_pic
1187 USE functions
1188 IMPLICIT NONE
1189
1190
1191
1192
1193 INTEGER :: IER
1194
1195 INTEGER :: L
1196
1197 INTEGER :: I, J, K, KM, I1, I2, J1, J2, K1, K2, IJK,&
1198 IM, JM, IJKB, IJKM, IJKP, IJK_WALL
1199
1200
1201
1202 = 1
1203 DO K1 = kmin3,kmax3
1204 DO I1 = imin3,imax3
1205 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1206 IJK = FUNIJK(I1,J1+1,K1)
1207 IJK_WALL = FUNIJK(I1,J1,K1)
1208 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1209 END DO
1210 END DO
1211 J1 = JMAX2
1212 DO K1 = kmin3, kmax3
1213 DO I1 = imin3, imax3
1214 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1215 IJK = FUNIJK(I1,J1-1,K1)
1216 IJK_WALL = FUNIJK(I1,J1,K1)
1217 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1218 END DO
1219 END DO
1220 I1 = 1
1221 DO K1 = kmin3, kmax3
1222 DO J1 = jmin3, jmax3
1223 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1224 IJK = FUNIJK(I1+1,J1,K1)
1225 IJK_WALL = FUNIJK(I1,J1,K1)
1226 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1227 END DO
1228 END DO
1229 I1 = IMAX2
1230 DO K1 = kmin3,kmax3
1231 DO J1 = jmin3,jmax3
1232 IF (.NOT.IS_ON_myPE_owns(I1,J1,K1)) CYCLE
1233 IJK = FUNIJK(I1-1,J1,K1)
1234 IJK_WALL = FUNIJK(I1,J1,K1)
1235 PIC_W_S(IJK_WALL,:) = PIC_W_S(IJK,:)
1236 END DO
1237 END DO
1238
1239
1240 END SUBROUTINE MPPIC_BC_W_S
1241
1242
1243
1244
1245
1246
1247
1248 SUBROUTINE WRITE_NODEDATA(bufin, funit)
1249 USE param
1250 USE param1
1251 USE parallel
1252 USE constant
1253 USE run
1254 USE geometry
1255 USE indices
1256 USE compar
1257 USE discretelement
1258 use desmpi
1259 USE functions
1260 USE fun_avg
1261
1262 IMPLICIT NONE
1263 integer, intent(in) :: funit
1264 double precision, dimension(:), intent(in) :: bufin
1265
1266 integer :: ijk, i, j,k
1267
1268 write(funit,*)'VARIABLES= ',' "I" ',' "J" ',' "K" ',' "DES_ROPS_NODE" '
1269
1270 write(funit, *)'ZONE F=POINT, I=', (IEND1-ISTART2)+1, ', J=', JEND1-JSTART2+1, ', K=', KEND1-KSTART2 + 1
1271
1272 DO K = KSTART2, KEND1
1273 DO J = JSTART2, JEND1
1274 DO I = ISTART2, IEND1
1275 IJK = funijk(I, J, K)
1276
1277
1278 WRITE(funit, '(3(2x, i10), 3x, g17.8)') I, J, K , DES_ROPS_NODE(IJK,1)
1279 ENDDO
1280 ENDDO
1281 ENDDO
1282 CLOSE(funit, status = 'keep')
1283 end SUBROUTINE WRITE_NODEDATA
1284
1285
1286
1287
1288
1289
1290
1291 SUBROUTINE WRITE_MPPIC_VEL_S
1292 USE param
1293 USE param1
1294 USE parallel
1295 USE constant
1296 USE run
1297 USE geometry
1298 USE indices
1299 USE compar
1300 USE fldvar, only : ep_g
1301 USE discretelement
1302 USE mfix_pic
1303 USE functions
1304 implicit none
1305 integer :: i, j, k, ijk, fluid_ind, LL, PC, IDIM
1306 double precision :: zcor
1307 character(LEN=100) :: filename
1308 logical finish
1309
1310 WRITE(filename,'(A,"_",I5.5,".dat")') TRIM(RUN_NAME)//'_U_S_',myPE
1311 OPEN(1000, file = TRIM(filename), form ='formatted', status='unknown')
1312 IF(DIMN.eq.2) then
1313 write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
1314 ' "EP_s " ', ' "U_S" ', ' "V_S" ',' "DES_U_s" ', ' "DES_V_s" '
1315 write(1000,*)'ZONE F=POINT, I=', (IEND3-ISTART3)+1, ', J=', JEND3-JSTART3+1, ', K=', KEND3-KSTART3 + 1
1316 else
1317 write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
1318 ' "EP_s " ', ' "U_S" ', ' "V_S" ', ' "W_S" ',' "DES_U_s" ', ' "DES_V_s" ', ' "DES_W_s" '
1319
1320 write(1000,*)'ZONE F=POINT, I=', (IEND3-ISTART3)+1, ', J=', JEND3-JSTART3+1, ', K=', KEND3-KSTART3 + 1
1321 ENDIF
1322 DO K=KSTART3, KEND3
1323 DO J=JSTART3, JEND3
1324 DO I=ISTART3, IEND3
1325 IJK = FUNIJK(I,J,K)
1326 IF(FLUID_AT(IJK)) THEN
1327 FLUID_IND = 1
1328 ELSE
1329 FLUID_IND = 0
1330 END IF
1331 IF(DIMN.EQ.2) THEN
1332 ZCOR = ZT(K)
1333 WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))') XE(I-1)+DX(I), YN(J-1)+DY(J),ZCOR, &
1334 1.D0 - EP_G(IJK), PIC_U_S(IJK,1), PIC_V_S(IJK,1), DES_U_S(IJK,1), &
1335 DES_V_S(IJK,1)
1336 ELSE
1337 ZCOR = ZT(K-1) + DZ(K)
1338 WRITE(1000,'(3(2X,G17.8),4( 2X, G17.8))') XE(I-1)+DX(I), YN(J-1)+DY(J),ZCOR, &
1339 1.D0 - EP_G(IJK), PIC_U_S(IJK,1), PIC_V_S(IJK,1), PIC_W_S(IJK,1), DES_U_S(IJK,1), &
1340 DES_V_S(IJK,1), DES_W_S(IJK,1)
1341 ENDIF
1342 ENDDO
1343 ENDDO
1344 ENDDO
1345 CLOSE(1000, STATUS='KEEP')
1346
1347 return
1348
1349 WRITE(FILENAME,'(A,"_",I5.5,".DAT")') TRIM(RUN_NAME)//'_PS_FORCE_',myPE
1350 OPEN(1000, file = TRIM(filename), form ='formatted', status='unknown')
1351
1352 IF(DIMN.eq.3) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ',' "Z" ', &
1353 ' "DELPX" ', '"DELPY"', '"DELPZ" ',' "US_part" ', '"VS_part"' , '"WS_part"', '"EP_s_part"'
1354 IF(DIMN.eq.2) write(1000,*)'VARIABLES= ',' "X" ',' "Y" ', &
1355 ' "DELPX" ', '"DELPY"', ' "US_part" ', '"VS_part"' , '"EP_S_part"'
1356
1357 PC = 1
1358 DO LL = 1, MAX_PIP
1359
1360 IF(PC .GT. PIP) EXIT
1361 IF(.NOT.PEA(LL,1)) CYCLE
1362 pc = pc+1
1363 IF(PEA(LL,4)) CYCLE
1364
1365 WRITE(1000,'(10( 2x, g17.8))') (DES_POS_NEW(IDIM, LL), IDIM = 1, DIMN), &
1366 (PS_GRAD(LL, IDIM) , IDIM = 1, DIMN), (AVGSOLVEL_P (IDIM, LL) , IDIM = 1, DIMN), 1-EPg_P(LL)
1367 ENDDO
1368 close(1000, status='keep')
1369
1370
1371
1372
1373 END SUBROUTINE WRITE_MPPIC_VEL_S
1374
1375
1376
1377
1378
1379