TURB_coRansWallDistOVFlo.F90

Go to the documentation of this file.

! *********************************************************************
! * Rocstar Simulation Suite                                          *
! * Copyright@2015, Illinois Rocstar LLC. All rights reserved.        *
! *                                                                   *
! * Illinois Rocstar LLC                                              *
! * Champaign, IL                                                     *
! * www.illinoisrocstar.com                                           *
! * sales@illinoisrocstar.com                                         *
! *                                                                   *
! * License: See LICENSE file in top level of distribution package or *
! * http://opensource.org/licenses/NCSA                               *
! *********************************************************************
! *********************************************************************
! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
! * NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
! * USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
! *********************************************************************
!******************************************************************************
!
! Purpose: compute distance of each cell centers to the nearest no-slip wall
!
! Description: distance to nearest wall (walldist) using this direct method
!              is computed by collecting coordinates of no-slip wall patches 
!              in all regions and stored in a buffer. The wall distance is 
!              the minimum of the distance between the coordinates in the 
!              buffer and each cell centroid. In this way, global wall distance 
!              is obtained, assuming open view (OV) from cell centers to 
!              the nearest wall.
!
! Input: region = data of current region.
!
! Output: region%levels%turb%lens = turbulence length scale.
!
! Notes: this method is computationally intensive, not suitable for moving
!        grid. More efficient method is provided by other routines.        
!
!******************************************************************************
!
! $Id: TURB_coRansWallDistOVFlo.F90,v 1.7 2009/08/26 12:28:53 mtcampbe Exp $
!
! Copyright: (c) 2003 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_coranswalldistov( regions )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
                            rflo_getnodeoffset
  USE turb_modinterfaces, ONLY : turb_ranswalldistovpatch, &
                                 turb_floextrapintcellscal 
  USE moderror
  USE modmpi
  USE modparameters
  USE turb_modparameters
  IMPLICIT NONE

#include "Indexing.h"

! ... parameters
  TYPE(t_region), POINTER :: regions(:)

! ... loop variables
  INTEGER :: ireg, i, j, k, n, ipatch, ilev

! ... local variables
  TYPE(t_global), POINTER :: global

#ifdef MPI
  INTEGER :: status(mpi_status_size)
#endif
  INTEGER :: ilevbeg,ilevend, ipcbeg,ipcend,jpcbeg,jpcend,kpcbeg,kpcend
  INTEGER :: icoff,ijcoff,inoff,ijnoff, ic,ijkn,ijkni,ijknj,ijknk
  INTEGER :: snddim, rcvdim, ncount, icount, accwalldim
  INTEGER, POINTER :: walldim(:)
  REAL(RFREAL) :: distance, xdist, ydist, zdist, loclens, cdes
  REAL(RFREAL) :: adsi, adsj, adsk, dsi(3), dsj(3), dsk(3)
  REAL(RFREAL), POINTER :: xyz(:,:), cofg(:,:), lens(:)

!******************************************************************************

  global => regions(1)%global
  CALL registerfunction( global,'TURB_CoRansWallDistOV',&
  'TURB_coRansWallDistOVFlo.F90' )

  IF (global%turbWallDim == 0) goto 888

! allocation and initialisation -----------------------------------------------

  ALLOCATE( walldim(global%nRegions) )  ! wall dimension/region, all processors

  IF (.NOT. (global%turbWorkUnused .eqv. .true.)) &
  CALL errorstop( global,err_turb_workspace,__line__ ) 
  global%turbWorkDim = ndir*global%turbWallDim
  ALLOCATE( global%turbWork1D(global%turbWorkDim) )

  global%turbWorkUnused = .false. 
  walldim(:) = 0
  ncount     = 0
  icount     = 0
  accwalldim = 0
  global%turbWorkDim = 0

! search for ns patch at finest level (iLev=1) and store wall xyz in
! global%turbWork1D

  DO ireg=1,global%nRegions

    IF (regions(ireg)%procid==global%myProcid .AND. &   ! region active and
        regions(ireg)%active==active) THEN              ! on my processor
      DO ipatch=1,regions(ireg)%nPatches
        CALL turb_ranswalldistovpatch( regions(ireg), &
                  regions(ireg)%levels(1)%patches(ipatch) )
      ENDDO ! iPatch
      walldim(ireg) = global%turbWorkDim - accwalldim
      accwalldim    = accwalldim + walldim(ireg)
    ENDIF

#ifdef MPI
! - distribut wall dimension acquired from each region/processor to all procs.

    CALL mpi_bcast( walldim(ireg),1,mpi_integer,regions(ireg)%procId, &
                    global%mpiComm,global%mpierr )
    IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )

! - screen output check
!    write(*,*)global%myProcid,iReg,global%turbWorkDim,wallDim(iReg), &
!          NDIR*global%turbWallDim

! - collect wall coordinates (xyz) from all processors to Master processor

    IF (global%myProcid == masterproc) THEN

      IF (walldim(ireg) > 0) THEN
        IF (regions(ireg)%procid /= global%myProcid) THEN

          CALL mpi_recv( global%turbWork1D(1+ncount),walldim(ireg),mpi_rfreal, &
               regions(ireg)%procId,ireg,global%mpiComm,status,global%mpierr )
          IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )

        ENDIF

        ncount = ncount + walldim(ireg)

! ----- screen output check        
!        write(*,*) 'ncount',ireg,ncount

      ENDIF
    ELSE   
      IF (regions(ireg)%procid == global%myProcid) THEN
        IF (walldim(ireg) > 0) THEN

! ------- screen output check
!          write(*,*)'send=',ireg

          CALL mpi_send( global%turbWork1D(1+icount),walldim(ireg),mpi_rfreal, &
                         masterproc,ireg,global%mpiComm,status,global%mpierr )
          IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )

          icount = icount + walldim(ireg)

        ENDIF ! wallDim
      ENDIF   ! region%procId
    ENDIF     ! myProcId
#endif
  ENDDO   ! iReg

! check consistency of total wall dimension

  IF (global%myProcid == masterproc) THEN 
#ifdef MPI
    global%turbWorkDim = ncount
#endif
    IF (global%turbWorkDim /= ndir*global%turbWallDim) &
    CALL errorstop( global,err_turb_ransinput,__line__, &
                 'work array dim. is not consistent with wall array dim.')
  ENDIF

! distribute wall coordinates from Master processor to all

#ifdef MPI
  CALL mpi_bcast( global%turbWork1D,ndir*global%turbWallDim,mpi_rfreal, &
                  masterproc,global%mpiComm,global%mpierr )
  IF (global%mpierr /=0 ) CALL errorstop( global,err_mpi_trouble,__line__ )
#endif

888  CONTINUE

! compute turbulence length scale for RaNS and modified if DES is selected

  DO ireg=1,global%nRegions
    IF (regions(ireg)%procid==global%myProcid .AND. &   ! region active and
        regions(ireg)%active==active) THEN              ! on my processor

      IF (regions(ireg)%mixtInput%moveGrid .eqv. .true.) THEN
        ilevbeg = regions(ireg)%currLevel
        ilevend = regions(ireg)%currLevel
      ELSE
        ilevbeg = 1
        ilevend = regions(ireg)%nGridLevels
      ENDIF

      cdes = regions(ireg)%turbInput%cDes
 
      DO ilev = ilevbeg, ilevend
        CALL rflo_getdimensphys( regions(ireg),ilev,ipcbeg,ipcend, &
                                 jpcbeg,jpcend,kpcbeg,kpcend )
        CALL rflo_getcelloffset( regions(ireg),ilev,icoff,ijcoff )
        CALL rflo_getnodeoffset( regions(ireg),ilev,inoff,ijnoff )

        xyz  => regions(ireg)%levels(ilev)%grid%xyz
        cofg => regions(ireg)%levels(ilev)%grid%cofg
        lens => regions(ireg)%levels(ilev)%turb%lens 
        lens =  global%refLength

! ----- RaNS wall distance providing no-slip wall exist

        DO k=kpcbeg,kpcend
          DO j=jpcbeg,jpcend
            DO i=ipcbeg,ipcend
              ic       = indijk(i,j,k,icoff,ijcoff) 
              lens(ic) = 1.e+32_rfreal
              DO n = 1,ndir*global%turbWallDim,ndir
                xdist     = cofg(xcoord,ic)-global%turbWork1D(n)
                ydist     = cofg(ycoord,ic)-global%turbWork1D(n+1)
                zdist     = cofg(zcoord,ic)-global%turbWork1D(n+2)
                distance  = sqrt( xdist**2 + ydist**2 + zdist**2 )
                lens(ic)  = min( lens(ic),distance )
              ENDDO ! n

! ----------- screen output check
!              IF(i==30 .and.k==2) write(*,*)ireg,j,lens(iC)
            ENDDO   ! i
          ENDDO     ! j
        ENDDO       ! k 

! ----- DES length scale and flows without no-slip walls

        IF ((regions(ireg)%mixtInput%turbModel == turb_model_dessa) .OR. &
            (regions(ireg)%mixtInput%turbModel == turb_model_hdessa) .OR. &
            (global%turbWallDim == 0)) THEN
          DO k=kpcbeg,kpcend
            DO j=jpcbeg,jpcend
              DO i=ipcbeg,ipcend
                ic       = indijk(i  ,j  ,k  ,icoff,ijcoff) 
                ijkn     = indijk(i  ,j  ,k  ,inoff,ijnoff)
                ijkni    = indijk(i+1,j  ,k  ,inoff,ijnoff)
                ijknj    = indijk(i  ,j+1,k  ,inoff,ijnoff)
                ijknk    = indijk(i  ,j  ,k+1,inoff,ijnoff)
                dsi(1:3) = xyz(xcoord:zcoord,ijkni)-xyz(xcoord:zcoord,ijkn)
                dsj(1:3) = xyz(xcoord:zcoord,ijknj)-xyz(xcoord:zcoord,ijkn)
                dsk(1:3) = xyz(xcoord:zcoord,ijknk)-xyz(xcoord:zcoord,ijkn)
                adsi     = sqrt( dsi(1)*dsi(1)+dsi(2)*dsi(2)+dsi(3)*dsi(3) )
                adsj     = sqrt( dsj(1)*dsj(1)+dsj(2)*dsj(2)+dsj(3)*dsj(3) )
                adsk     = sqrt( dsk(1)*dsk(1)+dsk(2)*dsk(2)+dsk(3)*dsk(3) )
                loclens  = max( adsi,adsj,adsk )
                lens(ic) = min( lens(ic),cdes*loclens )
              ENDDO   ! i
            ENDDO     ! j
          ENDDO       ! k 
        ENDIF         ! DES or noWall
      ENDDO           ! iLev

! --- extrapolate solution to dummy cells

      CALL turb_floextrapintcellscal( regions(ireg),lens )

    ENDIF             ! region on this processor and active
  ENDDO               ! iReg

! finalize --------------------------------------------------------------------

  DEALLOCATE( global%turbWork1D, walldim )
  global%turbWorkUnused = .true.

  CALL deregisterfunction( global )

END SUBROUTINE turb_coranswalldistov

!******************************************************************************
!
! RCS Revision history:
!
! $Log: TURB_coRansWallDistOVFlo.F90,v $
! Revision 1.7  2009/08/26 12:28:53  mtcampbe
! Ported to Hera.   Fixed logical expression syntax errors.  Replaced all
! IF (logical_variable)  with IF (logical_variable .eqv. .true.) as
! consistent with the specification.  Also changed: IF( ASSOCIATED(expr) )
! to IF ( ASSOCIATED(expr) .eqv. .true. ).   Intel compilers produce code
! which silently fails for some mal-formed expressions, so these changes
! are a net which should ensure that they are evaluated as intended.
!
! Revision 1.6  2008/12/06 08:44:43  mtcampbe
! Updated license.
!
! Revision 1.5  2008/11/19 22:17:55  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.4  2005/03/09 06:36:11  wasistho
! incorporated HDESSA
!
! Revision 1.3  2004/03/13 03:12:16  wasistho
! get rid of flo/flu identifier in TURB_Co.. routines
!
! Revision 1.2  2004/03/12 02:55:35  wasistho
! changed rocturb routine names
!
! Revision 1.1  2004/03/08 23:33:31  wasistho
! changed turb nomenclature
!
! Revision 1.8  2004/02/13 03:31:27  wasistho
! fixed bug in case a processor has multi regions
!
! Revision 1.7  2004/02/12 03:46:48  wasistho
! filled in RaNS lengthscale in dummy cells
!
! Revision 1.6  2004/01/21 03:47:10  wasistho
! modify k index in screen output check
!
! Revision 1.5  2003/10/20 00:45:42  wasistho
! initiate lens to reference length
!
! Revision 1.4  2003/10/09 20:48:29  wasistho
! added DES lengthscale coefficient CDES
!
! Revision 1.3  2003/10/07 23:55:30  wasistho
! bug fixed missing nodeOffsets
!
! Revision 1.2  2003/10/07 20:32:16  wasistho
! turbWork2D to turbWork1D
!
!
!******************************************************************************