mass_inflow_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: MPPIC_MI_BC                                             C
!  Purpose:                                                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE mass_inflow_pic

! Modules
!---------------------------------------------------------------------//
      USE bc
      USE constant
      USE cutcell
      USE discretelement
      use error_manager
      USE functions
      USE geometry
      USE mfix_pic
      USE mpi_utility
      USE param, only: dimension_m
      USE param1, only: half, zero
      use physprop, only: mmax, d_p0, ro_s0
      USE pic_bc
      USE randomno
      USE run
      IMPLICIT NONE

! Local Variables
!---------------------------------------------------------------------//
      INTEGER :: WDIR, IDIM, IPCOUNT, LAST_EMPTY_SPOT, NEW_SPOT
      INTEGER :: BCV, BCV_I, L, LC, PIP_ADD_COUNT, IPROC
      INTEGER ::  IFLUID, JFLUID, KFLUID, IJK, M
      DOUBLE PRECISION :: CORD_START(3), DOML(3),WALL_NORM(3)
      DOUBLE PRECISION :: AREA_INFLOW, VEL_INFLOW(dimn), STAT_WT

      LOGICAL :: DELETE_PART
      DOUBLE PRECISION :: EPS_INFLOW(dimension_m)
      DOUBLE PRECISION :: REAL_PARTS(dimension_m)
      DOUBLE PRECISION :: COMP_PARTS(dimension_m)
      DOUBLE PRECISION :: VEL_NORM_MAG, VOL_INFLOW, VOL_IJK

! Temp logical variables for checking constant npc and statwt specification
      LOGICAL :: CONST_NPC, CONST_STATWT

      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RANDPOS
      INTEGER :: CNP_CELL_COUNT
      DOUBLE PRECISION :: RNP_CELL_COUNT
      integer :: lglobal_id
      integer, dimension(0:numpes-1) :: add_count_all

      LOGICAL, parameter :: setDBG = .true.
      LOGICAL :: dFlag

      type :: ty_spotlist
         integer spot
         type(ty_spotlist),pointer :: next => null()
      end type ty_spotlist

      type(ty_spotlist),pointer :: &
           cur_spotlist => null(), &
           prev_spotlist => null(), &
           temp_spotlist => null()
!......................................................................!


      CALL init_err_msg("PIC_MI_BC")

      dflag = (dmp_log .AND. setdbg)
      pip_add_count = 0
      last_empty_spot = 0

      ALLOCATE(cur_spotlist); cur_spotlist%spot = -1

      DO bcv_i = 1, pic_bcmi
         bcv = pic_bcmi_map(bcv_i)

         wall_norm(1:3) =  pic_bcmi_normdir(bcv_i,1:3)

         !Find the direction of the normal for this wall
         wdir = 0
         DO idim = 1, dimn
            wdir = wdir + abs(wall_norm(idim))*idim
         end DO

         DO lc=pic_bcmi_ijkstart(bcv_i), pic_bcmi_ijkend(bcv_i)
            ijk = pic_bcmi_ijk(lc)

            IF(.NOT.fluid_at(ijk)) cycle

            ifluid = i_of(ijk)
            jfluid = j_of(ijk)
            kfluid = k_of(ijk)

            cord_start(1) = xe(ifluid) - pic_bcmi_offset(bcv_i,1)*dx(ifluid)

            cord_start(2) = yn(jfluid) - pic_bcmi_offset(bcv_i,2)*dy(jfluid)


            cord_start(3) = merge(zero, zt(kfluid) - pic_bcmi_offset(bcv_i,3)*dz(kfluid), no_k)

            doml(1) = dx(ifluid)
            doml(2) = dy(jfluid)
            doml(3) = merge(dz(1), dz(kfluid), no_k)

            area_inflow = doml(1)*doml(2)*doml(3)/doml(wdir)

            vol_ijk = doml(1)*doml(2)*doml(3)

            doml(wdir) = zero
            !set this to zero as the particles will
            !be seeded only on the BC plane

            DO m = 1, des_mmax+mmax

               IF(solids_model(m) /= 'PIC') cycle

               eps_inflow(m) = bc_rop_s(bcv, m)/ro_s0(m)
               vel_inflow(1) = bc_u_s(bcv, m)
               vel_inflow(2) = bc_v_s(bcv, m)
               vel_inflow(3) = bc_w_s(bcv, m)

               vel_norm_mag = abs(dot_product(vel_inflow(1:dimn), wall_norm(1:dimn)))
               vol_inflow   = area_inflow*vel_norm_mag*dtsolid

               real_parts(m) = 6.d0*eps_inflow(m)*vol_inflow/&
                               (pi*(d_p0(m)**3.d0))
               comp_parts(m) = zero

               const_npc    = (bc_pic_mi_const_npc  (bcv, m) .ne. 0)
               const_statwt = (bc_pic_mi_const_statwt(bcv, m) .ne. zero)
               IF(const_npc) THEN
                  IF(eps_inflow(m).GT.zero) &
                  comp_parts(m) = REAL(bc_pic_mi_const_npc(bcv, m))* &
                  VOL_INFLOW/VOL_IJK
               ELSEIF(const_statwt) THEN
                  comp_parts(m) = real_parts(m)/ &
                  bc_pic_mi_const_statwt(bcv, m)
               ENDIF

               pic_bcmi_rnp(lc,m) = pic_bcmi_rnp(lc,m) + real_parts(m)

               pic_bcmi_cnp(lc,m) = pic_bcmi_cnp(lc,m) + comp_parts(m)


               IF(pic_bcmi_cnp(lc,m).GE.1.d0) THEN
                  cnp_cell_count = int(pic_bcmi_cnp(lc,m))
                  pic_bcmi_cnp(lc,m) = pic_bcmi_cnp(lc,m) - cnp_cell_count

                  rnp_cell_count = pic_bcmi_rnp(lc,m)
                  pic_bcmi_rnp(lc,m)  = zero

                  !set pic_bcmi_rnp to zero to reflect that all real particles have been seeded
                  stat_wt = rnp_cell_count/REAL(cnp_cell_count)
                  ALLOCATE(randpos(cnp_cell_count*dimn))
                  CALL uni_rno(randpos(:))

                  DO ipcount = 1, cnp_cell_count


                     CALL pic_find_empty_spot(last_empty_spot, new_spot)

                     des_pos_old( new_spot,1:dimn) =  cord_start(1:dimn) &
                          & + randpos((ipcount-1)*dimn+1: &
                          & (ipcount-1)*dimn+dimn)*doml(1:dimn)
                     des_pos_new( new_spot,:) = des_pos_old(new_spot,:)
                     des_vel_old(new_spot,1:dimn) = vel_inflow(1:dimn)

                     des_vel_new( new_spot,:) = des_vel_old( new_spot,:)

                     des_radius(new_spot) = d_p0(m)*half

                     ro_sol(new_spot) =  ro_s0(m)

                     des_stat_wt(new_spot) = stat_wt

                     pijk(new_spot, 1) = ifluid
                     pijk(new_spot, 2) = jfluid
                     pijk(new_spot, 3) = kfluid
                     pijk(new_spot, 4) = ijk
                     pijk(new_spot, 5) = m

                     pvol(new_spot) = (4.0d0/3.0d0)*pi*des_radius(new_spot)**3
                     pmass(new_spot) = pvol(new_spot)*ro_sol(new_spot)

                     fc( new_spot,:) = zero
                     delete_part = .false.
                     IF(pic_bcmi_incl_cutcell(bcv_i)) &
                          CALL check_if_parcel_overlaps_stl &
                          (des_pos_new( new_spot,1:dimn), &
                          delete_part)

                     IF(.NOT.delete_part) THEN

                        pip = pip+1
                        pip_add_count = pip_add_count + 1
                        CALL set_normal(new_spot)
                        ! add to the list
                        ALLOCATE(temp_spotlist)
                        temp_spotlist%spot = new_spot
                        temp_spotlist%next => cur_spotlist
                        cur_spotlist => temp_spotlist
                        nullify(temp_spotlist)

                     ELSE
                        CALL set_nonexistent(new_spot)
                        last_empty_spot = new_spot - 1
                     ENDIF

                     !WRITE(*,'(A,2(2x,i5), 2x, A,2x, 3(2x,i2),2x, A, 3(2x,g17.8))') &
                     !   'NEW PART AT ', NEW_SPOT, MAX_PIP, 'I, J, K = ', IFLUID, JFLUID, KFLUID, 'POS =', DES_POS_NEW(NEW_SPOT,:)
                     !IF(DMP_LOG) WRITE(UNIT_LOG,'(A,2x,i5, 2x, A,2x, 3(2x,i2),2x, A, 3(2x,g17.8))') &
                     !    'NEW PART AT ', NEW_SPOT, 'I, J, K = ', IFLUID, JFLUID, KFLUID, 'POS =', DES_POS_NEW(NEW_SPOT,:)

                     !WRITE(*,*) 'WDIR, DOML = ', WDIR, DOML(:)
                  END DO
                  DEALLOCATE(randpos)

               end IF
            end DO
         end DO
      end DO

!Now assign global id to new particles added
      add_count_all(:) = 0
      add_count_all(mype) = pip_add_count
      call global_all_sum(add_count_all(0:numpes-1))
      lglobal_id = imax_global_id + sum(add_count_all(0:mype-1))

      do l = 1,pip_add_count
         lglobal_id = lglobal_id + 1
         iglobal_id(cur_spotlist%spot)= lglobal_id
         prev_spotlist=> cur_spotlist
         cur_spotlist => cur_spotlist%next
         deallocate(prev_spotlist)
      end do
      deallocate(cur_spotlist)
      imax_global_id = imax_global_id + sum(add_count_all(0:numpes-1))

      IF(pic_report_seeding_stats) then

         IF(sum(add_count_all(:)).GT.0) THEN
            WRITE(err_msg,'(/,2x,A,2x,i10)') &
               'TOTAL NUMBER OF PARCELS ADDED GLOBALLY = ',&
                sum(add_count_all(:))

            call flush_err_msg(header = .false., footer = .false.)

            DO iproc = 0, numpes-1
               IF(dmp_log) WRITE(unit_log, '(/,A,i4,2x,A,2x,i5)') &
                  'PARCELS ADDED ON PROC:', iproc,&
                  ' EQUAL TO', add_count_all(iproc)
            ENDDO
         ENDIF

      ENDIF

      CALL finl_err_msg
      RETURN
      END SUBROUTINE mass_inflow_pic



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PIC_FIND_EMPTY_SPOT                                     C
!  Purpose:                                                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE pic_find_empty_spot(LAST_INDEX, EMPTY_SPOT)

! Modules
!---------------------------------------------------------------------//
      USE funits
      USE mpi_utility
      USE error_manager
      USE discretelement, only: max_pip
      USE functions, only: is_nonexistent
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      INTEGER, INTENT(INOUT) :: LAST_INDEX
      INTEGER, INTENT(OUT) :: EMPTY_SPOT

! Local Variables
!---------------------------------------------------------------------//
      LOGICAL :: SPOT_FOUND
      INTEGER :: LL
!......................................................................!

      CALL init_err_msg("PIC_FIND_EMPTY_SPOT")

      IF(last_index.EQ.max_pip) THEN
         WRITE(err_msg,2001)

         CALL flush_err_msg(abort = .true.)
         call mfix_exit(mype)
      ENDIF
      spot_found = .false.

      DO ll = last_index+1, max_pip

         if(is_nonexistent(ll)) THEN
            empty_spot = ll
            last_index = ll
            spot_found = .true.
            EXIT
         ENDIF
      ENDDO

      IF(.NOT.spot_found) THEN
         WRITE(err_msg,2002)
         CALL flush_err_msg(abort = .true.)
      ENDIF

 2001 FORMAT(/,5x,  &
      & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5x, &
      & 'NO MORE EMPTY SPOT IN THE PARTICLE ARRAY TO ADD A NEW PARTICLE',/5x, &
      & 'TERMINAL ERROR: STOPPING')

 2002 FORMAT(/,5x,  &
      & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5x, &
      & 'COULD NOT FIND A SPOT FOR ADDING NEW PARTICLE',/5x, &
      & 'INCREASE THE SIZE OF THE INITIAL ARRAYS', 5x, &
      & 'TERMINAL ERROR: STOPPING')

      CALL finl_err_msg
      END SUBROUTINE pic_find_empty_spot





!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PIC_FIND_NEW_CELL                                     C
!  Purpose:                                                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE pic_find_new_cell(LL)

! Modules
!---------------------------------------------------------------------//
      USE discretelement
      use mpi_utility
      USE cutcell
      USE functions
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      INTEGER, INTENT(IN) :: LL

! Local Variables
!---------------------------------------------------------------------//
      DOUBLE PRECISION :: XPOS, YPOS, ZPOS
      INTEGER :: I, J, K
!......................................................................!

      i = pijk(ll,1)
      j = pijk(ll,2)
      k = pijk(ll,3)
      xpos = des_pos_new(ll,1)
      ypos = des_pos_new(ll,2)
      IF (dimn .EQ. 3) THEN
         zpos = des_pos_new(ll,3)
      ENDIF

      IF(xpos >= xe(i-1) .AND. xpos < xe(i)) THEN
         pijk(ll,1) = i
      ELSEIF(xpos >= xe(i)) THEN
         pijk(ll,1) = i+1
      ELSE
         pijk(ll,1) = i-1
      END IF

      IF(ypos >= yn(j-1) .AND. ypos < yn(j)) THEN
         pijk(ll,2) = j
      ELSEIF(ypos >= yn(j))THEN
         pijk(ll,2) = j+1
      ELSE
         pijk(ll,2) = j-1
      END IF

      IF(dimn.EQ.2) THEN
         pijk(ll,3) = 1
      ELSE
         IF(zpos >= zt(k-1) .AND. zpos < zt(k)) THEN
            pijk(ll,3) = k
         ELSEIF(zpos >= zt(k)) THEN
            pijk(ll,3) = k+1
         ELSE
            pijk(ll,3) = k-1
         END IF
      ENDIF

      i = pijk(ll,1)
      j = pijk(ll,2)
      k = pijk(ll,3)

      IF(i.EQ.iend1+1) then
         IF(xpos>=xe(iend1-1) .AND. xpos<=xe(iend1)) pijk(ll,1) = iend1
      ENDIF

      IF(j.EQ.jend1+1) then
         IF(ypos>=yn(jend1-1) .AND. ypos<=yn(jend1)) pijk(ll,2) = jend1
      ENDIF

      IF(dimn.EQ.3.AND.k.EQ.kend1+1) THEN
         IF(zpos>=zt(kend1-1) .AND. zpos<=zt(kend1)) pijk(ll,3) = kend1
      ENDIF

      pijk(ll,4) = funijk(pijk(ll,1), pijk(ll,2),pijk(ll,3))

      END SUBROUTINE pic_find_new_cell