cfassign.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CFASSIGN                                                C
!                                                                      C
!  Purpose: Assign the necessary values for DEM computation. For       C
!           example:                                                   C
!     - assigning DEM boundaries from the values entered for           C
!       MFIX input in mfix.datat                                       C
!     - assigning DEM gravity vector from MFIX input.                  C
!     - calculating DTSOLID based on particle properties: spring       C
!       coefficient, damping factor & mass                             C
!                                                                      C
!  Reviewer: Sreekanth Pannala                        Date: 09-Nov-06  C
!                                                                      C
!  Reviewer: Rahul Garg                               Date: 25-Mar-14  C
!  Comments: Breaking this subroutine into several subroutines for DEM C
!            and PIC models                                            C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE cfassign

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1
      USE constant, only : gravity_x, gravity_y, gravity_z
      USE discretelement
      USE mfix_pic
      use error_manager
! Flag: DEM solids present.
      use run, only: dem_solids
! Runtime flag specifying MPPIC solids
      use run, only: pic_solids

      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------

      CALL init_err_msg("CFASSIGN")

! The common assignments are done in this routine.
! The model spcific assignmets are moved to the specific subroutines

      grav(1) = gravity_x
      grav(2) = gravity_y
      grav(3) = gravity_z
      grav_mag = sqrt(dot_product(grav,grav))

      IF(dem_solids) CALL cfassign_dem
      IF(pic_solids) CALL cfassign_pic

! Finalize the error manager.
      CALL finl_err_msg

      RETURN
      END SUBROUTINE cfassign



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CFASSIGN_PIC                                            C
!                                                                      C
!  Purpose: PIC related cfassign source code moved here                C
!           example:                                                   C
!     - calculating DTSOLID based on particle response time            C
!                                                                      C
!                                                                      C
!  Reviewer: Rahul Garg                               Date: 25-Mar-14  C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE cfassign_pic


!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE discretelement, only: des_mmax, dtsolid
      use param1, only: large_number
      USE physprop, only: mu_g0, mmax, d_p0, ro_s0
      USE mfix_pic, only : dtpic_taup, des_tau_p
      USE mfix_pic, only: mppic_pdrag_implicit
      use error_manager
      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
      INTEGER :: M, MMAX_TOT
!-----------------------------------------------

      CALL init_err_msg("CFASSIGN_PIC")

      mmax_tot = des_mmax+mmax
      DO m = mmax+1, mmax_tot
! account for a possible offset index when using d_p0 and ro_s.
         des_tau_p(m) = ro_s0(m)*(d_p0(m)**2.d0)/(18.d0*mu_g0)
         WRITE(err_msg,'(/A,I2,A,G17.8)') &
         'TAU_P FOR ', m,'th SOLID PHASE= ', des_tau_p(m)
         CALL flush_err_msg (header = .false., footer = .false.)

      ENDDO

      dtsolid = minval(des_tau_p(mmax+1:mmax_tot))
      dtpic_taup = dtsolid      !maximum dt for point-particles based on taup
      if(mppic_pdrag_implicit) dtpic_taup = large_number

      WRITE(err_msg, 1000) dtsolid
      CALL flush_err_msg(header = .false.)

 1000 format('MPPIC: Point-particle ',&
      'approximation for particle-particle and ', /, &
      'particle-wall collisions', /, &
      'DTSOLID based on particle time response Taup', /, &
      'DTSOLID = ', 2x, e17.10)

! Finalize the error manager.
      CALL finl_err_msg

      END SUBROUTINE cfassign_pic




!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CFASSIGN_DEM                                            C
!                                                                      C
!  Purpose: DEM related cfassign source code moved here                C
!           example:                                                   C
!     - calculating DTSOLID based on particle properties: spring       C
!       coefficient, damping factor & mass                             C
!                                                                      C
!                                                                      C
!  Reviewer: Rahul Garg                               Date: 25-Mar-14  C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE cfassign_dem


!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1
      USE discretelement
      use error_manager
      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------

      CALL init_err_msg("CFASSIGN_DEM")

!      WRITE(ERR_MSG,'(A)') 'Setting collision model parameters for DEM'
!      CALL FLUSH_ERR_MSG (Footer = .false.)

! Finalize the error manager.
      CALL finl_err_msg
      END SUBROUTINE cfassign_dem

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
! subroutine: compute_volume_of_nodes                                      C
! Author: Rahul Garg                                                       C
! Dare: Dec 20, 2011                                                       C
! Purpose:  Due to the staggered grid, the interpolaiton of mean fields    C
! is always first done at the nodes (of the scalar cell) for any quantity. C
! For a quantity like drag force or ep_s, one needs to assign a geometric  C
! volume to a node. In the past this was done on-the-fly in drag_fgs.f.    C
! VCELL was the variable that was used and it used to be incorrecty set to C
! the volume of the scalar cell that the particle belongs to.              C
! This will be ok for uniform grid and will not hold for non-uniform grids.C
! In addition, at the edges (in 2-D, the node volume is only half of       C
! the standard cell volume. And for corner (in 2-D), it is only one fourth.C
! This was also done on-the-fly earlier                                    C
! But again the volume of the cell was used, which was not correct         C
! also not extendable to cut-cell. So this routine computes the geoemetric C
! volume of the nodes.                                                     C
!                                                                          C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE compute_volume_of_nodes

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE physprop
      USE fldvar
      USE run
      USE geometry
      USE indices
      USE bc
      USE compar
      USE sendrecv
      USE discretelement
      USE cutcell
      USE functions
      implicit none
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! ijk index of fluid grid and corresponding i, j, k indices
      integer :: ijk, iraw, jraw, kraw
! i, j, k and (i+1, j+1, k+1) indices corrected for any
! cyclic ghost boundaries on fluid grid
      INTEGER :: I, J, K, ip, jp, kp
      integer :: ipjk, ijpk, ipjpk, ijkp, ipjkp, ijpkp, ipjpkp
! volume of indicated grid
      double precision :: vol_ijk, vol_ipjk, vol_ijpk, vol_ipjpk
      double precision :: vol_ijkp, vol_ipjkp, vol_ijpkp, vol_ipjpkp
      double precision :: vol_node_count, vol_node_actual_count
! weighting factor
      double precision :: avg_factor
! not used?
      double precision :: vol_node_uncorr
!-----------------------------------------------

      avg_factor = merge(0.25d0, 0.125d0, no_k)

! compute the volume at the grid nodes
! grid nodes start from istart2 to iend1
      vol_node_count = merge(4., 8., no_k)

      DO ijk = ijkstart3,ijkend3
         des_vol_node(ijk) = zero
         iraw = i_of(ijk)
         jraw = j_of(ijk)
         kraw = k_of(ijk)



! start at 1 (ghost cell) and go to last fluid cell. why start at a
! ghost cell and not a fluid cell?
! See below

! Since we are interested in nodes making up the interpolation stencil,
! their numbering goes from 1 to iend1.
! Think of a case with IMAX = 3. Here the nodes where the interpolation will be
! done will run from 1 (=istart2) to 4 (=iend1).
         IF(iraw.LT.istart2 .OR. iraw.GT.iend1) cycle
         IF(jraw.LT.jstart2 .OR. jraw.GT.jend1) cycle
         IF(kraw.LT.kstart2 .OR. kraw.GT.kend1) cycle

! this will force indices of ghost cells on cyclic borders back to
! the corresponding fluid cells. since we are using i, j, k indices and
! not just a composite ijk index we need these to be shifted to account
! for periodicity
         i = imap_c(iraw)
         j = jmap_c(jraw)
         k = kmap_c(kraw)
         ip = imap_c(iraw+1)
         jp = jmap_c(jraw+1)

! using a function like ip_of(ijk) should work the same as getting funijk
! of the shifted i, j, k indices.  however, small differences will
! occur on the 'edges/corners'. so retaining the latter method at this
! time. see j. galvin for discussion
         ipjk = funijk(ip,j,k)
         ijpk = funijk(i,jp,k)
         ipjpk = funijk(ip,jp,k)

! the existing variable vol(ijk) is not used here for cut-cell reasons
! see r. garg for discussion
         vol_ijk   = dx(i) *dy(j) *dz(k)
         vol_ipjk  = dx(ip)*dy(j) *dz(k)
         vol_ijpk  = dx(i) *dy(jp)*dz(k)
         vol_ipjpk = dx(ip)*dy(jp)*dz(k)

         vol_node_uncorr = avg_factor*(vol_ijk + vol_ipjk + vol_ijpk + vol_ipjpk)
         vol_node_actual_count = vol_node_count

         IF(.NOT.fluid_at(ijk)) THEN
            vol_node_actual_count = vol_node_actual_count - 1
            vol_ijk  = zero
         ENDIF

         IF(.NOT.fluid_at(ipjk)) THEN
            vol_node_actual_count = vol_node_actual_count - 1
            vol_ipjk  = zero
         ENDIF

         IF(.NOT.fluid_at(ijpk)) THEN
            vol_node_actual_count = vol_node_actual_count - 1
            vol_ijpk  = zero
         ENDIF

         IF(.NOT.fluid_at(ipjpk)) THEN
            vol_node_actual_count = vol_node_actual_count - 1
            vol_ipjpk = zero
         ENDIF

! this will have non-zero values for non-fluid cells at the
! west/south/bottom borders but not for east/north/top borders?
         des_vol_node(ijk) = avg_factor*(vol_ijk + vol_ipjk + &
            vol_ijpk + vol_ipjpk)

         IF (do_k) THEN
            kp     = kmap_c(kraw+1)
            ijkp   = funijk(i, j, kp)
            ipjkp  = funijk(ip,j, kp)
            ijpkp  = funijk(i, jp,kp)
            ipjpkp = funijk(ip,jp,kp)

            vol_ijkp   = dx(i) *dy(j) *dz(kp)
            vol_ipjkp  = dx(ip)*dy(j) *dz(kp)
            vol_ijpkp  = dx(i) *dy(jp)*dz(kp)
            vol_ipjpkp = dx(ip)*dy(jp)*dz(kp)

            vol_node_uncorr = avg_factor*(vol_node_uncorr + vol_ijkp + &
               vol_ipjkp + vol_ijpkp + vol_ipjpkp)

            IF(.NOT.fluid_at(ijkp)) THEN
               vol_node_actual_count = vol_node_actual_count - 1
               vol_ijkp   = zero
            ENDIF

            IF(.NOT.fluid_at(ipjkp)) THEN
               vol_node_actual_count = vol_node_actual_count - 1
               vol_ipjkp  = zero
            ENDIF

            IF(.NOT.fluid_at(ijpkp)) THEN
               vol_node_actual_count = vol_node_actual_count - 1
               vol_ijpkp  = zero
            ENDIF

            IF(.NOT.fluid_at(ipjpkp)) THEN
               vol_node_actual_count = vol_node_actual_count - 1
               vol_ipjpkp = zero
            ENDIF

            des_vol_node(ijk) = des_vol_node(ijk) + avg_factor*&
               (vol_ijkp + vol_ipjpkp + vol_ijpkp + vol_ipjkp)

         ENDIF


      ENDDO   ! end do ijk=ijkstart3,ijkend3

      RETURN
      RETURN
      END SUBROUTINE compute_volume_of_nodes