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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: CALC_EPG_DES                                            !
3     !  Author: R.Garg                                     Date: ??-???-??  !
4     !                                                                      !
5     !  Purpose: Calculate the gas phase volume fraction (and in turn the   !
6     !  gas phase bulk density) from the sum of the solids volume fractions.!
7     !                                                                      !
8     !  NOTE: This routine uses a global communication to notify all ranks  !
9     !  of potential errors. Therefore all ranks can call MFIX_EXIT and     !
10     !  prevent dead-lock. Communications may be reduced by passing the     !
11     !  flag back to the caller and combining with other error checks.      !
12     !                                                                      !
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
14           SUBROUTINE CALC_EPG_DES
15     
16     ! Global Variables:
17     !---------------------------------------------------------------------//
18     ! Flag: Discrete and continuum solids co-exist
19           use discretelement, only: DES_CONTINUUM_COUPLED
20     ! Number of discrete solids phases
21           use discretelement, only: DES_MMAX
22     ! Global ID of particles
23           use discretelement, only: iGLOBAL_ID
24     ! Particle positions
25           use discretelement, only: DES_POS_NEW
26     ! Number of continuum solids phases
27           use physprop, only: MMAX
28     ! Number of particles in indexed fluid cell
29           use discretelement, only: PINC
30     ! List of particles in each cell.
31           use derived_types, only: PIC
32     ! Gas phae volume fraction, density, and build density
33           use fldvar, only: EP_G, RO_G, ROP_G
34     ! Bulk density of continuum solids phases
35           use fldvar, only: EP_S
36     ! Volume of scalar grid cell.
37           use geometry, only: VOL
38     ! Flag: Status of indexed cell
39           use cutcell, only: CUT_CELL_AT
40     ! Flag: Indexed cell contains fluid
41           USE functions, only: FLUID_AT
42     ! Fluid grid loop bounds.
43           use compar, only: IJKStart3, IJKEnd3
44     ! Flag: Fluid exists at indexed cell
45           use functions, only: FLUID_AT
46     ! The I, J, and K values that comprise an IJK
47           use indices, only: I_OF, J_OF, K_OF
48     ! Rank ID of current process
49           use compar, only: myPE
50     ! Global communication function to sum to all ranks.
51           use mpi_utility, only: GLOBAL_ALL_SUM
52     ! Flag for PIC simulation
53           use mfix_pic, only: MPPIC
54     
55     ! Global Parameters:
56     !---------------------------------------------------------------------//
57           USE param1, only: ZERO, ONE
58     
59           use error_manager
60     
61           IMPLICIT NONE
62     
63     ! Local Variables:
64     !---------------------------------------------------------------------//
65     ! Loop indices
66           INTEGER :: IJK, M, LC
67     ! Total solids volume fraction
68           DOUBLE PRECISION :: SUM_EPS
69     ! Packed
70           DOUBLE PRECISION :: PACKED_EPS
71     ! Integer Error Flag
72           INTEGER :: IER
73     !......................................................................!
74     
75     ! Initialize error flag.
76           IER = 0
77     
78     ! Set a max solids volume fraction for MPPIC. The value is arbitrarily
79     ! larger than EP_STAR. However, the value should be large enough so
80     ! that it is rarely used. This was added as a crude work around for
81     ! poor initial conditions that start cells overpacked.
82           PACKED_EPS = merge(0.9d0, ONE, MPPIC)
83     
84     ! Calculate gas volume fraction from solids volume fraction:
85     !---------------------------------------------------------------------//
86     !$omp parallel do if(ijkend3 .ge. 2000) default(none) reduction(+:IER) &
87     !$omp shared(IJKSTART3, IJKEND3, DES_CONTINUUM_COUPLED, DES_MMAX, MMAX,&
88     !$omp        EP_G, RO_G, ROP_G, MPPIC, PACKED_EPS) &
89     !$omp private(IJK, SUM_EPs, M)
90           DO IJK = IJKSTART3, IJKEND3
91     ! Skip wall cells.
92              IF(.NOT.FLUID_AT(IJK)) CYCLE
93     ! Initialize EP_g and the accumulator.
94              EP_G(IJK) = ONE
95              SUM_EPS = ZERO
96     ! Sum the DES solids volume fraction.  If hybrid TFM solids contributions
97     ! will be included
98              DO M = 1, DES_MMAX+MMAX
99                 SUM_EPS = SUM_EPS + EP_S(IJK,M)
100              ENDDO
101     
102     ! Calculate the gas phase volume fraction.
103              EP_G(IJK) = ONE - min(PACKED_EPS, SUM_EPS)
104     
105     ! Calculate the gas phase bulk density.
106              ROP_G(IJK) = RO_G(IJK) * EP_G(IJK)
107     ! Flag an error if gas volume fraction is unphysical.
108              IF(DES_CONTINUUM_COUPLED) THEN
109                 IF(EP_G(IJK) <= ZERO .OR. EP_G(IJK) > ONE) IER = IER + 1
110              ENDIF
111           ENDDO
112     !omp end parallel do
113     
114     
115           CALL GLOBAL_ALL_SUM(IER)
116           IF(IER == 0) RETURN
117     
118     
119     ! Report any errors. Volume fraction errors are fatal.
120     !---------------------------------------------------------------------//
121           CALL INIT_ERR_MSG("CALC_EPG_DES")
122           CALL OPEN_PE_LOG(IER)
123     
124           WRITE(ERR_MSG, 1100)
125           CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
126     
127      1100 FORMAT('Error 1100: Unphysical gas phase volume fraction ',      &
128              'calculated. A .vtp',/'file will be written and the code ',   &
129              'will exit. Fluid cell details:')
130     
131              DO IJK=IJKSTART3, IJKEND3
132                 IF(.NOT.FLUID_AT(IJK)) CYCLE
133                 IF(EP_G(IJK) > ZERO .AND. EP_G(IJK) <= ONE) CYCLE
134     
135                 WRITE(ERR_MSG,1101) trim(iVal(IJK)), trim(iVal(I_OF(IJK))),&
136                    trim(iVal(J_OF(IJK))), trim(iVal(K_OF(IJK))),EP_G(IJK), &
137                    CUT_CELL_AT(IJK), trim(iVal(PINC(IJK))), VOL(IJK)
138                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
139     
140                 WRITE(ERR_MSG,1102)
141                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
142                 DO LC=1,PINC(IJK)
143                    M=PIC(IJK)%P(LC)
144                    WRITE(ERR_MSG,1103) iGlobal_ID(M), trim(iVal(           &
145                       DES_POS_NEW(M,1))), trim(iVal(DES_POS_NEW(M,2))),    &
146                       trim(iVal(DES_POS_NEW(M,3)))
147                    CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
148                 ENDDO
149              ENDDO
150     
151      1101 FORMAT(/3x,'Fluid Cell IJK: ',A,6x,'I/J/K: (',A,',',A,',',A,')',/&
152              T6,'EP_G = ',g11.4,T30,'CUT_CELL_AT = ',L1,/T6,'PINC: ',A,T30,&
153              'VOL = ',g11.4)
154     
155      1102 FORMAT(/T6,'Global ID',T30,'Position')
156     
157      1103 FORMAT(T6,I9,3x,'(',A,', ',A,', ',A,')')
158     
159           WRITE(ERR_MSG, 1104)
160           CALL FLUSH_ERR_MSG(HEADER=.FALSE.)
161      1104 FORMAT('This is a fatal error. A particle output file (vtp) ',   &
162              'will be written',/'to aid debugging.')
163     
164           CALL WRITE_DES_DATA
165           CALL MFIX_EXIT(myPE)
166     
167           END SUBROUTINE CALC_EPG_DES
168