dif_phi_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_des(M, DIF, A_M, B_M)

!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE toleranc
      USE run
      USE geometry
      USE compar
      USE sendrecv
      USE indices
      USE fun_avg
      USE functions
      USE cutcell

      IMPLICIT NONE

! 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)

      DOUBLE PRECISION, INTENT(IN) :: DIF(dimension_3)

! Fluid Cell indices
      INTEGER :: I, J, K, IJK
      INTEGER :: IMJK, IM, IJKW, IPJK, IJKE
      INTEGER :: IJMK, JM, IJKS, IJPK, IJKN
      INTEGER :: IJKM, KM, IJKB, IJKP, IJKT
!
! Diffusion parameter
      DOUBLE PRECISION :: D_f


!  Calculate convection-diffusion fluxes through each of the faces
      DO ijk = ijkstart3, ijkend3

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

         IF(.NOT.fluid_at(ijk)) cycle

         ipjk = ip_of(ijk)
         ijpk = jp_of(ijk)

         ijke = east_of(ijk)
         ijkn = north_of(ijk)

! East face (i+1/2, j, k)
         d_f = avg_x_h(dif(ijk),dif(ijke),i)*odx_e(i)*ayz(ijk)
         IF(cut_treatment_at(ijk)) THEN
            IF(cut_cell_at(ijk).AND.(.NOT.fluid_at(ipjk))) THEN
               d_f = avg_x_h(dif(ijk),dif(ijke),i)*odx_e(i)*dy(j)*dz(k)
            ENDIF
         ENDIF

         a_m(ijk,east,m) = d_f
         a_m(ipjk,west,m) = d_f

! West face (i-1/2, j, k)
         imjk = im_of(ijk)
         IF (.NOT.fluid_at(imjk)) THEN
            im = im1(i)
            ijkw = west_of(ijk)
            d_f = avg_x_h(dif(ijkw),dif(ijk),im)*odx_e(im)*ayz(imjk)
            IF(cut_treatment_at(ijk)) THEN
               IF(cut_cell_at(ijk).AND.(.NOT.fluid_at(imjk))) THEN
                  d_f = avg_x_h(dif(ijkw),dif(ijk),im)*odx_e(im)*dy(j)*dz(k)
               ENDIF
            ENDIF
            a_m(ijk,west,m) = d_f
         ENDIF


! North face (i, j+1/2, k)
         d_f = avg_y_h(dif(ijk),dif(ijkn),j)*ody_n(j)*axz(ijk)
         IF(cut_treatment_at(ijk)) THEN
            IF(cut_cell_at(ijk).AND.(.NOT.fluid_at(ijpk))) THEN
               d_f = avg_y_h(dif(ijk),dif(ijkn),j)*ody_n(j)*dx(i)*dz(k)
            ENDIF
         ENDIF
         a_m(ijk,north,m) = d_f
         a_m(ijpk,south,m) = d_f

! South face (i, j-1/2, k)
         ijmk = jm_of(ijk)
         IF (.NOT.fluid_at(ijmk)) THEN
            jm = jm1(j)
            ijks = south_of(ijk)
            d_f = avg_y_h(dif(ijks),dif(ijk),jm)*ody_n(jm)*axz(ijmk)
            IF(cut_treatment_at(ijk)) THEN
               IF(cut_cell_at(ijk).AND.(.NOT.fluid_at(ijmk))) THEN
                  d_f = avg_y_h(dif(ijks),dif(ijk),jm)*ody_n(jm)*dx(i)*dz(k)
               ENDIF
            ENDIF
            a_m(ijk,south,m) = d_f
         ENDIF


! Top face (i, j, k+1/2)
         IF (do_k) THEN
            ijkp = kp_of(ijk)
            ijkt = top_of(ijk)
            d_f = avg_z_h(dif(ijk),dif(ijkt),k)*ox(i)*odz_t(k)*axy(ijk)
            IF(cut_treatment_at(ijk)) THEN
               IF(cut_cell_at(ijk).AND.(.NOT.fluid_at(ijkp))) THEN
                  d_f = avg_z_h(dif(ijk),dif(ijkt),k)*ox(i)*odz_t(k)*dx(i)*dy(j)
               ENDIF
            ENDIF
            a_m(ijk,top,m) = d_f
            a_m(ijkp,bottom,m) = d_f


! Bottom face (i, j, k-1/2)
            ijkm = km_of(ijk)
            IF (.NOT.fluid_at(ijkm)) THEN
               km = km1(k)
               ijkb = bottom_of(ijk)
               d_f = avg_z_h(dif(ijkb),dif(ijk),km)*ox(i)*odz_t(km)*axy(ijkm)
               IF(cut_treatment_at(ijk)) THEN
                  IF(cut_cell_at(ijk).AND.(.NOT.fluid_at(ijkm))) THEN
                     d_f = avg_z_h(dif(ijkb),dif(ijk),km)*ox(i)*odz_t(km)*dx(i)*dy(j)
                  ENDIF
               ENDIF
               a_m(ijk,bottom,m) = d_f
            ENDIF
         ENDIF
      END DO

      CALL dif_phi_is(dif, a_m, m)


      RETURN
      END SUBROUTINE dif_phi_des