dif_phi_bc_des.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DIFFUSE_MEAN_FIELDS                                     !
!  Author: J.Musser                                   Date: 11-NOV-14  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE dif_phi_bc_des(PHI, M, A_M, B_M)

      USE param
      USE param1
      USE geometry
      USE indices
      USE bc
      USE compar
      USE cutcell, only : cartesian_grid, cg_safe_mode
      USE fun_avg
      USE functions

      IMPLICIT NONE

!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------


      DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)
! Phase index
      INTEGER, INTENT(IN) :: M
! 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)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Boundary condition index
      INTEGER :: L
! Indices
      INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK, &
                 IM, JM, KM
!-----------------------------------------------

! Set up the default walls (i.e., bc_type='dummy' or undefined/default
! boundaries) as non-conducting...
! ---------------------------------------------------------------->>>
! when setting up default walls do not use cutcells to avoid conflict

      IF(.NOT.cartesian_grid) THEN

! west yz plane
         i1 = imin2
         DO k1 = kmin3, kmax3
         DO j1 = jmin3, jmax3
            IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
            IF (dead_cell_at(i1,j1,k1)) cycle  ! skip dead cells
            ijk = funijk(i1,j1,k1)
            IF (default_wall_at(ijk)) THEN
! Cutting the neighbor link between fluid cell and wall cell
               a_m(ip_of(ijk),west,m) = zero
! Setting the wall value equal to the adjacent fluid cell value
               a_m(ijk,:,m) = zero
               a_m(ijk,east,m) = one
               a_m(ijk,0,m) = -one
               b_m(ijk,m) = zero
            ENDIF
         ENDDO
         ENDDO

! east yz plane
         i1 = imax2
         DO k1 = kmin3, kmax3
         DO j1 = jmin3, jmax3
            IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
            IF (dead_cell_at(i1,j1,k1)) cycle  ! skip dead cells
            ijk = funijk(i1,j1,k1)
            IF (default_wall_at(ijk)) THEN
! Cutting the neighbor link between fluid cell and wall cell
               a_m(im_of(ijk),east,m) = zero
! Setting the wall value equal to the adjacent fluid cell value
               a_m(ijk,:,m) = zero
               a_m(ijk,west,m) = one
               a_m(ijk,0,m) = -one
               b_m(ijk,m) = zero
            ENDIF
         ENDDO
         ENDDO

! south xz plane
         j1 = 1
         DO k1 = kmin3, kmax3
         DO i1 = imin3, imax3
            IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
            IF (dead_cell_at(i1,j1,k1)) cycle  ! skip dead cells
            ijk = funijk(i1,j1,k1)
            IF (default_wall_at(ijk)) THEN
! Cutting the neighbor link between fluid cell and wall cell
               a_m(jp_of(ijk),south,m) = zero
! Setting the wall value equal to the adjacent fluid cell value
               a_m(ijk,:,m) = zero
               a_m(ijk,north,m) = one
               a_m(ijk,0,m) = -one
               b_m(ijk,m) = zero
            ENDIF
         ENDDO
         ENDDO

! north xz plane
         j1 = jmax2
         DO k1 = kmin3, kmax3
         DO i1 = imin3, imax3
            IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
            IF (dead_cell_at(i1,j1,k1)) cycle  ! skip dead cells
            ijk = funijk(i1,j1,k1)
            IF (default_wall_at(ijk)) THEN
! Cutting the neighbor link between fluid cell and wall cell
               a_m(jm_of(ijk),north,m) = zero
! Setting the wall value equal to the adjacent fluid cell value
               a_m(ijk,:,m) = zero
               a_m(ijk,south,m) = one
               a_m(ijk,0,m) = -one
               b_m(ijk,m) = zero
            ENDIF
         ENDDO
         ENDDO

         IF (do_k) THEN
! bottom xy plane
            k1 = 1
            DO j1 = jmin3, jmax3
            DO i1 = imin3, imax3
               IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
               IF (dead_cell_at(i1,j1,k1)) cycle  ! skip dead cells
               ijk = funijk(i1,j1,k1)

               IF (default_wall_at(ijk)) THEN
! Cutting the neighbor link between fluid cell and wall cell
                  a_m(kp_of(ijk),bottom,m) = zero
! Setting the wall value equal to the adjacent fluid cell value (set
! the boundary cell value equal to adjacent fluid cell value)
                  a_m(ijk,:,m) = zero
                  a_m(ijk,top,m) = one
                  a_m(ijk,0,m) = -one
                  b_m(ijk,m) = zero
               ENDIF
            ENDDO
            ENDDO

! top xy plane
            k1 = kmax2
            DO j1 = jmin3, jmax3
            DO i1 = imin3, imax3
               IF (.NOT.is_on_mype_plus2layers(i1,j1,k1)) cycle
               IF (dead_cell_at(i1,j1,k1)) cycle  ! skip dead cells
               ijk = funijk(i1,j1,k1)
               IF (default_wall_at(ijk)) THEN
! Cutting the neighbor link between fluid cell and wall cell
                  a_m(km_of(ijk),top,m) = zero
! Setting the wall value equal to the adjacent fluid cell value
                  a_m(ijk,:,m) = zero
                  a_m(ijk,bottom,m) = one
                  a_m(ijk,0,m) = -one
                  b_m(ijk,m) = zero
               ENDIF
            ENDDO
            ENDDO
         ENDIF

      ENDIF

! Set user defined wall boundary conditions
! ---------------------------------------------------------------->>>
      DO l = 1, dimension_bc
         IF (bc_defined(l)) THEN

            i1 = bc_i_w(l)
            i2 = bc_i_e(l)
            j1 = bc_j_s(l)
            j2 = bc_j_n(l)
            k1 = bc_k_b(l)
            k2 = bc_k_t(l)

            DO k = k1, k2
            DO j = j1, j2
            DO i = i1, i2

               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)


! Set the boundary cell equal to the known value in that
! cell
               a_m(ijk,:,m) = zero
               b_m(ijk,m) = zero
! second modify the matrix equation according to the user specified
! boundary condition
               IF (fluid_at(east_of(ijk))) THEN
                  a_m(ijk,0,m) = -odx_e(i)
                  a_m(ijk,east,m) =  odx_e(i)

               ELSEIF (fluid_at(west_of(ijk))) THEN
                  im = im1(i)
                  a_m(ijk,west,m) =  odx_e(im)
                  a_m(ijk,0,m) = -odx_e(im)

               ELSEIF (fluid_at(north_of(ijk))) THEN
                  a_m(ijk,0,m) = -ody_n(j)
                  a_m(ijk,north,m) =  ody_n(j)

               ELSEIF (fluid_at(south_of(ijk))) THEN
                  jm = jm1(j)
                  a_m(ijk,south,m) =  ody_n(jm)
                  a_m(ijk,0,m) = -ody_n(jm)

               ELSEIF (fluid_at(top_of(ijk))) THEN
                  a_m(ijk,0,m) = -ox(i)*odz_t(k)
                  a_m(ijk,top,m) =  ox(i)*odz_t(k)

               ELSEIF (fluid_at(bottom_of(ijk))) THEN
                  km = km1(k)
                  a_m(ijk,bottom,m) =  ox(i)*odz_t(km)
                  a_m(ijk,0,m) = -ox(i)*odz_t(km)

               ENDIF
            ENDDO
            ENDDO
            ENDDO
         ENDIF   ! end if (bc_defined)
      ENDDO   ! end L do loop over dimension_bc

! modifications for cartesian grid implementation
      IF(cartesian_grid .AND. .NOT.(cg_safe_mode(1)==1)) &
         CALL dif_phi_bc_des_cg(phi, 0, a_m, b_m)

      RETURN
      END SUBROUTINE dif_phi_bc_des


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: BC_PHI_CG                                               C
!  Purpose: Modify boundary conditions for cartesian grid cut-cell     C
!           implementation                                             C
!                                                                      C
!  Author: Jeff Dietiker                                               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE dif_phi_bc_des_cg(PHI, M, A_M, B_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE indices
      USE bc
      USE compar
      USE cutcell
      USE fun_avg
      USE functions

      IMPLICIT NONE

!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! The field variable being solved for:
!     e.g., T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
!     e_Turb_G
      DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)

! Phase index
      INTEGER, INTENT(IN) :: M
! 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)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IJK, IM, JM, KM

      LOGICAL :: ALONG_GLOBAL_GHOST_LAYER

!-----------------------------------------------

      DO ijk = ijkstart3, ijkend3

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

         along_global_ghost_layer = (i < imin1) .OR. ( i>imax1 )       &
            .OR. (j < jmin1) .OR. (j > jmax1)

         IF(do_k) along_global_ghost_layer = along_global_ghost_layer  &
            .OR. (k < kmin1) .OR. (k > kmax1)

! Setting the value in the boundary cell equal to what is known
         IF(blocked_cell_at(ijk)) THEN
            a_m(ijk,:,m) = zero
            a_m(ijk,0,m) = -one
            b_m(ijk,m) = -phi(ijk)
         ENDIF

         IF(blocked_cell_at(ijk).OR.along_global_ghost_layer) THEN


            IF (cut_cell_at(ip_of(ijk))) THEN
               IF( bc_id(ip_of(ijk)) > 0 ) THEN
                  a_m(ijk,0,m) = -odx_e(i)
                  a_m(ijk,east,m) =  odx_e(i)
                  b_m(ijk,m) =  zero
               ENDIF

            ELSEIF (cut_cell_at(im_of(ijk))) THEN
               IF(bc_id(im_of(ijk)) > 0 ) THEN
                  im = im1(i)
                  a_m(ijk,west,m) =  odx_e(im)
                  a_m(ijk,0,m) = -odx_e(im)
                  b_m(ijk,m) = zero
               ENDIF

            ELSEIF (cut_cell_at(jp_of(ijk))) THEN
               IF(bc_id(jp_of(ijk)) > 0 ) THEN
                  a_m(ijk,0,m) = -ody_n(j)
                  a_m(ijk,north,m) =  ody_n(j)
                  b_m(ijk,m) = zero
               ENDIF

            ELSEIF (cut_cell_at(jm_of(ijk))) THEN
               IF(bc_id(jm_of(ijk)) > 0 ) THEN
                  jm = jm1(j)
                  a_m(ijk,south,m) =  ody_n(jm)
                  a_m(ijk,0,m) = -ody_n(jm)
                  b_m(ijk,m) = zero
               ENDIF

            ELSEIF (cut_cell_at(kp_of(ijk))) THEN
               IF(bc_id(kp_of(ijk)) > 0 ) THEN
                  a_m(ijk,0,m) = -ox(i)*odz_t(k)
                  a_m(ijk,top,m) =  ox(i)*odz_t(k)
                  b_m(ijk,m)  = zero
               ENDIF

            ELSEIF (cut_cell_at(km_of(ijk))) THEN
               IF(bc_id(km_of(ijk)) > 0 ) THEN
                  km = km1(k)
                  a_m(ijk,bottom,m) =  ox(i)*odz_t(km)
                  a_m(ijk,0,m) = -ox(i)*odz_t(km)
                  b_m(ijk,m) = zero
               ENDIF
            ENDIF
         ENDIF
      ENDDO

      RETURN

      END SUBROUTINE dif_phi_bc_des_cg