gas_drag.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GAS_DRAG_W                                              !
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
!                                                                      !
!  Purpose: Account for the equal and opposite drag force on the gas   !
!           phase due to particles by introducing the drag as a        !
!           source term.  Face centered.                               !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE gas_drag_u(A_M, B_M, IER)

! Global Variables:
!---------------------------------------------------------------------//
! Runtime Flag: Interpolate DES values
      use particle_filter, only: des_interp_scheme_enum, des_interp_garg
! Flag: Gas sees the effect of particles in gas/solids flows.
      use discretelement, only: des_oneway_coupled
! Coefficient at cell corners added to the gas momentum A matix.
      use discretelement, only: drag_am, f_gds
! Coefficient at cell corners added to gas momentum B vector.
      use discretelement, only: drag_bm
! Volume of X-momentum cell
      use geometry, only: vol_u
! Flag to calculate Z direction
      use geometry, only: do_k

! Global Parameters:
!---------------------------------------------------------------------//
! Size of IJK arrays and solids phase..
      use param, only: dimension_3, dimension_m
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! The I, J, and K values that comprise an IJK
      use indices, only: i_of, j_of, k_of
! Flag: Fluid exists at indexed cell
      use functions, only: fluid_at
! IJK of cell to east.
      use functions, only: east_of
! IJK function for I,J,K that includes mapped indices.
      use compar, only: funijk_map_c
! Function for averaging to a scalar cell's east face.
      use fun_avg, only: avg_x
! Domain index bounds.
      use compar, only: istart2, jstart2, kstart2
      use compar, only: iend2, jend2, kend2

      IMPLICIT NONE

! Dummy Arguments:
!---------------------------------------------------------------------//
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_M(dimension_3,-3:3,0:dimension_m)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_M(dimension_3, 0:dimension_m)
! Error index
      INTEGER, INTENT(INOUT) :: IER


! Local Variables:
!---------------------------------------------------------------------//
! Grid cell indices
      INTEGER :: I, J, K, IJK, IJMK, IJKM, IJMKM, IJKE
! temporary variables for matrix A_M and vector B_M
      DOUBLE PRECISION :: tmp_A, tmp_B
! Averaging factor
      DOUBLE PRECISION :: AVG_FACTOR
!......................................................................!

! Initialize error flag.
      ier = 0

! Skip this routine if the gas/solids are only one-way coupled.
      IF(des_oneway_coupled) RETURN

! Average the interpoalted drag force from the cell corners to the cell face.
      IF(des_interp_scheme_enum == des_interp_garg)THEN

         avg_factor = merge(0.25d0, 0.5d0, do_k)

!!$omp parallel do schedule(guided,50) default(none)                    &
!!$omp shared(IJKSTART3, IJKEND3, FUNIJK_MAP_C, I_OF, J_OF,             &
!!$omp    K_OF, DO_K, AVG_FACTOR, DRAG_AM, DRAG_BM, A_M, B_M, VOL_U)    &
!!$omp private(IJK, I, J, K, IJMK, IJKM, IJMKM, tmp_A, tmp_B)
         DO ijk = ijkstart3, ijkend3

            IF(.NOT.fluid_at(ijk)) cycle

            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            IF (i.LT.istart2 .OR. i.GT.iend2) cycle
            IF (j.LT.jstart2 .OR. j.GT.jend2) cycle
            IF (k.LT.kstart2 .OR. k.GT.kend2) cycle

            ijmk = funijk_map_c(i, j-1, k)

            tmp_a = -avg_factor*(drag_am(ijk) + drag_am(ijmk))
            tmp_b = -avg_factor*(drag_bm(ijk,1) + drag_bm(ijmk,1))

            IF(do_k) THEN
               ijkm = funijk_map_c(i, j, k-1)
               ijmkm = funijk_map_c(i, j-1, k-1)
               tmp_a = tmp_a - avg_factor*                             &
                  (drag_am(ijkm) + drag_am(ijmkm))
               tmp_b = tmp_b - avg_factor*                             &
                  (drag_bm(ijkm,1) + drag_bm(ijmkm,1))
            ENDIF

            a_m(ijk,0,0) = a_m(ijk,0,0) + tmp_a*vol_u(ijk)
            b_m(ijk,0) = b_m(ijk,0) + tmp_b*vol_u(ijk)

         ENDDO
!!$omp end parallel do

      ELSE

!$omp parallel do default(none) &
!$omp private(IJK, I, IJKE, tmp_A, tmp_B) &
!$omp shared(IJKSTART3,IJKEND3,I_OF, F_GDS, DRAG_BM, A_M, B_M, VOL_U)
         DO ijk = ijkstart3, ijkend3
            IF(fluid_at(ijk)) THEN
               i = i_of(ijk)
               ijke = east_of(ijk)

               tmp_a = avg_x(f_gds(ijk), f_gds(ijke), i)
               tmp_b = avg_x(drag_bm(ijk,1), drag_bm(ijke,1), i)

               a_m(ijk,0,0) = a_m(ijk,0,0) - vol_u(ijk) * tmp_a
               b_m(ijk,0) = b_m(ijk,0) - vol_u(ijk) * tmp_b
            ENDIF
         ENDDO
!$omp end parallel do
      ENDIF

      END SUBROUTINE gas_drag_u

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GAS_DRAG_V                                              !
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
!                                                                      !
!  Purpose: Account for the equal and opposite drag force on the gas   !
!           phase due to particles by introducing the drag as a        !
!           source term.  Face centered.                               !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE gas_drag_v(A_M, B_M, IER)


! Global Variables:
!---------------------------------------------------------------------//
! Runtime Flag: Interpolate DES values
      use particle_filter, only: des_interp_scheme_enum, des_interp_garg
! Flag: Gas sees the effect of particles in gas/solids flows.
      use discretelement, only: des_oneway_coupled
! Coefficient at cell corners added to the gas momentum A matix.
      use discretelement, only: drag_am, f_gds
! Coefficient at cell corners added to gas momentum B vector.
      use discretelement, only: drag_bm
! Volume of Y-momentum cell
      use geometry, only: vol_v
! Flag to calculate Z direction
      use geometry, only: do_k

! Global Parameters:
!---------------------------------------------------------------------//
! Size of IJK arrays and solids phase..
      use param, only: dimension_3, dimension_m
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! The I, J, and K values that comprise an IJK
      use indices, only: i_of, j_of, k_of
! Flag: Fluid exists at indexed cell
      use functions, only: fluid_at
! IJK of cell to north.
      use functions, only: north_of
! IJK function for I,J,K that includes mapped indices.
      use compar, only: funijk_map_c
! Function for averaging to a scalar cell's north face.
      use fun_avg, only: avg_y
! Domain index bounds.
      use compar, only: istart2, jstart2, kstart2
      use compar, only: iend2, jend2, kend2


      IMPLICIT NONE

! Dummy Arguments:
!---------------------------------------------------------------------//
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_M(dimension_3, -3:3, 0:dimension_m)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_M(dimension_3, 0:dimension_m)
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local Variables:
!---------------------------------------------------------------------//
! Grid cell indices
      INTEGER :: I, J, K, IJK, IMJK, IJKM, IMJKM, IJKN
! temporary variables for matrix A_M and vector B_M
      DOUBLE PRECISION tmp_A, tmp_B
! Averaging factor
      DOUBLE PRECISION :: AVG_FACTOR
!......................................................................!

! Initialize error flag.
      ier = 0

! Skip this routine if the gas/solids are only one-way coupled.
      IF(des_oneway_coupled) RETURN

      IF(des_interp_scheme_enum == des_interp_garg)THEN

         avg_factor = merge(0.25d0, 0.5d0, do_k)

!!$omp parallel do schedule (guided,50) default(none)                   &
!!$omp shared(IJKSTART3, IJKEND3, FUNIJK_MAP_C, I_OF, J_OF,             &
!!$omp    K_OF, DO_K, AVG_FACTOR, DRAG_AM, DRAG_BM, A_M, B_M, VOL_V)    &
!!$omp private(IJK, I, J, K, IMJK, IJKM, IMJKM, tmp_A, tmp_B)
         DO ijk = ijkstart3, ijkend3
            IF(.NOT.fluid_at(ijk)) cycle

            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            IF (i.LT.istart2 .OR. i.GT.iend2) cycle
            IF (j.LT.jstart2 .OR. j.GT.jend2) cycle
            IF (k.LT.kstart2 .OR. k.GT.kend2) cycle

            imjk = funijk_map_c(i-1,j,k)

            tmp_a = -avg_factor*(drag_am(ijk) + drag_am(imjk))
            tmp_b = -avg_factor*(drag_bm(ijk,2) + drag_bm(imjk,2))

            IF(do_k) THEN

               ijkm = funijk_map_c(i,j,k-1)
               imjkm = funijk_map_c(i-1,j,k-1)

               tmp_a = tmp_a - avg_factor*                             &
                  (drag_am(ijkm) + drag_am(imjkm))
               tmp_b = tmp_b - avg_factor*                             &
                  (drag_bm(ijkm,2) + drag_bm(imjkm,2))
            ENDIF

            a_m(ijk,0,0) = a_m(ijk,0,0) + tmp_a*vol_v(ijk)
            b_m(ijk,0) = b_m(ijk,0) + tmp_b*vol_v(ijk)

         ENDDO
!!$omp end parallel do

      ELSE

!$omp parallel do default(none) &
!$omp private(IJK, J, IJKN, tmp_a, tmp_B) &
!$omp shared(IJKSTART3, IJKEND3, J_OF, F_GDS, DRAG_AM, DRAG_BM, A_M, B_M, VOL_V)
         DO ijk = ijkstart3, ijkend3
            IF(fluid_at(ijk)) THEN
               j = j_of(ijk)
               ijkn = north_of(ijk)

               tmp_a = avg_y(f_gds(ijk), f_gds(ijkn), j)
               tmp_b = avg_y(drag_bm(ijk,2), drag_bm(ijkn,2), j)

               a_m(ijk,0,0) = a_m(ijk,0,0) - vol_v(ijk) * tmp_a
               b_m(ijk,0) = b_m(ijk,0) - vol_v(ijk) * tmp_b
            ENDIF
         ENDDO
!$omp end parallel do
      ENDIF

      END SUBROUTINE gas_drag_v

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: GAS_DRAG_W                                              !
!  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
!                                                                      !
!  Purpose: Account for the equal and opposite drag force on the gas   !
!           phase due to particles by introducing the drag as a        !
!           source term.  Face centered.                               !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE gas_drag_w(A_M, B_M, IER)

! Global Variables:
!---------------------------------------------------------------------//
! Runtime Flag: Interpolate DES values
      use particle_filter, only: des_interp_scheme_enum, des_interp_garg
! Flag: Gas sees the effect of particles in gas/solids flows.
      use discretelement, only: des_oneway_coupled
! Coefficient at cell corners added to the gas momentum A matix.
      use discretelement, only: drag_am, f_gds
! Coefficient at cell corners added to gas momentum B vector.
      use discretelement, only: drag_bm
! Volume of Z-momentum cell
      use geometry, only: vol_w

! Global Parameters:
!---------------------------------------------------------------------//
! Size of IJK arrays and solids phase..
      use param, only: dimension_3, dimension_m
! Fluid grid loop bounds.
      use compar, only: ijkstart3, ijkend3
! The I, J, and K values that comprise an IJK
      use indices, only: i_of, j_of, k_of
! Flag: Fluid exists at indexed cell
      use functions, only: fluid_at
! IJK of cell to top.
      use functions, only: top_of
! IJK function for I,J,K that includes mapped indices.
      use compar, only: funijk_map_c
! Function for averaging to a scalar cell's north face.
      use fun_avg, only: avg_z
! Domain index bounds.
      use compar, only: istart2, jstart2, kstart2
      use compar, only: iend2, jend2, kend2

      IMPLICIT NONE

! Dummy Arguments:
!---------------------------------------------------------------------//
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_M(dimension_3,-3:3,0:dimension_m)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_M(dimension_3, 0:dimension_m)
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local Variables:
!---------------------------------------------------------------------//
! Grid cell indices
      INTEGER :: I, J, K, IJK, IMJK, IJMK, IMJMK, IJKT
! temporary variables for matrix A_M and vector B_M
      DOUBLE PRECISION tmp_A, tmp_B
! Averaging factor
! (=0.25 in 3D and =0.5 in 2D)
      DOUBLE PRECISION :: AVG_FACTOR
!......................................................................!

! Initialize error flag.
      ier = 0

! Skip this routine if the gas/solids are only one-way coupled.
      IF(des_oneway_coupled) RETURN

      IF(des_interp_scheme_enum == des_interp_garg)THEN

         avg_factor = 0.25d0

!!$omp parallel do schedule (guided,50) default(none)                   &
!!$omp shared(IJKSTART3, IJKEND3, FUNIJK_MAP_C, I_OF, J_OF,             &
!!$omp    K_OF, AVG_FACTOR, DRAG_AM, DRAG_BM, A_M, B_M, VOL_W)          &
!!$omp private(IJK, I, J, K, IMJK, IJMK, IMJMK, tmp_A, tmp_B)
         DO ijk = ijkstart3, ijkend3
            IF(.NOT.fluid_at(ijk)) cycle

            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)

            IF (i.LT.istart2 .OR. i.GT.iend2) cycle
            IF (j.LT.jstart2 .OR. j.GT.jend2) cycle
            IF (k.LT.kstart2 .OR. k.GT.kend2) cycle

            imjk = funijk_map_c(i-1,j,k)
            ijmk = funijk_map_c(i,j-1,k)
            imjmk = funijk_map_c(i-1,j-1,k)

            tmp_a = -avg_factor*(drag_am(ijk) + drag_am(imjk) +        &
               drag_am(ijmk) + drag_am(imjmk))

            tmp_b = -avg_factor*(drag_bm(ijk,3) + drag_bm(imjk,3) +    &
               drag_bm(ijmk,3) + drag_bm(imjmk,3))

            a_m(ijk,0,0) = a_m(ijk,0,0) + tmp_a*vol_w(ijk)
            b_m(ijk,0) = b_m(ijk,0) + tmp_b*vol_w(ijk)

         ENDDO
!!$omp end parallel do

      ELSE

!$omp parallel do default(none) &
!$omp private(IJK, K, IJKT, tmp_A, tmp_B) &
!$omp shared(IJKSTART3,IJKEND3,K_OF, F_GDS, DRAG_BM, A_M, B_M, VOL_W)
         DO ijk = ijkstart3, ijkend3
            IF(fluid_at(ijk)) THEN
               k = k_of(ijk)
               ijkt = top_of(ijk)

               tmp_a = avg_z(f_gds(ijk), f_gds(ijkt), k)
               tmp_b = avg_z(drag_bm(ijk,3), drag_bm(ijkt,3), k)

               a_m(ijk,0,0) = a_m(ijk,0,0) - vol_w(ijk) * tmp_a
               b_m(ijk,0) = b_m(ijk,0) - vol_w(ijk) * tmp_b
            ENDIF
         ENDDO
!$omp end parallel do

      ENDIF

      RETURN
      END SUBROUTINE gas_drag_w