set_geometry.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: SET_GEOMETRY                                           !
!  Author: M. Syamlal                                 Date: 21-JAN-92  !
!                                                                      !
!  Purpose: Calculate X, X_E,  oX, oX_E                                !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE set_geometry

! Global Variables:
!---------------------------------------------------------------------//
! Domain decomposition and dimensions
      use geometry, only: dx, xlength, odx, odx_e
      use geometry, only: dy, odz, odz_t
      use geometry, only: dz, zlength, ody, ody_n
      use geometry, only: x, x_e, ox, ox_e, xmin
      use geometry, only: z, z_t
! Domain indices.
      use geometry, only: do_i, imin1, imax1, imax2, imax3, imin3, imax
      use geometry, only: do_j, jmin1, jmax1, jmax2, jmax3, jmin3
      use geometry, only: do_k, kmin1, kmax1, kmax2, kmax3, kmin3
! Averaging factors.
      use geometry, only: fx_e, fx_e_bar, fx, fx_bar
      use geometry, only: fy_n, fy_n_bar
      use geometry, only: fz_t, fz_t_bar
! Cyclic domain flags.
      use geometry, only: cyclic
      use geometry, only: cyclic_x, cyclic_x_pd, cyclic_x_mf
      use geometry, only: cyclic_y, cyclic_y_pd, cyclic_y_mf
      use geometry, only: cyclic_z, cyclic_z_pd, cyclic_z_mf
! Flag for cylindrical coordinates.
      use geometry, only: cylindrical
! For cylindrical_2d simulations
      use geometry, only: cylindrical_2d, i_cyl_num, i_cyl_transition
      use geometry, only: cyl_x, cyl_x_e
! MPI-Domain decompoint and rank flags.
      use compar, only: nodesi, nodesj, nodesk
! Flag for specificed constant mass flux.
      use bc, only: flux_g

! Global Parameters:
!---------------------------------------------------------------------//
      use param1, only: zero, half, one, undefined

! Module procedures
!---------------------------------------------------------------------//
      use mpi_utility, only: global_all_sum
      use error_manager
      use toleranc

      IMPLICIT NONE

! Local Variables:
!---------------------------------------------------------------------//
! Generic loop indices
      INTEGER :: I, J, K
! X-direction dimension of U-momentum cell
      DOUBLE PRECISION :: DX_E
! Y-direction dimension of V-momentum cell
      DOUBLE PRECISION :: DY_N
! Z-direction dimension of W-momentum cell
      DOUBLE PRECISION :: DZ_T
! Local variables for cylindrical_2d simulation
      integer i_cyl_min, i_cyl_max
      double precision l_ver, l_ab, rrr, ddy
!......................................................................!

! Initialize the error manager.
      CALL init_err_msg("SET_GEOMETRY")

! Allocate geometry arrays.
      CALL allocate_arrays_geometry

!  Determine the cyclic direction with a specified mass flux
      cyclic_x_mf = (flux_g /= undefined .AND. cyclic_x_pd)
      cyclic_y_mf = (flux_g /= undefined .AND. cyclic_y_pd)
      cyclic_z_mf = (flux_g /= undefined .AND. cyclic_z_pd)

! Force the cyclic flag if cyclic with pressure drop.
      IF (cyclic_x_pd) cyclic_x = .true.
      IF (cyclic_y_pd) cyclic_y = .true.
      IF (cyclic_z_pd) cyclic_z = .true.
      cyclic = cyclic_x .OR. cyclic_y .OR. cyclic_z

      IF(cylindrical .AND. compare(zlength,8.d0*atan(one)) .AND. do_k) &
         cyclic_z = .true.

      IF(cyclic_x) THEN
         dx(1) = dx(imax1)
         dx(imax2) = dx(imin1)
         IF(nodesi.NE.1) THEN
            dx(imin3) = dx(imax1-1)
            dx(imax3) = dx(imin1+1)
         ENDIF
      ENDIF

      IF(cyclic_y) THEN
         dy(1) = dy(jmax1)
         dy(jmax2) = dy(jmin1)
         IF(nodesj.NE.1) THEN
            dy(jmin3) = dy(jmax1-1)
            dy(jmax3) = dy(jmin1+1)
         ENDIF
      ENDIF

      IF (cyclic_z) THEN
         dz(1) = dz(kmax1)
         dz(kmax2) = dz(kmin1)
         IF(nodesk.NE.1) THEN
            dz(kmin3) = dz(kmax1-1)
            dz(kmax3) = dz(kmin1+1)
         ENDIF
      ENDIF


! Initialize the X-axis variables for CYLINDRICAL coordinates.
      IF(cylindrical) THEN
         x(imin3:imax3)    = undefined
         x_e(imin3:imax3)  = undefined
         ox(imin3:imax3)   = undefined
         ox_e(imin3:imax3) = undefined
         odx(imin3:imax3)  = undefined

         IF(xmin == zero) THEN
            odx(1) = one/dx(1)
            ox(1) = undefined
            ox_e(1) = undefined
            IF (do_i) THEN
               x(1) = -half*dx(1)
               x_e(1) = 0.0
            ELSE
               x(1) = half*dx(1)
               x_e(1) = dx(1)
            ENDIF
         ELSE
            IF (do_i) THEN
               x_e(1) = xmin
               x(1) = xmin - half*dx(1)
            ELSE
               x_e(1) = xmin + dx(1)
               x(1) = xmin + half*dx(1)
            ENDIF
            ox(1) = one/x(1)
            ox_e(1) = one/x_e(1)
            odx(1) = one/dx(1)
         ENDIF

         IF (do_i) THEN
            DO i = imin1, imax2
               x(i) = x(i-1) + (dx(i-1)+dx(i))/2.
               x_e(i) = x_e(i-1) + dx(i)
               ox(i) = one/x(i)
               ox_e(i) = one/x_e(i)
               odx(i) = one/dx(i)
            END DO
         ENDIF

! Initialize the X-axis variables for CARTESIAN coordinates.
      ELSE

         x(imin3:imax3) = one
         x_e(imin3:imax3) = one
         ox(imin3:imax3) = one
         ox_e(imin3:imax3) = one
         odx(imin3:imax3) = one/dx(imin3:imax3)

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++<<
!     Implementation of the 2.5D model by Li et al. CES 123 (2015) 236-246.
!     A computational domain of two wedges connected by a thin plate is used.
!     The model is invoked by setting CYLINDRICAL_2D to .TRUE.
!     The half width of the plate is determined by I_CYL_NUM
!     Please refer to the paper for details of the model
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++<<
         IF(cylindrical_2d)THEN

            IF (xmin == zero) THEN
               IF (do_i) THEN
                  cyl_x(1) = -half*dx(1) -half*xlength
                  cyl_x_e(1) = 0.d0 -half*xlength
               ELSE
                  cyl_x(1) = half*dx(1) -half*xlength
                  cyl_x_e(1) = dx(1) -half*xlength
               ENDIF
            ELSE
               IF (do_i) THEN
                  cyl_x_e(1) = xmin -half*xlength
                  cyl_x(1) = xmin - half*dx(1) -half*xlength
               ELSE
                  cyl_x_e(1) = xmin + dx(1) -half*xlength
                  cyl_x(1) = xmin + half*dx(1) -half*xlength
               ENDIF
            ENDIF

            IF (do_i) THEN
               DO i = imin1, imax2
                  cyl_x(i) = cyl_x(i-1) + (dx(i-1)+dx(i))/2.
                  cyl_x_e(i) = cyl_x_e(i-1) + dx(i)
               END DO

               DO i = imin3, imax3
                  cyl_x(i) = cyl_x(i) *2.d0*sin(dz(1)/2.d0)
                  cyl_x_e(i) = cyl_x_e(i) *2.d0*sin(dz(1)/2.d0)
                  cyl_x(i) = abs(cyl_x(i))
                  cyl_x_e(i) = abs(cyl_x_e(i))
               END DO

               if(mod(imax,2).eq.1)then     ! odd
                  i_cyl_min = (imax+1)/2 + 1 - i_cyl_num
                  i_cyl_max = (imax+1)/2 + 1 + i_cyl_num

                  do i=i_cyl_min, i_cyl_max
                     cyl_x(i)=dx(i)*(i_cyl_num + half) *2.d0*sin(dz(1)/2.d0)
                  enddo
                  do i=i_cyl_min ,i_cyl_max -1
                     cyl_x_e(i)=dx(i)*(i_cyl_num + half) *2.d0*sin(dz(1)/2.d0)
                  enddo
               else                    ! even
                  i_cyl_min = (imax)/2 + 2 - i_cyl_num
                  i_cyl_max = (imax)/2 + 1 + i_cyl_num

                  do i=i_cyl_min, i_cyl_max
                     cyl_x(i)=dx(i)* i_cyl_num *2.d0*sin(dz(1)/2.d0)
                  enddo
                  do i=i_cyl_min, i_cyl_max -1
                     cyl_x_e(i)=dx(i)* i_cyl_num *2.d0*sin(dz(1)/2.d0)
                  enddo
               endif

!  To smooth the transition from cylindrical to 2D domain
!  Only one or two cells transition are implemented.

               if(i_cyl_transition .eq. 1)then
                  l_ver=dx(i_cyl_min)*tan(dz(1)/2.0)
                  l_ab = sqrt(l_ver*l_ver + (dx(i_cyl_min-1)+dx(i_cyl_min))**2)
                  rrr = half*l_ab/sin(dz(1)/4.0)
                  ddy = rrr - sqrt(rrr**2 - dx(i_cyl_min)**2)

                  cyl_x_e(i_cyl_min - 1) = cyl_x_e(i_cyl_min) + 2*ddy
                  cyl_x_e(i_cyl_max ) = cyl_x_e(i_cyl_max - 1) + 2*ddy
               elseif(i_cyl_transition .eq. 2) then
                  l_ver=2.d0*dx(i_cyl_min)*tan(dz(1)/2.0)
                  l_ab = sqrt(l_ver*l_ver + (4.d0*dx(i_cyl_min))**2)
                  rrr = half*l_ab/sin(dz(1)/4.0)
                  ddy = rrr - sqrt(rrr**2 - dx(i_cyl_min)**2)

                  cyl_x_e(i_cyl_min ) = cyl_x_e(i_cyl_min + 1) + 2*ddy
                  cyl_x_e(i_cyl_max - 1) = cyl_x_e(i_cyl_max - 2) + 2*ddy

                  ddy = rrr - sqrt(rrr**2 - (2.d0*dx(i_cyl_min))**2)
                  cyl_x_e(i_cyl_min - 1 ) = cyl_x_e(i_cyl_min + 1) + 2*ddy
                  cyl_x_e(i_cyl_max ) = cyl_x_e(i_cyl_max - 2) + 2*ddy

                  ddy = rrr - sqrt(rrr**2 - (3.d0*dx(i_cyl_min))**2)
                  cyl_x_e(i_cyl_min - 2 ) = cyl_x_e(i_cyl_min + 1) + 2*ddy
                  cyl_x_e(i_cyl_max + 1) = cyl_x_e(i_cyl_max - 2) + 2*ddy
               else ! no smooth

               endif

               do i=i_cyl_min - i_cyl_transition, i_cyl_max + i_cyl_transition
                  cyl_x(i) = half * (cyl_x_e(i-1) + cyl_x_e(i))
               enddo

!  To make sure all locations are positive
               do i= imin3,imax3
                  cyl_x(i)=abs(cyl_x(i))
                  cyl_x_e(i)=abs(cyl_x_e(i))
               enddo
            ENDIF  !IF (DO_I)
         ENDIF  !IF(CYLINDRICAL_2D)THEN

      ENDIF


! Initialize the Y-Axis variables.
      ody(jmin3:jmax3) = one/dy(jmin3:jmax3)


! Initialize the Z-Axis variables.
      DO k = 1, kmax3
         IF (k == 1) THEN
            z(k) = zero - half*dz(k)
            z_t(k) = zero
            IF(nodesk.NE.1) THEN
               z(k-1) =z_t(k) - half*dz(k-1)
               z_t(k-1) = z_t(k) - dz(k-1)
            ENDIF
              ELSE
            z(k) = z_t(k-1) + half*dz(k)
            z_t(k) = z_t(k-1) + dz(k)
         ENDIF
         odz(k) = one/dz(k)

         IF(nodesk.NE.1 .AND. k==1) odz(k-1) = one/dz(k-1)
      END DO


      dx_e = half*(dx(1)+dx(imin1))
      dy_n = half*(dy(1)+dy(jmin1))
      dz_t = half*(dz(1)+dz(kmin1))

      odx_e(1) = one/dx_e
      ody_n(1) = one/dy_n
      odz_t(1) = one/dz_t

      fx(1) = half
      fx_bar(1) = half
      fx_e(1) = half
      fx_e_bar(1) = half
      fy_n(1) = half
      fy_n_bar(1) = half
      fz_t(1) = half
      fz_t_bar(1) = half

      IF(nodesi.NE.1) odx_e(imin3) = one/dx_e
      IF(nodesj.NE.1) ody_n(jmin3) = one/dy_n
      IF(nodesk.NE.1) odz_t(kmin3) = one/dz_t

      IF(nodesi.NE.1) THEN
         fx(imin3) = half
         fx_bar(imin3) = half
         fx_e(imin3) = half
         fx_e_bar(imin3) = half
      ENDIF

      IF(nodesj.NE.1) THEN
         fy_n(jmin3) = half
         fy_n_bar(jmin3) = half
      ENDIF

      IF(nodesk.NE.1) THEN
         fz_t(kmin3) = half
         fz_t_bar(kmin3) = half
      ENDIF


! Look at 2 through IMAX1 U-momentum cells
      IF (do_i) THEN
         DO i = imin1, imax1
            dx_e = half*(dx(i+1)+dx(i))
            odx_e(i) = one/dx_e
            fx(i) = half
            fx_bar(i) = one - fx(i)
            fx_e(i) = dx(i+1)/(dx(i+1)+dx(i))
            fx_e_bar(i) = one - fx_e(i)
         END DO
      ENDIF

! Look at 2 through JMAX1 V-momentum cells
      IF (do_j) THEN
         DO j = jmin1, jmax1
            dy_n = half*(dy(j+1)+dy(j))
            ody_n(j) = one/dy_n
            fy_n(j) = dy(j+1)/(dy(j+1)+dy(j))
            fy_n_bar(j) = one - fy_n(j)
         END DO
      ENDIF

! Look at 2 through KMAX1 W-momentum cells
      IF (do_k) THEN
         DO k = kmin1, kmax1
            dz_t = half*(dz(k+1)+dz(k))
            odz_t(k) = one/dz_t
            fz_t(k) = dz(k+1)/(dz(k+1)+dz(k))
            fz_t_bar(k) = one - fz_t(k)
         END DO
      ENDIF

! Look at last U-, V-, and W-momentum cells
      dx_e = dx(imax2)
      dy_n = dy(jmax2)
      dz_t = dz(kmax2)
      odx_e(imax2) = one/dx_e
      ody_n(jmax2) = one/dy_n
      odz_t(kmax2) = one/dz_t
      fx(imax2) = half
      fx_bar(imax2) = half
      fx_e(imax2) = half
      fx_e_bar(imax2) = half

      fy_n(jmax2) = half
      fy_n_bar(jmax2) = half

      fz_t(kmax2) = half
      fz_t_bar(kmax2) = half
      fz_t(kmax3) = half
      fz_t_bar(kmax3) = half

      odx_e(imax3) = one/dx_e
      ody_n(jmax3) = one/dy_n
      odz_t(kmax3) = one/dz_t
      fx(imax3) = half
      fx_bar(imax3) = half
      fx_e(imax3) = half
      fx_e_bar(imax3) = half

      fy_n(jmax3) = half
      fy_n_bar(jmax3) = half


      IF(cyclic_x) THEN
         fx_e(1) = fx_e(imax1)
         fx_e_bar(1) = fx_e_bar(imax1)
          IF(nodesi.NE.1) THEN
            fx_e(imin3) = fx_e(imax1-1)
            fx_e_bar(imin3) = fx_e_bar(imax1-1)
         ENDIF
      ENDIF
      IF (cyclic_y) THEN
         fy_n(1) = fy_n(jmax1)
         fy_n_bar(1) = fy_n_bar(jmax1)
         IF(nodesj.NE.1) THEN
            fy_n(jmin3) = fy_n(jmax1-1)
            fy_n_bar(jmin3) = fy_n_bar(jmax1-1)
         ENDIF
      ENDIF
      IF (cyclic_z) THEN
         fz_t(1) = fz_t(kmax1)
         fz_t_bar(1) = fz_t_bar(kmax1)
         IF(nodesk.NE.1) THEN
            fz_t(kmin3) = fz_t(kmax1-1)
            fz_t_bar(kmin3) = fz_t_bar(kmax1-1)
         ENDIF
      ENDIF

      CALL finl_err_msg

      RETURN
      END SUBROUTINE set_geometry