set_bc1.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_BC1                                                 C
!  Purpose: Set transient flow boundary conditions                     C
!                                                                      C
!  Author: M. Syamlal                                 Date: 29-JAN-92  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc1

! Modules
!---------------------------------------------------------------------//
      use bc
      USE param, only: dimension_bc
      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! index for boundary condition
      INTEGER :: L


! Set the boundary conditions
      DO l = 1, dimension_bc
         IF (bc_defined(l)) THEN

            SELECT CASE(bc_type_enum(l))
            CASE (p_outflow)
               CALL set_outflow(l)
               CALL set_bc1_report_outflow(l)
            CASE (mass_outflow)
               CALL set_outflow(l)
               CALL set_bc1_adjust_outflow(l)
            CASE (mass_inflow)
               CALL set_bc1_jet(l)
            CASE (p_inflow)
               CALL set_outflow(l)
            CASE (outflow)
               CALL set_outflow(l)
               CALL set_bc1_report_outflow(l)
            END SELECT
         ENDIF   ! end if (bc_defined(l))
      ENDDO    ! end do loop (l=1,dimension_bc)

      RETURN
      END SUBROUTINE set_bc1



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc1_jet                                             C
!  Purpose: update transient jet conditions                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc1_jet(BCV)

! Modules
!---------------------------------------------------------------------//
      use bc, only: bc_plane
      use bc, only: bc_k_b, bc_k_t
      use bc, only: bc_j_s, bc_j_n
      use bc, only: bc_i_w, bc_i_e
      use bc, only: bc_jet_g
      use bc, only: bc_time
      use bc, only: bc_jet_gh, bc_dt_h
      use bc, only: bc_jet_gl, bc_dt_l

      use fldvar, only: u_g, v_g, w_g

      use run, only: time, dt

      use param1, only: undefined

      use functions, only: im_of, jm_of, km_of
      use functions, only: is_on_mype_plus2layers
      use functions, only: funijk
      use compar, only: dead_cell_at

      IMPLICIT NONE
! Dummy arguments
!---------------------------------------------------------------------//
! index for boundary condition
      INTEGER, INTENT(IN) :: BCV

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: I, J, K, IJK
! IJK index for setting velocity bc
      INTEGER :: IJK2
!---------------------------------------------------------------------//

      IF (time + 0.1d0*dt>=bc_time(bcv) .AND. &
          bc_jet_g(bcv)/=undefined) THEN

         IF (bc_jet_g(bcv) == bc_jet_gh(bcv)) THEN
            bc_jet_g(bcv) = bc_jet_gl(bcv)
            bc_time(bcv) = time + bc_dt_l(bcv)
         ELSEIF (bc_jet_g(bcv) == bc_jet_gl(bcv)) THEN
            bc_jet_g(bcv) = bc_jet_gh(bcv)
            bc_time(bcv) = time + bc_dt_h(bcv)
         ELSE
            bc_jet_g(bcv) = bc_jet_gh(bcv)
            bc_time(bcv) = time + bc_dt_h(bcv)
         ENDIF

         DO k = bc_k_b(bcv), bc_k_t(bcv)
         DO j = bc_j_s(bcv), bc_j_n(bcv)
         DO i = bc_i_w(bcv), bc_i_e(bcv)
            IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
            IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
            ijk = funijk(i,j,k)
! Why is the velocity of the boundary cell not always set (in case of
! w, s, or b plane then velocity of the adjacent fluid cell is set)?
! It should not really matter for MI...
            SELECT CASE (trim(bc_plane(bcv)))
            CASE ('W')
               ijk2 = im_of(ijk)
               u_g(ijk2) = bc_jet_g(bcv)
            CASE ('E')
               u_g(ijk) = bc_jet_g(bcv)
            CASE ('S')
               ijk2 = jm_of(ijk)
               v_g(ijk2) = bc_jet_g(bcv)
            CASE ('N')
               v_g(ijk) = bc_jet_g(bcv)
            CASE ('B')
               ijk2 = km_of(ijk)
               w_g(ijk2) = bc_jet_g(bcv)
            CASE ('T')
               w_g(ijk) = bc_jet_g(bcv)
            END SELECT
         ENDDO
         ENDDO
         ENDDO
      ENDIF
      RETURN
      END SUBROUTINE set_bc1_jet


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc1_report_outflow                                  C
!  Purpose: print out outflow conditions                               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc1_report_outflow(BCV)

! Modules
!---------------------------------------------------------------------//
      use bc, only: bc_dt_0, bc_time
      use bc, only: bc_out_n
      use bc, only: bc_mout_g, bc_mout_s
      use bc, only: bc_vout_g, bc_vout_s
      use funits, only: dmp_log, unit_log
      use param1, only: undefined, zero
      use physprop, only: smax
      use machine, only: start_log, end_log
      use run, only: time, dt, tstop
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! index for boundary condition
      INTEGER, INTENT(IN) :: BCV

! Local variables
!---------------------------------------------------------------------//
! Solids phase index
      INTEGER :: M

      IF (bc_dt_0(bcv) == undefined) RETURN

      CALL calc_outflow(bcv)

! Calculate and accumulate the actual mass and volume outflow
      IF (time + 0.1d0*dt>=bc_time(bcv) .OR. &
          time+0.1d0*dt>=tstop) THEN
         bc_time(bcv) = time + bc_dt_0(bcv)

! Average and print out the flow rates
         bc_mout_g(bcv) = abs(bc_mout_g(bcv))/bc_out_n(bcv)
         bc_vout_g(bcv) = abs(bc_vout_g(bcv))/bc_out_n(bcv)
         CALL start_log
         IF(dmp_log)WRITE (unit_log, 1000) bcv, time
         IF(dmp_log)WRITE (unit_log, 1100) bc_mout_g(bcv), &
            bc_vout_g(bcv)
         bc_mout_g(bcv) = zero
         bc_vout_g(bcv) = zero
         DO m = 1, smax
            bc_mout_s(bcv,m) = abs(bc_mout_s(bcv,m))/bc_out_n(bcv)
            bc_vout_s(bcv,m) = abs(bc_vout_s(bcv,m))/bc_out_n(bcv)
            IF(dmp_log)WRITE (unit_log, 1200) m, &
               bc_mout_s(bcv,m), bc_vout_s(bcv,m)
            bc_mout_s(bcv,m) = zero
            bc_vout_s(bcv,m) = zero
         ENDDO
         CALL end_log
         bc_out_n(bcv) = 0
      ENDIF

 1000 FORMAT(/,1x,'Average outflow rates at BC No. ',i2,'  At Time = ',g12.5)
 1100 FORMAT(3x,'Gas : Mass flow = ',g12.5,'     Volumetric flow = ',g12.5)
 1200 FORMAT(3x,'Solids-',i1,' : Mass flow = ',g12.5,&
         '     Volumetric flow = ',g12.5)

      RETURN
      END SUBROUTINE set_bc1_report_outflow


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc1_adjust_outflow                                  C
!  Purpose: Adjust velocities to get specified mass or volumetric      C
!  flow rate based on average outflow rate                             C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc1_adjust_outflow(BCV)

! Modules
!---------------------------------------------------------------------//
      use bc, only: bc_dt_0, bc_time
      use bc, only: bc_i_w, bc_i_e
      use bc, only: bc_j_s, bc_j_n
      use bc, only: bc_k_b, bc_k_t
      use bc, only: bc_massflow_g, bc_massflow_s
      use bc, only: bc_mout_g, bc_mout_s
      use bc, only: bc_out_n
      use bc, only: bc_plane
      use bc, only: bc_u_g, bc_v_g, bc_w_g
      use bc, only: bc_u_s, bc_v_s, bc_w_s
      use bc, only: bc_volflow_g, bc_volflow_s
      use bc, only: bc_vout_g, bc_vout_s
      use compar, only: dead_cell_at
      use fldvar, only: u_g, v_g, w_g
      use fldvar, only: u_s, v_s, w_s, rop_s
      use functions, only: funijk
      use functions, only: im_of, jm_of, km_of
      use functions, only: is_on_mype_plus2layers
      use funits, only: dmp_log, unit_log
      use machine, only: start_log, end_log
      use param1, only: undefined, zero, small_number
      use physprop, only: smax, mmax
      use run, only: kt_type_enum, ghd_2007
      use run, only: time, dt, tstop

      IMPLICIT NONE
! Dummy arguments
!---------------------------------------------------------------------//
! index for boundary condition
      INTEGER, INTENT(IN) :: BCV

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: I, J, K, IJK
! IJK index for setting velocity bc
      INTEGER :: IJK2
! Solids phase index
      INTEGER :: M
! local solids velocity for mixture (for ghd)
      DOUBLE PRECISION :: lvel_s
!---------------------------------------------------------------------//

      CALL calc_outflow(bcv)

! Calculate and accumulate the actual mass and volume outflow
      IF (time + 0.1d0*dt>=bc_time(bcv) .OR. &
          time+0.1d0*dt>=tstop) THEN
         bc_time(bcv) = time + bc_dt_0(bcv)

! Average and print out the flow rates
         bc_mout_g(bcv) = abs(bc_mout_g(bcv))/bc_out_n(bcv)
         bc_vout_g(bcv) = abs(bc_vout_g(bcv))/bc_out_n(bcv)
         CALL start_log
         IF(dmp_log)WRITE (unit_log, 1000) bcv, time
         IF(dmp_log)WRITE (unit_log, 1100) bc_mout_g(bcv), &
            bc_vout_g(bcv)
         DO m = 1, smax
            bc_mout_s(bcv,m) = abs(bc_mout_s(bcv,m))/bc_out_n(bcv)
            bc_vout_s(bcv,m) = abs(bc_vout_s(bcv,m))/bc_out_n(bcv)
            IF(dmp_log)WRITE (unit_log, 1200) m, &
               bc_mout_s(bcv,m), bc_vout_s(bcv,m)
         ENDDO
         CALL end_log
         bc_out_n(bcv) = 0

! Now that we know the mass and volume outflow update the bc velocities
! (gas phase)
         IF (bc_massflow_g(bcv) /= undefined) THEN
            IF (bc_mout_g(bcv) > small_number) THEN
               SELECT CASE (trim(bc_plane(bcv)))
               CASE ('W', 'E')
                  bc_u_g(bcv) = bc_u_g(bcv)*bc_massflow_g(bcv)/&
                     bc_mout_g(bcv)
               CASE ('S', 'N')
                  bc_v_g(bcv) = bc_v_g(bcv)*bc_massflow_g(bcv)/&
                     bc_mout_g(bcv)
               CASE ('B', 'T')
                  bc_w_g(bcv) = bc_w_g(bcv)*bc_massflow_g(bcv)/&
                     bc_mout_g(bcv)
               END SELECT
            ENDIF
         ELSEIF (bc_volflow_g(bcv) /= undefined) THEN
            IF (bc_vout_g(bcv) > small_number) THEN
               SELECT CASE (trim(bc_plane(bcv)))
               CASE ('W', 'E')
                  bc_u_g(bcv) = bc_u_g(bcv)*bc_volflow_g(bcv)/&
                     bc_vout_g(bcv)
               CASE ('S', 'N')
                  bc_v_g(bcv) = bc_v_g(bcv)*bc_volflow_g(bcv)/&
                     bc_vout_g(bcv)
               CASE ('B', 'T')
                  bc_w_g(bcv) = bc_w_g(bcv)*bc_volflow_g(bcv)/&
                     bc_vout_g(bcv)
               END SELECT
            ENDIF
         ENDIF
! zero out counter for new cycle
         bc_mout_g(bcv) = zero
         bc_vout_g(bcv) = zero

! Now that we know the mass and volume outflow update the bc velocities
! (solids phase)
         DO m = 1, smax
            IF (bc_massflow_s(bcv,m) /= undefined) THEN
               IF (bc_mout_s(bcv,m) > small_number) THEN
                  SELECT CASE (trim(bc_plane(bcv)))
                  CASE ('W', 'E')
                     bc_u_s(bcv,m) = bc_u_s(bcv,m)*&
                        bc_massflow_s(bcv,m)/bc_mout_s(bcv,m)
                  CASE ('S', 'N')
                     bc_v_s(bcv,m) = bc_v_s(bcv,m)*&
                        bc_massflow_s(bcv,m)/bc_mout_s(bcv,m)
                  CASE ('B', 'T')
                     bc_w_s(bcv,m) = bc_w_s(bcv,m)*&
                        bc_massflow_s(bcv,m)/bc_mout_s(bcv,m)
                  END SELECT
               ENDIF
            ELSEIF (bc_volflow_s(bcv,m) /= undefined) THEN
               IF (bc_vout_s(bcv,m) > small_number) THEN
                  SELECT CASE (trim(bc_plane(bcv)))
                  CASE ('W', 'E')
                     bc_u_s(bcv,m) = bc_u_s(bcv,m)*&
                        bc_volflow_s(bcv,m)/bc_vout_s(bcv,m)
                  CASE ('S', 'N')
                     bc_v_s(bcv,m) = bc_v_s(bcv,m)*&
                        bc_volflow_s(bcv,m)/bc_vout_s(bcv,m)
                  CASE ('B', 'T')
                     bc_w_s(bcv,m) = bc_w_s(bcv,m)*&
                        bc_volflow_s(bcv,m)/bc_vout_s(bcv,m)
                 END SELECT
               ENDIF
            ENDIF
! zero out counter for new cycle
            bc_mout_s(bcv,m) = zero
            bc_vout_s(bcv,m) = zero
         ENDDO

! Apply updated boundary velocities - Define the field variables at the
! boundaries according to user specifications with modifications from
! the above calculations.
! If the boundary plane is W, S, or B (i.e., the fluid cell is on the
! west, south or bottom of the boundary cell) then define the velocity
! of the adjacent fluid cell according to the boundary velocity rather
! than the velocity of the boundary cell.
! Why not set the velocity in the boundary cell itself?  Based on the
! momentum bc routine it should not really matter for MO.
         DO k = bc_k_b(bcv), bc_k_t(bcv)
         DO j = bc_j_s(bcv), bc_j_n(bcv)
         DO i = bc_i_w(bcv), bc_i_e(bcv)
            IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
            IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
            ijk = funijk(i,j,k)
            SELECT CASE (trim(bc_plane(bcv)))
            CASE ('W')
               ijk2 = im_of(ijk)
               u_g(ijk2) = bc_u_g(bcv)
            CASE ('E')
               u_g(ijk) = bc_u_g(bcv)
            CASE ('S')
               ijk2 = jm_of(ijk)
               v_g(ijk2) = bc_v_g(bcv)
            CASE ('N')
               v_g(ijk) = bc_v_g(bcv)
            CASE ('B')
               ijk2 = km_of(ijk)
               w_g(ijk2) = bc_w_g(bcv)
            CASE ('T')
               w_g(ijk) = bc_w_g(bcv)
            END SELECT
            DO m = 1, smax
               SELECT CASE (trim(bc_plane(bcv)))
               CASE ('W')
                  ijk2 = im_of(ijk)
                  u_s(ijk2,m) = bc_u_s(bcv,m)
               CASE ('E')
                  u_s(ijk,m) = bc_u_s(bcv,m)
               CASE ('S')
                  ijk2 = jm_of(ijk)
                  v_s(ijk2,m) = bc_v_s(bcv,m)
               CASE ('N')
                  v_s(ijk,m) = bc_v_s(bcv,m)
               CASE ('B')
                  ijk2 = km_of(ijk)
                  w_s(ijk2,m) = bc_w_s(bcv,m)
               CASE ('T')
                  w_s(ijk,m) = bc_w_s(bcv,m)
               END SELECT
            ENDDO

! compute mixutre velocity BC for GHD theory
            IF(kt_type_enum == ghd_2007) THEN
               lvel_s = zero

! bulk density is set by set_outflow and is set in bc according to
! neighboring fluid appropriately. so we don't need to reference
! bulk density with ijk vs ijk2 index (ijk value will = ijk2 value).
               DO m = 1, smax
                  SELECT CASE (trim(bc_plane(bcv)))
                  CASE ('W', 'E')
                     lvel_s = lvel_s + &
                        bc_u_s(bcv,m)*rop_s(ijk,m)
                  CASE ('S', 'N')
                     lvel_s = lvel_s + &
                        bc_v_s(bcv,m)*rop_s(ijk,m)
                  CASE ('B', 'T')
                     lvel_s = lvel_s + &
                        bc_w_s(bcv,m)*rop_s(ijk,m)
                  END SELECT
               ENDDO

               IF (rop_s(ijk,mmax) > 0) THEN
                  lvel_s = lvel_s /rop_s(ijk,mmax)
               ELSE
                  lvel_s = zero
               ENDIF

               SELECT CASE (trim(bc_plane(bcv)))
               CASE ('W')
                  ijk2 = im_of(ijk)
                  u_s(ijk2,mmax) =  lvel_s
               CASE ('E')
                  u_s(ijk,mmax) = lvel_s
               CASE ('S')
                  ijk2 = jm_of(ijk)
                  v_s(ijk2,mmax) = lvel_s
               CASE ('N')
                  v_s(ijk,mmax) = lvel_s
               CASE ('B')
                  ijk2 = km_of(ijk)
                  w_s(ijk2,mmax) = lvel_s
               CASE ('T')
                  w_s(ijk,mmax) = lvel_s
               END SELECT

            ENDIF   ! end if (kt_type_enum==ghd_2007)
         ENDDO
         ENDDO
         ENDDO
      ENDIF   ! if time to update outflow condition

 1000 FORMAT(/,1x,'Average outflow rates at BC No. ',i2,'  At Time = ',g12.5)
 1100 FORMAT(3x,'Gas : Mass flow = ',g12.5,'     Volumetric flow = ',g12.5)
 1200 FORMAT(3x,'Solids-',i1,' : Mass flow = ',g12.5,&
         '     Volumetric flow = ',g12.5)

      RETURN
      END SUBROUTINE set_bc1_adjust_outflow