TURB_floWlmMetric.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 metric coefficients needed in wlm computations.
!
! Description: Metric quantities stored in data structure are,
!              - transformation matrix coeffs. from body fitted to Cartesian
!              - grid spacing in body fitted coord. in wall parallel 
!                directions at first (or user selected) points from the wall
!              - wall normal distance to first points from the wall.
!
! Input: region  = data of current region.
!        patch   = current patch.
!
! Output: metric quantities described above.
!
! Notes: none.
!
!******************************************************************************
!
! $Id: TURB_floWlmMetric.F90,v 1.5 2008/12/06 08:44:44 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_flowlmmetric( region,patch )

  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modinterfaces, ONLY : rflo_getpatchindices,rflo_getpatchdirection, &
                            rflo_getnodeoffset, rflo_getcelloffset
  USE moderror
  USE modparameters
  USE turb_modparameters
  IMPLICIT NONE

#include "Indexing.h"

! ... parameters
  TYPE(t_region) :: region
  TYPE(t_patch)  :: patch

! ... loop variables
  INTEGER :: i, j, k, l

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

  INTEGER :: ilev, lbound, nref, isub, jsub, ksub
  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend
  INTEGER :: inode, jnode, knode, idir, jdir, kdir, ir, jr, kr
  INTEGER :: icoff, ijcoff, inoff, ijnoff, ix,jx,kx, iz,jz,kz, ixz,jxz,kxz
  INTEGER :: ijkn1b, ijkn1, ijkn1x, ijkn1z, ijkval
  INTEGER :: ijkn2b, ijkn2, ijkn2x, ijkn2z
  INTEGER :: ijkn1bx, ijkn1bz, ijkn1bxz, ijkn2bx, ijkn2bz, ijkn2bxz
  INTEGER :: ijkn0b, ijkn0bx, ijkn0bz, ijkn0bxz

  REAL(RFREAL)          :: sxi(3), set(3), szt(3), fcenb(3), fcene(3), xyz(3,4)
  REAL(RFREAL)          :: modsxi, modset, modszt, sgn
  REAL(RFREAL), POINTER :: sfxi(:,:), sfet(:,:), sfzt(:,:)

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

  rcsidentstring = '$RCSfile: TURB_floWlmMetric.F90,v $ $Revision: 1.5 $'

  global => region%global
  CALL registerfunction( global,'TURB_FloWlmMetric',&
  'TURB_floWlmMetric.F90' )

! get dimensions, parameters and pointers -------------------------------------

  ilev    =  region%currLevel
  lbound  =  patch%lbound
  nref    =  patch%valBola%switches(wlm_input_refpoint) - 1

  CALL rflo_getpatchindices( region,patch,ilev,ibeg,iend,jbeg,jend,kbeg,kend )
  CALL rflo_getpatchdirection( patch,idir,jdir,kdir )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )

  sgn  = +1._rfreal
  inode = 0
  jnode = 0
  knode = 0
  IF (lbound==2 .OR. lbound==4 .OR. lbound==6) THEN
    sgn   = -1._rfreal
    inode = -idir
    jnode = -jdir
    knode = -kdir
  ENDIF
  
  IF (lbound==1 .OR. lbound==2) THEN
    ir = idir*nref
    jr = 0
    kr = 0
    ix = 0
    jx = 0
    kx = 1
    iz = 0
    jz = 1
    kz = 0
    sfxi => region%levels(ilev)%grid%sk
    sfet => region%levels(ilev)%grid%si
    sfzt => region%levels(ilev)%grid%sj
  ELSEIF (lbound==3 .OR. lbound==4) THEN
    ir = 0
    jr = jdir*nref
    kr = 0
    ix = 1
    jx = 0
    kx = 0
    iz = 0
    jz = 0
    kz = 1
    sfxi => region%levels(ilev)%grid%si
    sfet => region%levels(ilev)%grid%sj
    sfzt => region%levels(ilev)%grid%sk
  ELSEIF (lbound==5 .OR. lbound==6) THEN
    ir = 0
    jr = 0
    kr = kdir*nref
    ix = 0
    jx = 1
    kx = 0
    iz = 1
    jz = 0
    kz = 0
    sfxi => region%levels(ilev)%grid%sj
    sfet => region%levels(ilev)%grid%sk
    sfzt => region%levels(ilev)%grid%si
  ENDIF
  ixz = ix+iz
  jxz = jx+jz
  kxz = kx+kz

! Get mapping coefficients : define xi,eta,zeta normal vectors from i,j,k 
! face vectors at cell center. The fact that i,j,k face vectors points to 
! negative i,j,k directions are taken into account. Next, get finite 
! difference spacings

  DO k=kbeg,kend
    DO j=jbeg,jend
      DO i=ibeg,iend

! ----- get mapping coefficients

        isub = i+ir
        jsub = j+jr
        ksub = k+kr

        ijkn1b = indijk(isub+inode,jsub+jnode,ksub+knode,inoff,ijnoff) !bnd nodes
        ijkn1  = indijk(isub      ,jsub      ,ksub      ,inoff,ijnoff)  
        ijkn1x = indijk(isub+ix   ,jsub+jx   ,ksub+kx   ,inoff,ijnoff)  
        ijkn1z = indijk(isub+iz   ,jsub+jz   ,ksub+kz   ,inoff,ijnoff)  

        ijkn2b = &
        indijk(isub+inode+idir,jsub+jnode+jdir,ksub+knode+kdir,inoff,ijnoff) 

        IF (lbound==1 .OR. lbound==2) THEN
          ijkval = indij(j-jbeg ,k-kbeg ,jend-jbeg+1)
        ELSEIF (lbound==3 .OR. lbound==4) THEN
          ijkval = indij(k-kbeg ,i-ibeg ,kend-kbeg+1)
        ELSEIF (lbound==5 .OR. lbound==6) THEN
          ijkval = indij(i-ibeg ,j-jbeg ,iend-ibeg+1)
        ENDIF

        DO l=xcoord,zcoord
          sxi(l) = -0.5_rfreal*(sfxi(l,ijkn1)+sfxi(l,ijkn1x))
          set(l) = -sgn*0.5_rfreal*(sfet(l,ijkn1b)+sfet(l,ijkn2b))
          szt(l) = -0.5_rfreal*(sfzt(l,ijkn1)+sfzt(l,ijkn1z))
        ENDDO

        modsxi = sqrt( sxi(1)*sxi(1)+sxi(2)*sxi(2)+sxi(3)*sxi(3) )
        modset = sqrt( set(1)*set(1)+set(2)*set(2)+set(3)*set(3) )
        modszt = sqrt( szt(1)*szt(1)+szt(2)*szt(2)+szt(3)*szt(3) )

        patch%valBola%vals(ijkval,wlm_vals_xix) = sxi(1)/modsxi
        patch%valBola%vals(ijkval,wlm_vals_etx) = set(1)/modset
        patch%valBola%vals(ijkval,wlm_vals_ztx) = szt(1)/modszt

        patch%valBola%vals(ijkval,wlm_vals_xiy) = sxi(2)/modsxi
        patch%valBola%vals(ijkval,wlm_vals_ety) = set(2)/modset
        patch%valBola%vals(ijkval,wlm_vals_zty) = szt(2)/modszt

        patch%valBola%vals(ijkval,wlm_vals_xiz) = sxi(3)/modsxi
        patch%valBola%vals(ijkval,wlm_vals_etz) = set(3)/modset
        patch%valBola%vals(ijkval,wlm_vals_ztz) = szt(3)/modszt

! ----- get finite difference spacings dXi, dZeta, respectively

        isub = i+ir+inode
        jsub = j+jr+jnode
        ksub = k+kr+knode

        ijkn1bx = indijk(isub+ix  ,jsub+jx  ,ksub+kx  ,inoff,ijnoff) 
        ijkn1bz = indijk(isub+iz  ,jsub+jz  ,ksub+kz  ,inoff,ijnoff) 
        ijkn1bxz= indijk(isub+ixz ,jsub+jxz ,ksub+kxz ,inoff,ijnoff)

        ijkn2bx = &
        indijk(isub+idir+ix  ,jsub+jdir+jx  ,ksub+kdir+kx  ,inoff,ijnoff)  
        ijkn2bz = &
        indijk(isub+idir+iz  ,jsub+jdir+jz  ,ksub+kdir+kz  ,inoff,ijnoff)  
        ijkn2bxz= &
        indijk(isub+idir+ixz ,jsub+jdir+jxz ,ksub+kdir+kxz ,inoff,ijnoff)  

        DO l=xcoord,zcoord
          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn1b)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn2b)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn1bz)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn2bz)
          fcenb(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4))

          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn1bx)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn2bx)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn1bxz)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn2bxz)
          fcene(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4))
        ENDDO
        patch%valBola%vals(ijkval,wlm_vals_dxi) = sqrt( (fcene(1)-fcenb(1))**2+ &
                                                        (fcene(2)-fcenb(2))**2+ &
                                                        (fcene(3)-fcenb(3))**2 )

        DO l=xcoord,zcoord
          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn1b)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn2b)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn1bx)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn2bx)
          fcenb(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4))

          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn1bz)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn2bz)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn1bxz)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn2bxz)
          fcene(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4))
        ENDDO
        patch%valBola%vals(ijkval,wlm_vals_dzt) = sqrt( (fcene(1)-fcenb(1))**2+ &
                                                        (fcene(2)-fcenb(2))**2+ &
                                                        (fcene(3)-fcenb(3))**2 )

! ----- get wall distance wDist

        isub = i+inode
        jsub = j+jnode
        ksub = k+knode

        ijkn0b   = indijk(isub     ,jsub     ,ksub     ,inoff,ijnoff)
        ijkn0bx  = indijk(isub+ix  ,jsub+jx  ,ksub+kx  ,inoff,ijnoff)
        ijkn0bz  = indijk(isub+iz  ,jsub+jz  ,ksub+kz  ,inoff,ijnoff)
        ijkn0bxz = indijk(isub+ixz ,jsub+jxz ,ksub+kxz ,inoff,ijnoff)

        DO l=xcoord,zcoord
          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn1b)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn1bx)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn1bz)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn1bxz)
          fcenb(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4))
          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn2b)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn2bx)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn2bz)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn2bxz)
          fcene(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4)) 
          fcene(l) = 0.50_rfreal*(fcene(l)+fcenb(l))
          xyz(l,1) = region%levels(ilev)%grid%xyz(l,ijkn0b)
          xyz(l,2) = region%levels(ilev)%grid%xyz(l,ijkn0bx)
          xyz(l,3) = region%levels(ilev)%grid%xyz(l,ijkn0bz)
          xyz(l,4) = region%levels(ilev)%grid%xyz(l,ijkn0bxz)
          fcenb(l) = 0.25_rfreal*(xyz(l,1) + xyz(l,2) + xyz(l,3) + xyz(l,4))
          set(l)   = -sgn*0.5_rfreal*(sfet(l,ijkn1b)+sfet(l,ijkn2b))
        ENDDO
        modset = sqrt( set(1)*set(1)+set(2)*set(2)+set(3)*set(3) )
        patch%valBola%vals(ijkval,wlm_vals_wdist) = &
                      abs( (fcene(1)-fcenb(1))*set(1)+ &
                           (fcene(2)-fcenb(2))*set(2)+ &
                           (fcene(3)-fcenb(3))*set(3) )/modset

      ENDDO   ! i
    ENDDO     ! j
  ENDDO       ! k

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

  CALL deregisterfunction( global )

END SUBROUTINE turb_flowlmmetric

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: TURB_floWlmMetric.F90,v $
! Revision 1.5  2008/12/06 08:44:44  mtcampbe
! Updated license.
!
! Revision 1.4  2008/11/19 22:17:56  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.3  2004/03/19 02:53:16  wasistho
! prepared for RFLU
!
! Revision 1.2  2004/03/12 02:55:36  wasistho
! changed rocturb routine names
!
! Revision 1.1  2004/03/08 23:35:46  wasistho
! changed turb nomenclature
!
! Revision 1.3  2004/02/14 03:42:33  wasistho
! added new WLM parameter: reference point
!
!
!******************************************************************************