calc_epg_des.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_EPG_DES                                            !
!  Author: R.Garg                                     Date: ??-???-??  !
!                                                                      !
!  Purpose: Calculate the gas phase volume fraction (and in turn the   !
!  gas phase bulk density) from the sum of the solids volume fractions.!
!                                                                      !
!  NOTE: This routine uses a global communication to notify all ranks  !
!  of potential errors. Therefore all ranks can call MFIX_EXIT and     !
!  prevent dead-lock. Communications may be reduced by passing the     !
!  flag back to the caller and combining with other error checks.      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_epg_des

! Global Variables:
!---------------------------------------------------------------------//
! Flag: Discrete and continuum solids co-exist
      use discretelement, only: des_continuum_coupled
! Number of discrete solids phases
      use discretelement, only: des_mmax
! Global ID of particles
      use discretelement, only: iglobal_id
! Particle positions
      use discretelement, only: des_pos_new
! Number of continuum solids phases
      use physprop, only: mmax
! Number of particles in indexed fluid cell
      use discretelement, only: pinc
! List of particles in each cell.
      use derived_types, only: pic
! Gas phae volume fraction, density, and build density
      use fldvar, only: ep_g, ro_g, rop_g
! Bulk density of continuum solids phases
      use fldvar, only: ep_s
! Volume of scalar grid cell.
      use geometry, only: vol
! Flag: Status of indexed cell
      use cutcell, only: cut_cell_at
! Flag: Indexed cell contains fluid
      USE functions, only: fluid_at
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! Flag: Fluid exists at indexed cell
      use functions, only: fluid_at
! The I, J, and K values that comprise an IJK
      use indices, only: i_of, j_of, k_of
! Rank ID of current process
      use compar, only: mype
! Global communication function to sum to all ranks.
      use mpi_utility, only: global_all_sum
! Flag for PIC simulation
      use mfix_pic, only: mppic

! Global Parameters:
!---------------------------------------------------------------------//
      USE param1, only: zero, one

      use error_manager

      IMPLICIT NONE

! Local Variables:
!---------------------------------------------------------------------//
! Loop indices
      INTEGER :: IJK, M, LC
! Total solids volume fraction
      DOUBLE PRECISION :: SUM_EPS
! Packed
      DOUBLE PRECISION :: PACKED_EPS
! Integer Error Flag
      INTEGER :: IER
!......................................................................!

! Initialize error flag.
      ier = 0

! Set a max solids volume fraction for MPPIC. The value is arbitrarily
! larger than EP_STAR. However, the value should be large enough so
! that it is rarely used. This was added as a crude work around for
! poor initial conditions that start cells overpacked.
      packed_eps = merge(0.9d0, one, mppic)

! Calculate gas volume fraction from solids volume fraction:
!---------------------------------------------------------------------//
!$omp parallel do if(ijkend3 .ge. 2000) default(none) reduction(+:IER) &
!$omp shared(IJKSTART3, IJKEND3, DES_CONTINUUM_COUPLED, DES_MMAX, MMAX,&
!$omp        EP_G, RO_G, ROP_G, MPPIC, PACKED_EPS) &
!$omp private(IJK, SUM_EPs, M)
      DO ijk = ijkstart3, ijkend3
! Skip wall cells.
         IF(.NOT.fluid_at(ijk)) cycle
! Initialize EP_g and the accumulator.
         ep_g(ijk) = one
         sum_eps = zero
! Sum the DES solids volume fraction.  If hybrid TFM solids contributions
! will be included
         DO m = 1, des_mmax+mmax
            sum_eps = sum_eps + ep_s(ijk,m)
         ENDDO

! Calculate the gas phase volume fraction.
         ep_g(ijk) = one - min(packed_eps, sum_eps)

! Calculate the gas phase bulk density.
         rop_g(ijk) = ro_g(ijk) * ep_g(ijk)
! Flag an error if gas volume fraction is unphysical.
         IF(des_continuum_coupled) THEN
            IF(ep_g(ijk) <= zero .OR. ep_g(ijk) > one) ier = ier + 1
         ENDIF
      ENDDO
!omp end parallel do


      CALL global_all_sum(ier)
      IF(ier == 0) RETURN


! Report any errors. Volume fraction errors are fatal.
!---------------------------------------------------------------------//
      CALL init_err_msg("CALC_EPG_DES")
      CALL open_pe_log(ier)

      WRITE(err_msg, 1100)
      CALL flush_err_msg(footer=.false.)

 1100 FORMAT('Error 1100: Unphysical gas phase volume fraction ',      &
         'calculated. A .vtp',/'file will be written and the code ',   &
         'will exit. Fluid cell details:')

         DO ijk=ijkstart3, ijkend3
            IF(.NOT.fluid_at(ijk)) cycle
            IF(ep_g(ijk) > zero .AND. ep_g(ijk) <= one) cycle

            WRITE(err_msg,1101) trim(ival(ijk)), trim(ival(i_of(ijk))),&
               trim(ival(j_of(ijk))), trim(ival(k_of(ijk))),ep_g(ijk), &
               cut_cell_at(ijk), trim(ival(pinc(ijk))), vol(ijk)
            CALL flush_err_msg(header=.false., footer=.false.)

            WRITE(err_msg,1102)
            CALL flush_err_msg(header=.false., footer=.false.)
            DO lc=1,pinc(ijk)
               m=pic(ijk)%P(lc)
               WRITE(err_msg,1103) iglobal_id(m), trim(ival(           &
                  des_pos_new(m,1))), trim(ival(des_pos_new(m,2))),    &
                  trim(ival(des_pos_new(m,3)))
               CALL flush_err_msg(header=.false., footer=.false.)
            ENDDO
         ENDDO

 1101 FORMAT(/3x,'Fluid Cell IJK: ',a,6x,'I/J/K: (',a,',',a,',',a,')',/&
         t6,'EP_G = ',g11.4,t30,'CUT_CELL_AT = ',l1,/t6,'PINC: ',a,t30,&
         'VOL = ',g11.4)

 1102 FORMAT(/t6,'Global ID',t30,'Position')

 1103 FORMAT(t6,i9,3x,'(',a,', ',a,', ',a,')')

      WRITE(err_msg, 1104)
      CALL flush_err_msg(header=.false.)
 1104 FORMAT('This is a fatal error. A particle output file (vtp) ',   &
         'will be written',/'to aid debugging.')

      CALL write_des_data
      CALL mfix_exit(mype)

      END SUBROUTINE calc_epg_des