set_bc_dem_mi.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SET_BC_DEM_MI                                           !
!  Author: J.Musser                                   Date: 23-Nov-09  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE set_bc_dem_mi

! Modules
!---------------------------------------------------------------------//
      USE bc
      USE compar
      USE constant
      USE des_allocate
      USE des_bc
      USE discretelement
      USE error_manager
      USE funits
      USE geometry
      USE indices
      USE param
      USE param1
      USE physprop
      USE run
      USE toleranc
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
      INTEGER :: BCV
      INTEGER BCV_I      ! BC loop counter
      INTEGER M, MM           ! Mass phase loop counter
      INTEGER RANGE_TOP, RANGE_BOT ! Dummy values
      INTEGER PHASE_CNT        ! Number of solid phases at bc
      INTEGER PHASE_LIST(dim_m) ! List of phases used in current bc

! the number of particles injected in a solids time step
      DOUBLE PRECISION NPMpSEC(dim_m) ! For solid phase m
      DOUBLE PRECISION NPpSEC
      DOUBLE PRECISION NPpDT        ! Total for BC
      DOUBLE PRECISION SCALED_VAL
      DOUBLE PRECISION MAX_DIA ! Max diameter of incoming particles at bc

      DOUBLE PRECISION :: EPs_ERR
      DOUBLE PRECISION :: VOL_FLOW

      LOGICAL, parameter :: setDBG = .false.
      LOGICAL :: dFlag
      LOGICAL :: FATAL

! Temp inlet velocity for solids phase M
      DOUBLE PRECISION VEL_TMP(dim_m)
      DOUBLE PRECISION EPs_TMP(dim_m)

! Minimum/maximum solids velocity at inlet.  Also used in the iterative
! steps as the starting and ending velocities
      DOUBLE PRECISION  MAX_VEL

      DOUBLE PRECISION  MINIPV, MAXIPV
      INTEGER :: OCCUPANTS

!......................................................................!

      CALL init_err_msg("SET_BC_DEM_MI")

      dflag = (dmp_log .AND. setdbg)
      if(dflag) write(*,"(2/,2x,'DEM inlet count: ',I4)") dem_bcmi

! Allocate the MI data structures for NEW and RES2 cases. Allocation for
! RES1 cases is done prior to reading the RES file.
      IF(run_type /= 'RESTART_1') CALL allocate_dem_mi

! Loop over BCs that flagged for DEM mass inflow.
      DO bcv_i = 1, dem_bcmi

! Get the user defined BC ID.
         bcv = dem_bcmi_map(bcv_i)

         if(dflag) write(*,"(2/,'Setting DEM_MI:',I3)") bcv_i

! The number of mass phases at this inlet.  While a system may be
! polydisperse, the inlet could consist of a single mass phase
         phase_cnt = 0
! The mass phase indices of incoming particles at this inlet
         phase_list(:) = -1
! The max diameter of incoming particles at this inlet
         max_dia = zero

! Determine if the inlet is mono or polydisperse
         DO m=1, mmax + des_mmax
            IF(solids_model(m) /= 'DEM') cycle
            IF(bc_rop_s(bcv,m) == undefined) cycle
            IF(compare(bc_rop_s(bcv,m),zero)) cycle
            phase_cnt = phase_cnt + 1
            phase_list(phase_cnt) = m
            max_dia = max(max_dia,d_p0(m))
         ENDDO
! Set the polydispersity flag.
         dem_mi(bcv_i)%POLYDISPERSE = (phase_cnt > 1)

! Layout the feed pattern.
         CALL layout_mi_dem(bcv, bcv_i, max_dia)

! Initialize
         max_vel = zero
         npmpsec(:) = zero

! Calculate the individual velocities for each solid phase
         DO mm = 1, phase_cnt
            m = phase_list(mm)

! Pull off the BC velocity normal to the flow plane.
            SELECT CASE(bc_plane(bcv))
            CASE('N','S'); vel_tmp(m) = abs(bc_v_s(bcv,m))
            CASE('E','W'); vel_tmp(m) = abs(bc_u_s(bcv,m))
            CASE('T','B'); vel_tmp(m) = abs(bc_w_s(bcv,m))
            END SELECT

! Check for min/max inlet velocity
            max_vel = max(abs(vel_tmp(m)), max_vel)

! Calculate the number of particles of mass phase M are injected per
! second for each solid phase present at the boundary.

! Use the mass flow rate if defined.
            IF(bc_massflow_s(bcv,m) == zero) THEN
               npmpsec(m) = zero

            ELSEIF(bc_massflow_s(bcv,m) /= undefined) THEN
               npmpsec(m) = bc_massflow_s(bcv,m) / &
                  (ro_s0(m)*(pi/6.d0*d_p0(m)**3))

! Otherwise use the volumetric flow rate if defined.
            ELSEIF(bc_volflow_s(bcv,m) == zero) THEN
               npmpsec(m) = zero

            ELSEIF(bc_volflow_s(bcv,m) /= undefined) THEN
               npmpsec(m) = bc_volflow_s(bcv,m) / (pi/6.d0*d_p0(m)**3)

! Calculate volumetric flow rate to convert to particle count. BC_AREA
! was already corrected for cut cells and velocity was recalculated
! to ensure user-specified mass or volumetric flow rates.
            ELSE
               vol_flow = vel_tmp(m) * bc_area(bcv) * bc_ep_s(bcv,m)
! Calculate the number of particles of mass phase M are injected per
! second for each solid phase present at the boundary
               npmpsec(m) = vol_flow / (pi/6.d0*d_p0(m)**3)
! Write some debugging information if needed.
            ENDIF
            if(dflag) write(*,1100) m, vel_tmp(m), npmpsec(m)
         ENDDO

! Total number of particles at BCV injected per second
         nppsec = sum(npmpsec)

 1100 FORMAT(/2x,'Conversion Info: Phase ',i2,/4x,'Velocity: ',g12.5,/ &
         4x,'NPMpSEC = ',f11.1)

         if(dflag) write(*,"(/2x,'Max Velocity:',3x,g12.5)") max_vel
         if(dflag) write(*,"( 2x,'NPpSEC:',3x,F11.1)") nppsec

! The number of total particles per solid time step DTSOLID
         nppdt = nppsec * dtsolid

! Inject one particle every PI_FACTOR solids time steps.
         IF(nppdt == zero)THEN
            pi_count(bcv_i) = 0
            pi_factor(bcv_i) = undefined_i
         ELSEIF(nppdt .LT. 1.0)THEN
            pi_count(bcv_i) = 1
            pi_factor(bcv_i) = floor(real(1.d0/nppdt))
! Inject PI_COUNT particles every soilds time step.
         ELSE
            pi_count(bcv_i) = ceiling(real(nppdt))
            pi_factor(bcv_i) = 1
         ENDIF

         occupants = dem_mi(bcv_i)%OCCUPANTS

         IF(pi_count(bcv_i) > 0) THEN
! Calculate the minimum inlet velocity. The cutoff is associated with
! square packing of disks on a plane.
            minipv = max_dia/( dtsolid*dble(pi_factor(bcv_i)) * dble(  &
               floor( real(occupants)/real(pi_count(bcv_i)))))
! Calculate the velocity needed to ensure that half the inlet is free.
! Inlets with velocities greater than this value can be randomly seeded,
! otherwise, particles are seeded in according to the grid.
            maxipv = max_dia/( dtsolid*dble(pi_factor(bcv_i)) * dble(  &
             floor(ceiling(real(occupants)/2.0)/real(PI_COUNT(bcv_i)))))
         ELSE
            minipv = -undefined
            maxipv =  undefined
         ENDIF

         if(dflag) write(*,"(/2x,'MaxIPV:',3x,g12.5)") maxipv
         if(dflag) write(*,"( 2x,'MinIPV:',3x,g12.5)") minipv

         IF(max_vel < minipv) THEN
            WRITE(err_msg,1110) bcv, max_vel, minipv
            CALL flush_err_msg(abort=.true.)
         ENDIF

 1110 FORMAT('Error 1110: Solids velocity for BC ',i3,' is too low ', &
         'to satisfy DEM',/'inlet packing restrictions. Potential ',  &
         'solutions:',//,' > If the BC velocities (BC_U_s, BC_V_s, ', &
         'BC_W_s) are defined, specify',/3x,'a larger value for the ',&
         'velocity normal to the flow plane.',//' > If MASSFLOW or ', &
         'VOLFLOW are defined, decrease the solids volume',/3x,       &
         'fraction to increase solids velocity.',//2x,'Max user-',    &
         'specified BC velocity:   ',g12.5,/2x,'Minimum required ',   &
         'solids Velocity: ',g12.5)

! Set all BC solids velocities to the largest velocity and recalculate
! BC_EP_s to determine the magnitude of the change.
         eps_err = zero
         eps_tmp = zero
         DO mm = 1, phase_cnt
            m = phase_list(mm)

            IF(bc_massflow_s(bcv,m) == zero .OR. &
               bc_volflow_s(bcv,m) == zero) THEN
! If no solids, set the tmp EPs to zero and clear solid velocity
               eps_tmp(m) = 0.0d0
               SELECT CASE(bc_plane(bcv))
               CASE('N','S'); bc_v_s(bcv,m) =  zero
               CASE('E','W'); bc_u_s(bcv,m) =  zero
               CASE('T','B'); bc_w_s(bcv,m) =  zero
               END SELECT

            ELSE
               eps_tmp(m) = bc_ep_s(bcv,m) * (vel_tmp(m) / max_vel)
! Over-write the current BC value.
               SELECT CASE(bc_plane(bcv))
               CASE('N'); bc_v_s(bcv,m) =  abs(max_vel)
               CASE('S'); bc_v_s(bcv,m) = -abs(max_vel)
               CASE('E'); bc_u_s(bcv,m) =  abs(max_vel)
               CASE('W'); bc_u_s(bcv,m) = -abs(max_vel)
               CASE('T'); bc_w_s(bcv,m) =  abs(max_vel)
               CASE('B'); bc_w_s(bcv,m) = -abs(max_vel)
               END SELECT

            ENDIF
            eps_err = eps_err + (bc_ep_s(bcv,m) - eps_tmp(m))
         ENDDO

! If the net change in solids volume fraction is greatere than 0.01,
! flag this as an error and exit. >> Let the user fix the input.
         IF(.NOT.compare(eps_err,small_number)) THEN
            IF(eps_err > 0.01) THEN
               WRITE(err_msg,1200) bcv, max_vel, eps_err
               CALL flush_err_msg(footer=.false.)
               fatal = .true.

! Report the amount of changes imposed on the BC in setting a
! uniform inlet velocity.
            ELSE
               WRITE(err_msg,1205) bcv, max_vel, eps_err
               CALL flush_err_msg(footer=.false.)
               fatal = .false.
            ENDIF

            WRITE(err_msg, 1210)
            CALL flush_err_msg(header=.false., footer=.false.)
            DO mm = 1, phase_cnt
            m = phase_list(mm)
               WRITE(err_msg,1211) m, bc_ep_s(bcv,m), eps_tmp(m), &
                  (bc_ep_s(bcv,m)-eps_tmp(m))
               CALL flush_err_msg(header=.false., footer=.false.)
            ENDDO
            CALL flush_err_msg(header=.false., abort=fatal)
         ENDIF


 1200 FORMAT('Error 1200: Unable to impose a uniform solids velocity ',&
         'on BC',i3,'.',/'Setting all solids to the highest velocity ',&
         'results in a large error',/'in the solids volume fraction. ',&
         'Please correct the mfix.dat file.',2/5x,'Max Inlet Velocity',&
         ':',1x,es11.4,/5x,'Total BC_EP_s Error:',es11.4)


 1205 FORMAT('Warning 1201: Uniform solids velocity imposed on BC',i3, &
         '.',2/,5x,'Uniform Inlet Velocity:',1x,es11.4,/5x,'Total ',   &
         'BC_EP_s Error:',4x,es11.4,/' ')

 1210 FORMAT(/,5x,'|',11('-'),'|',3(14('-'),'|'),/5x,'|',3x,'Phase',3x,&
         '|',4x,'BC_EP_s',3x,'|',2x,'Calculated',2x,'|',3x,'ABS ',   &
         'Error',2x,'|',/5x,'|',11('-'),'|',3(14('-'),'|'))

 1211 FORMAT(5x,'|',4x,i2,5x,'| ',1x,es11.4,1x,'|',1x,es11.4,2x,'|',   &
         1x,es11.4,1x,' |',/5x,'|',11('-'),'|',3(14('-'),'|'))


! For polydisperse inlets, construct the DES_POLY_LAYOUT array
         IF(dem_mi(bcv_i)%POLYDISPERSE) THEN
            range_bot = 1
            DO mm=1,phase_cnt - 1
               m = phase_list(mm)
               scaled_val = dble(numfrac_limit)*(npmpsec(m)/nppsec)
               range_top = floor(scaled_val) + (range_bot-1)
               dem_bc_poly_layout(bcv_i,range_bot:range_top) = m
               range_bot = range_top+1
            ENDDO

            m = phase_list(phase_cnt)
            dem_bc_poly_layout(bcv_i,range_bot:numfrac_limit) = m
! For monodisperse inlets, store the single mass phase used
         ELSE
            dem_bc_poly_layout(bcv_i,:) = phase_list(1)
         ENDIF


! Calculate des mass inlet time; time between injection.  If the run
! type is RESTART_1, DES_MI_TIME will be picked up from the restart file
! with an updated value.
         IF(run_type /= 'RESTART_1') &
            dem_mi_time(bcv_i) = time + dble(pi_factor(bcv_i)) * dtsolid


         WRITE(err_msg,1000) bcv, nppdt, int(nppdt/dtsolid), &
            pi_factor(bcv_i), pi_count(bcv_i), dem_mi_time(bcv_i)
         CALL flush_err_msg(header=.false., footer=.false.)

 1000 FORMAT(2/,2x,'For mass inlet BC: ', i3,/,&
         4x,'No. particles injected per solids time step = ', es15.8,/,&
         4x,'No. particles injected per second = ', i10,/,&
         4x,'PI_FACTOR = ', i10,' PI_COUNT = ', i5,/,&
         4x,'start DES_MI_TIME = ', es15.8,/'  ')

      ENDDO

!
      CALL set_dem_bcmi_ijk


      CALL finl_err_msg


      RETURN
      END SUBROUTINE set_bc_dem_mi


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SET_BC_DEM                                              !
!  Author: J.Musser                                   Date: 13-Jul-09  !
!                                                                      !
!  Purpose: Check the data provided for the des mass inflow boundary   !
!  condition and flag errors if the data is improper.  This module is  !
!  also used to convert the proveded information into the format       !
!  necessary for the dependent subrountines to function properly.      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE set_dem_bcmi_ijk

! Modules
!---------------------------------------------------------------------//
      use bc, only: bc_plane
      use bc, only: bc_x_w, bc_x_e, bc_y_s, bc_y_n, bc_z_b, bc_z_t
      use des_bc, only: dem_bcmi, dem_bcmi_map, dem_bcmi_ijk
      use des_bc, only: dem_bcmi_ijkstart, dem_bcmi_ijkend
      use desgrid, only: dg_imin1, dg_imax1
      use desgrid, only: dg_jmin1, dg_jmax1
      use desgrid, only: dg_kmin1, dg_kmax1
      use desgrid, only: dg_funijk
      use desgrid, only: iofpos, jofpos, kofpos
      use desgrid, only: dg_is_on_mype_plus1layers
      use error_manager
      use funits, only: dmp_log
      use mpi_utility
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
      INTEGER, ALLOCATABLE :: LOC_DEM_BCMI_IJK(:)
      INTEGER :: BCV, BCV_I
      INTEGER :: LC
      INTEGER :: MAX_CELLS
      INTEGER :: BND1, BND2
      LOGICAL, parameter :: setDBG = .false.
      LOGICAL :: dFlag
      INTEGER :: I,J,K,IJK
      INTEGER :: I_w, I_e, J_s, J_n, K_b, K_t
!......................................................................!

      CALL init_err_msg("SET_DEM_BCMI_IJK")

      dflag = (dmp_log .AND. setdbg)

      if(dflag) write(*,"(2/,2x,'From: SET_DEM_BCMI_IJK')")

! Loop over all inflow BCs to get an approximate count of the number
! of fluid cells that are adjacent to them.
      max_cells = 0
      DO bcv_i=1, dem_bcmi
         bcv = dem_bcmi_map(bcv_i)

! Set the search area a little bigger than the inlet area.
         if(dflag) WRITE(*,"(/2x,'Adding cells for BCV: ',I3)") bcv
         SELECT CASE (bc_plane(bcv))
         CASE('N','S')
            bnd1 = min(iofpos(bc_x_e(bcv))+1,dg_imax1) - &
               max(iofpos(bc_x_w(bcv))-1,dg_imin1)
            bnd2 = min(kofpos(bc_z_t(bcv))+1,dg_kmax1) - &
               max(kofpos(bc_z_b(bcv))-1,dg_kmin1)

         CASE('E','W')
            bnd1 = min(jofpos(bc_y_n(bcv))+1,dg_jmax1) - &
               max(jofpos(bc_y_s(bcv))-1,dg_jmin1)
            bnd2 = min(kofpos(bc_z_t(bcv))+1,dg_kmax1) - &
               max(kofpos(bc_z_b(bcv))-1,dg_kmin1)

         CASE('T','B')
            bnd1 = min(iofpos(bc_x_e(bcv))+1,dg_imax1) - &
               max(iofpos(bc_x_w(bcv))-1,dg_imin1)
            bnd2 = min(jofpos(bc_y_n(bcv))+1,dg_jmax1) - &
               max(jofpos(bc_y_s(bcv))-1,dg_jmin1)
         END SELECT

         max_cells = max_cells + (bnd1 + 1)*(bnd2 + 1)
         if(dflag) WRITE(*,"(4x,'Plane:   ',A)") bc_plane(bcv)
         if(dflag) WRITE(*,"(4x,'Cells: ', I8)") (bnd1+1)*(bnd2+1)
      ENDDO

      if(dflag) write(*,"(2x,'Max Cells: ',I8)") max_cells

! Allocate an array to hold the IJK values. This array should be
! more than enough to store all the IJKs.
      allocate( loc_dem_bcmi_ijk(max_cells) )


! Loop over the IJKs for each BC and store only the IJKs that you
! own as well as the start/end array positions for each BC.
      lc = 1
      DO bcv_i = 1, dem_bcmi

         dem_bcmi_ijkstart(bcv_i) = lc
         bcv = dem_bcmi_map(bcv_i)

         if(dflag) write(*,"(/2x,'Searching for fluid cells:',I3)") bcv

         i_w = max(iofpos(bc_x_w(bcv))-1,dg_imin1)
         i_e = min(iofpos(bc_x_e(bcv))+1,dg_imax1)
         j_s = max(jofpos(bc_y_s(bcv))-1,dg_jmin1)
         j_n = min(jofpos(bc_y_n(bcv))+1,dg_jmax1)
         k_b = max(kofpos(bc_z_b(bcv))-1,dg_kmin1)
         k_t = min(kofpos(bc_z_t(bcv))+1,dg_kmax1)

! Depending on the flow plane, the 'common' index needs set to reference
! the fluid cell.
         SELECT CASE (bc_plane(bcv))
         CASE('N')
            j_s = jofpos(bc_y_s(bcv))
            j_n = jofpos(bc_y_n(bcv))
         CASE('S')
            j_s = jofpos(bc_y_s(bcv))
            j_n = jofpos(bc_y_n(bcv))
         CASE('E')
            i_w = iofpos(bc_x_w(bcv))
            i_e = iofpos(bc_x_e(bcv))
         CASE('W')
            i_w = iofpos(bc_x_w(bcv))
            i_e = iofpos(bc_x_e(bcv))
         CASE('T')
            k_b = kofpos(bc_z_b(bcv))
            k_t = kofpos(bc_z_t(bcv))
         CASE('B')
            k_b = kofpos(bc_z_b(bcv))
            k_t = kofpos(bc_z_t(bcv))
         END SELECT

         if(dflag) then
            write(*,"(4x,'Search bounds: ')")
            write(*,"(6x,'I_w/I_e:',2(2x,I6))") i_w, i_e
            write(*,"(6x,'J_s/J_n:',2(2x,I6))") j_s, j_n
            write(*,"(6x,'K_b/K_t:',2(2x,I6))") k_b, k_t
         endif

! Store the IJKs.
         DO k = k_b, k_t
         DO j = j_s, j_n
         DO i = i_w, i_e
! Skip cells that this rank does not own or are considered dead.
            IF(.NOT.dg_is_on_mype_plus1layers(i,j,k))cycle

            ijk = dg_funijk(i,j,k)
            loc_dem_bcmi_ijk(lc) = ijk
            lc = lc+1
         ENDDO
         ENDDO
         ENDDO

         dem_bcmi_ijkend(bcv_i) = lc-1

         IF(dflag) write(*,1111) bcv, bcv_i,                           &
            dem_bcmi_ijkstart(bcv_i), dem_bcmi_ijkend(bcv_i)

      ENDDO

 1111 FORMAT(/2x,'DEM Mass Inflow:',/4x,'BC:',i4,3x,'MAP:',i4,         &
         /4x,'IJKSTART:',i6,/4x,'IJKEND:  ',i6)


! Allocate the global store arrary array. This changes across MPI ranks.
      IF(lc > 1) THEN
         allocate( dem_bcmi_ijk(lc-1) )
         dem_bcmi_ijk(1:lc-1) = loc_dem_bcmi_ijk(1:lc-1)
      ELSE
         allocate( dem_bcmi_ijk(1) )
         dem_bcmi_ijk(1) = loc_dem_bcmi_ijk(1)
      ENDIF

      deallocate(loc_dem_bcmi_ijk)

      CALL finl_err_msg

      RETURN
      END SUBROUTINE set_dem_bcmi_ijk