File: RELATIVE:/../../../mfix.git/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: Coupled gas-solids DEM simulation
19           use discretelement, only: DES_CONTINUUM_HYBRID
20     ! Flag: Discrete and continuum solids co-exist
21           use discretelement, only: DES_CONTINUUM_COUPLED
22     ! Number of discrete solids phases
23           use discretelement, only: DES_MMAX
24     ! Global ID of particles
25           use discretelement, only: iGLOBAL_ID
26     ! Particle positions
27           use discretelement, only: DES_POS_NEW
28     ! Number of continuum solids phases
29           use physprop, only: SMAX
30     ! Discrete particle material and bulk densities
31           use discretelement, only: DES_RO_S, DES_ROP_s
32     ! Number of particles in indexed fluid cell
33           use discretelement, only: PINC
34     ! List of particles in each cell.
35           use discretelement, only: PIC
36     ! Gas phae volume fraction, density, and build density
37           use fldvar, only: EP_G, RO_G, ROP_G
38     ! Bulk density of continuum solids phases
39           use fldvar, only: EP_S
40     ! Volume of scalar grid cell.
41           use geometry, only: VOL
42     ! Flag: Status of indexed cell
43           use cutcell, only: CUT_CELL_AT
44     ! Flag: Indexed cell contains fluid
45           USE functions, only: FLUID_AT
46     ! Fluid grid loop bounds.
47           use compar, only: IJKStart3, IJKEnd3
48     ! Flag: Fluid exists at indexed cell
49           use functions, only: FLUID_AT
50     ! The I, J, and K values that comprise an IJK
51           use indices, only: I_OF, J_OF, K_OF
52     ! Rank ID of current process
53           use compar, only: myPE
54     ! Global communication function to sum to all ranks.
55           use mpi_utility, only: GLOBAL_ALL_SUM
56     
57     ! Global Parameters:
58     !---------------------------------------------------------------------//
59           USE param1, only: ZERO, ONE
60     
61           use error_manager
62     
63           IMPLICIT NONE
64     
65     ! Local Variables:
66     !---------------------------------------------------------------------//
67     ! Loop indices
68           INTEGER :: IJK, M, LC
69     ! Total solids volume fraction
70           DOUBLE PRECISION SUM_EPS
71     ! Integer Error Flag
72           INTEGER :: IER
73     !......................................................................!
74     
75     ! Initialize error flag.
76           IER = 0
77     
78     ! Calculate gas volume fraction from solids volume fraction:
79     !---------------------------------------------------------------------//
80     !$omp parallel do if(ijkend3 .ge. 2000) default(none) reduction(+:IER) &
81     !$omp shared(IJKSTART3, IJKEND3, DES_CONTINUUM_COUPLED, DES_MMAX, SMAX,&
82     !$omp    EP_G, RO_G, ROP_G, DES_ROP_S, DES_RO_S, DES_CONTINUUM_HYBRID) &
83     !$omp private(IJK, SUM_EPs, M)
84           DO IJK = IJKSTART3, IJKEND3
85     ! Skip wall cells.
86              IF(.NOT.FLUID_AT(IJK)) CYCLE
87     ! Initialize EP_g and the accumulator.
88              EP_G(IJK) = ONE
89              SUM_EPS = ZERO
90     ! Sum the DES solids volume fraction.
91              DO M = 1, DES_MMAX
92                 SUM_EPS = SUM_EPS + DES_ROP_S(IJK,M)/DES_RO_S(M)
93              ENDDO
94     ! Add in TFM solids contributions.
95              IF(DES_CONTINUUM_HYBRID) THEN
96                 DO M = 1,SMAX
97                    SUM_EPS = SUM_EPS + EP_S(IJK,M)
98                 ENDDO
99              ENDIF
100     ! Calculate the gas phase volume fraction and bulk density.
101              EP_G(IJK) = ONE - SUM_EPS
102              ROP_G(IJK) = RO_G(IJK) * EP_G(IJK)
103     ! Flag an error if gas volume fraction is unphysical.
104              IF(DES_CONTINUUM_COUPLED) THEN
105                 IF(EP_G(IJK) <= ZERO .OR. EP_G(IJK) > ONE) IER = IER + 1
106              ENDIF
107           ENDDO
108     !omp end parallel do
109     
110     
111           CALL GLOBAL_ALL_SUM(IER)
112           IF(IER == 0) RETURN
113     
114     
115     ! Report any errors. Volume fraction errors are fatal.
116     !---------------------------------------------------------------------//
117           CALL INIT_ERR_MSG("CALC_EPG_DES")
118           CALL OPEN_PE_LOG(IER)
119     
120           WRITE(ERR_MSG, 1100)
121           CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
122     
123      1100 FORMAT('Error 1100: Unphysical gas phase volume fraction ',      &
124              'calculated. A .vtp',/'file will be written and the code ',   &
125              'will exit. Fluid cell details:')
126     
127              DO IJK=IJKSTART3, IJKEND3
128                 IF(.NOT.FLUID_AT(IJK)) CYCLE
129                 IF(EP_G(IJK) > ZERO .AND. EP_G(IJK) <= ONE) CYCLE
130     
131                 WRITE(ERR_MSG,1101) trim(iVal(IJK)), trim(iVal(I_OF(IJK))),&
132                    trim(iVal(J_OF(IJK))), trim(iVal(K_OF(IJK))),EP_G(IJK), &
133                    CUT_CELL_AT(IJK), trim(iVal(PINC(IJK))), VOL(IJK)
134                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
135     
136                 WRITE(ERR_MSG,1102)
137                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
138                 DO LC=1,PINC(IJK)
139                    M=PIC(IJK)%P(LC)
140                    WRITE(ERR_MSG,1103) iGlobal_ID(M), trim(iVal(           &
141                       DES_POS_NEW(1,M))), trim(iVal(DES_POS_NEW(2,M))),    &
142                       trim(iVal(DES_POS_NEW(3,M)))
143                    CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
144                 ENDDO
145              ENDDO
146     
147      1101 FORMAT(/3x,'Fluid Cell IJK: ',A,6x,'I/J/K: (',A,',',A,',',A,')',/&
148              T6,'EP_G = ',g11.4,T30,'CUT_CELL_AT = ',L1,/T6,'PINC: ',A,T30,&
149              'VOL = ',g11.4)
150     
151      1102 FORMAT(/T6,'Global ID',T30,'Position')
152     
153      1103 FORMAT(T6,I9,3x,'(',A,', ',A,', ',A,')')
154     
155           WRITE(ERR_MSG, 1104)
156           CALL FLUSH_ERR_MSG(HEADER=.FALSE.)
157      1104 FORMAT('This is a fatal error. A particle output file (vtp) ',   &
158              'will be written',/'to aid debugging.')
159     
160           CALL WRITE_DES_DATA
161           CALL MFIX_EXIT(myPE)
162     
163           END SUBROUTINE CALC_EPG_DES
164