TURB_VisFluxEddyPatch.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 viscous fluxes based on viscosity principle: mu*S_ij through 
!          current patch with mu defined at cell faces
!
! Description: this routine works in the same way as TURB_VisFluxEddy
!              but applied on region patches
!
! Input: region  = data of current region.
!        patch   = current patch.
!
! Output: mixt%diss = total viscous fluxes added to dissipation
!
! Notes: This flux routine is used for the fixed and dynamic Smagorinsky model.
!        This subroutine is similar to subroutine RFLO_ViscousFluxPatch for the 
!        computation of viscous-fluxes at region patches based on mu defined 
!        at cell centers 
!
!******************************************************************************
!
! $Id: TURB_VisFluxEddyPatch.F90,v 1.14 2008/12/06 08:44:42 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_visfluxeddypatch( region,patch )

!$startcopy TURB_vFluxHybridPatch
  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modturbulence, ONLY : t_turb
#ifdef RFLO
  USE modinterfaces, ONLY : rflo_getpatchindices, rflo_getpatchdirection, &
                            rflo_getcelloffset, rflo_getnodeoffset

#include "Indexing.h"
#endif
  USE moderror
  USE modparameters
  USE turb_modparameters
  IMPLICIT NONE

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

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

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

  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend
  INTEGER :: ijk, bctype, indcp, ijkcb0, ijkcd, ijknb, ijknbg
  INTEGER :: n1, n2, noff, j2d, aerocoeff

  REAL(RFREAL)          :: beta,rprt,cpprt,muel,tcol,tcot,engmodel
  REAL(RFREAL)          :: tx,ty,tz,tgradf, bvalfactor
  REAL(RFREAL)          :: velf(3),fd(4),sf(3),sij(3,3)
  REAL(RFREAL)          :: velcb0(3),tfd(4),ucitauij,taukk, ac0,ac1
  REAL(RFREAL)          :: rref, vref, rcfref, rchref
  REAL(RFREAL), POINTER :: diss(:,:),tv(:,:),gv(:,:),muet(:,:),srate(:,:)
  REAL(RFREAL), POINTER :: trace(:),vol(:)

#ifdef RFLO
  INTEGER :: inode, jnode, knode, idir, jdir, kdir
  INTEGER :: ilev, lbound, icoff, ijcoff, inoff, ijnoff, acid0, acid1
  REAL(RFREAL)          :: sgn
  REAL(RFREAL), POINTER :: avgco(:,:), sface(:,:), dv(:,:), grad(:,:)
#endif
#ifdef RFLU
  INTEGER :: ifgbeg
  REAL(RFREAL)          :: rdens0
  REAL(RFREAL), POINTER :: fn(:,:), cv(:,:), grad(:,:,:) 
#endif
!$endcopy

  REAL(RFREAL)          :: muf, modstrain

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

  rcsidentstring = '$RCSfile: TURB_VisFluxEddyPatch.F90,v $'

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

!$startcopy TURB_vFluxHybridPatch
! get dimensions and pointers -------------------------------------------------

  bctype = patch%bcType

#ifdef RFLO
  ilev    =  region%currLevel
  lbound  =  patch%lbound
  vol     => region%levels(ilev)%grid%vol
  dv      => region%levels(ilev)%mixt%dv
  tv      => region%levels(ilev)%mixt%tv
  gv      => region%levels(ilev)%mixt%gv
  diss    => region%levels(ilev)%mixt%diss  
  muet    => region%levels(ilev)%turb%mueT
  trace   => region%levels(ilev)%turb%trace
  turb    => region%levels(ilev)%turb
#endif  
#ifdef RFLU
  vol     => region%grid%vol
  cv      => region%mixt%cv
  tv      => region%mixt%tv
  gv      => region%mixt%gv
  diss    => region%mixt%diss  
  muet    => region%turb%bMueT
  trace   => region%turb%trace
  turb    => region%turb
#endif  

! get coefficients -----------------------------------------------------------

  IF (region%turbInput%engModel==0) THEN
    engmodel = 0._rfreal
  ELSE
    engmodel = 1._rfreal
  ENDIF
  beta  = region%mixtInput%betrk(region%irkStep)

  IF (bctype>=bc_noslipwall .AND. bctype<=bc_noslipwall+bc_range) THEN
    bvalfactor = 0._rfreal
  ELSE
    bvalfactor = 1._rfreal
  ENDIF
  
#ifdef RFLO
  rprt      = 1._rfreal/region%levels(ilev)%mixt%prTurb
  indcp     = region%levels(ilev)%mixt%indCp
  aerocoeff = global%aeroCoeffs
  noff      = abs(patch%l1end-patch%l1beg) + 1
#endif
#ifdef RFLU
  rprt      = 1._rfreal/region%mixtInput%prTurb
  indcp     = region%mixtInput%indCp
#endif

  rref   = global%refDensity
  vref   = global%refVelocity
  
  rcfref = 2.0_rfreal/(rref*vref*vref)
  rchref = 2.0_rfreal/(rref*vref*vref*vref)

#ifdef RFLO

! take the right face vector and make it point outwards

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

  sgn   = +1._rfreal
  inode = 0
  jnode = 0
  knode = 0
  acid0 = 2
  acid1 = 1
  IF (lbound==2 .OR. lbound==4 .OR. lbound==6) THEN
    sgn   = -1._rfreal
    inode = -idir
    jnode = -jdir
    knode = -kdir
    acid0 = 1
    acid1 = 2
  ENDIF
!$endcopy

! get the appropriate face vector

  IF (lbound==1 .OR. lbound==2) THEN
    ijk   =  diri
    avgco => region%levels(ilev)%grid%c2fCoI
    sface => region%levels(ilev)%grid%si
    grad  => region%levels(ilev)%mixt%gradi
    srate => region%levels(ilev)%turb%mISij
  ELSE IF (lbound==3 .OR. lbound==4) THEN
    ijk   =  dirj
    avgco => region%levels(ilev)%grid%c2fCoJ
    sface => region%levels(ilev)%grid%sj
    grad  => region%levels(ilev)%mixt%gradj
    srate => region%levels(ilev)%turb%mJSij
  ELSE
    ijk   =  dirk
    avgco => region%levels(ilev)%grid%c2fCoK
    sface => region%levels(ilev)%grid%sk
    grad  => region%levels(ilev)%mixt%gradk
    srate => region%levels(ilev)%turb%mKSij
  ENDIF

!$startcopy TURB_vFluxHybridPatch

! non-conforming region interface

  IF (bctype>=bc_regionint .AND. bctype<=bc_regionint+bc_range) THEN

  ELSE IF (bctype>=bc_regnonconf .AND. bctype<=bc_regnonconf+bc_range) THEN

! everything else

  ELSE

! flux in the direction normal to the patch

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

          ijkcb0 = indijk(i      ,j      ,k      ,icoff,ijcoff)  ! bnd cells
          ijkcd  = indijk(i-idir ,j-jdir ,k-kdir ,icoff,ijcoff)  ! dummy
          ijknb  = indijk(i+inode,j+jnode,k+knode,inoff,ijnoff)  ! bnd nodes
          ijknbg = ijknb

          sf(1)  = sgn*sface(xcoord,ijknb)
          sf(2)  = sgn*sface(ycoord,ijknb)
          sf(3)  = sgn*sface(zcoord,ijknb)
          ac0 = avgco(acid0,ijknb)
          ac1 = avgco(acid1,ijknb)
#endif
#ifdef RFLU
! specific Rocflu, check the state of cv first
  IF (region%mixt%cvState /= cv_mixt_state_cons) &
                            CALL errorstop(global,err_reached_default,__line__)
  ijk    = diri
  ibeg   = 1
  iend   = patch%nBFaces
! TEMPORARY : removing usage of bf2bg from everywhere
!  ifgBeg = patch%bf2bg(BF2BG_BEG)
  ac0    = 0.5_rfreal
  ac1    = 0.5_rfreal

  fn    => patch%fn
! TEMPORARY 
  grad  => patch%mixt%gradFace
  srate => region%turb%bmISij  

  DO ic=ibeg,iend
          ijkcb0 = patch%bf2c(ic)
          ijkcd  = ijkcb0
          ijknb  = ic
          ijknbg = ic + ifgbeg-1

          sf(1)  = fn(xcoord,ijknb)*fn(xyzmag,ijknb)
          sf(2)  = fn(ycoord,ijknb)*fn(xyzmag,ijknb)
          sf(3)  = fn(zcoord,ijknb)*fn(xyzmag,ijknb)           
          rdens0 = 1._rfreal/cv(cv_mixt_dens,ijkcb0)
#endif

!pert TURB
          cpprt = (ac0*gv(gv_mixt_cp,ijkcb0*indcp) + &
                   ac1*gv(gv_mixt_cp,ijkcd*indcp))*rprt
          muel  = ac0*tv(tv_mixt_muel,ijkcb0)+ac1*tv(tv_mixt_muel,ijkcd)
          tcol  = ac0*tv(tv_mixt_tcol,ijkcb0)+ac1*tv(tv_mixt_tcol,ijkcd)
          tcot  = cpprt*muet(ijk,ijknbg)
!endpert

#ifdef RFLO
          velcb0(1)= dv(dv_mixt_uvel,ijkcb0)
          velcb0(2)= dv(dv_mixt_vvel,ijkcb0)
          velcb0(3)= dv(dv_mixt_wvel,ijkcb0)
          velf(1)= ac0*dv(dv_mixt_uvel,ijkcb0)+ac1*dv(dv_mixt_uvel,ijkcd)
          velf(2)= ac0*dv(dv_mixt_vvel,ijkcb0)+ac1*dv(dv_mixt_vvel,ijkcd)
          velf(3)= ac0*dv(dv_mixt_wvel,ijkcb0)+ac1*dv(dv_mixt_wvel,ijkcd)
          tx  = grad(gr_mixt_tx,ijknb)
          ty  = grad(gr_mixt_ty,ijknb)
          tz  = grad(gr_mixt_tz,ijknb)
#endif
#ifdef RFLU
          velcb0(1)= cv(cv_mixt_xmom,ijkcb0)*rdens0
          velcb0(2)= cv(cv_mixt_ymom,ijkcb0)*rdens0
          velcb0(3)= cv(cv_mixt_zmom,ijkcb0)*rdens0
          velf(1)  = velcb0(1)
          velf(2)  = velcb0(2)
          velf(3)  = velcb0(3)
          velf(1:3)= bvalfactor*velf(1:3)
          tx  = grad(xcoord,grf_mixt_temp,ijknbg)
          ty  = grad(ycoord,grf_mixt_temp,ijknbg)
          tz  = grad(zcoord,grf_mixt_temp,ijknbg)
#endif

          tgradf = tx*sf(1)+ty*sf(2)+tz*sf(3)

          sij(1,1) = srate(e11,ijknbg)
          sij(1,2) = srate(e12,ijknbg)
          sij(1,3) = srate(e13,ijknbg)

          sij(2,1) = sij(1,2)
          sij(2,2) = srate(e22,ijknbg)
          sij(2,3) = srate(e23,ijknbg)

          sij(3,1) = sij(1,3)
          sij(3,2) = sij(2,3)
          sij(3,3) = srate(e33,ijknbg)

          fd(1) = muel*(sij(1,1)*sf(1)+sij(1,2)*sf(2)+sij(1,3)*sf(3))
          fd(2) = muel*(sij(2,1)*sf(1)+sij(2,2)*sf(2)+sij(2,3)*sf(3))
          fd(3) = muel*(sij(3,1)*sf(1)+sij(3,2)*sf(2)+sij(3,3)*sf(3))
          fd(4) = dot_product(fd(1:3),velf(1:3)) + tcol*tgradf
!wlmCheckprobe--------------------------------------------------------------
!  IF (bcType>=BC_NOSLIPWALL .AND. bcType<=BC_NOSLIPWALL+BC_RANGE) THEN
!    write(*,*) region%procId, i, j, k, mueL*sij(1,2), mueL*sij(3,2), &
!               mueL*sij(2,1), mueL*sij(2,3), tCoL*tgradf
!  ENDIF
!---------------------------------------------------------------------------
!$endcopy
! Yoshizawa k-model for eddy viscosity type LES
          modstrain = sqrt(sij(1,1)*sij(1,1) &
                          +sij(2,2)*sij(2,2) &
                          +sij(1,1)*sij(2,2) &
                          +sij(1,2)*sij(1,2) &
                          +sij(1,3)*sij(1,3) &
                          +sij(2,3)*sij(2,3))
          taukk = 2._rfreal*muet(ijk,ijknbg)*modstrain
          trace(ijkcb0) = trace(ijkcb0) + taukk

          muf      = -muet(ijk,ijknbg)/muel
          tfd(1)   = muf*fd(1)
          tfd(2)   = muf*fd(2)
          tfd(3)   = muf*fd(3)
          ucitauij = dot_product(tfd(1:3),velcb0(1:3))
          tfd(4)   = ucitauij - tcot*tgradf

!$startcopy TURB_vFluxHybridPatch
          fd(1) = fd(1)-tfd(1)
          fd(2) = fd(2)-tfd(2)
          fd(3) = fd(3)-tfd(3)
          fd(4) = fd(4)-engmodel*tfd(4)

#ifdef STATS
          turb%stwork(1,ijkcb0)=turb%stwork(1,ijkcb0)+ucitauij    ! alp_1
          turb%stwork(2,ijkcb0)=turb%stwork(2,ijkcb0)+tcot*tgradf ! alp_2+alp_3
#endif
          diss(cv_mixt_xmom,ijkcb0) = diss(cv_mixt_xmom,ijkcb0)+fd(1)*beta
          diss(cv_mixt_ymom,ijkcb0) = diss(cv_mixt_ymom,ijkcb0)+fd(2)*beta
          diss(cv_mixt_zmom,ijkcb0) = diss(cv_mixt_zmom,ijkcb0)+fd(3)*beta
          diss(cv_mixt_ener,ijkcb0) = diss(cv_mixt_ener,ijkcb0)+fd(4)*beta
          global%esg1Sum  = global%esg1Sum + vol(ijkcb0)*ucitauij

! ------- Set friction and heat-transfer coefficients

#ifdef RFLO
          IF      (lbound==1 .OR. lbound==2) THEN
            n1 = j - jbeg
            n2 = k - kbeg
          ELSE IF (lbound==3 .OR. lbound==4) THEN
            n1 = k - kbeg
            n2 = i - ibeg
          ELSE IF (lbound==5 .OR. lbound==6) THEN
            n1 = i - ibeg
            n2 = j - jbeg
          ENDIF
          j2d  = aerocoeff * indij(n1,n2,noff)
          patch%cf(xcoord,j2d) = rcfref*fd(1)
          patch%cf(ycoord,j2d) = rcfref*fd(2)
          patch%cf(zcoord,j2d) = rcfref*fd(3)
          patch%ch(       j2d) = rchref*fd(4)

        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k
  ENDIF         ! bcType
#endif
#ifdef RFLU
          patch%cf(xcoord,ic) = rcfref*fd(1)
          patch%cf(ycoord,ic) = rcfref*fd(2)
          patch%cf(zcoord,ic) = rcfref*fd(3)
          patch%ch(       ic) = rchref*fd(4)               
  ENDDO         ! iC
#endif

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

  CALL deregisterfunction( global )
!$endcopy
END SUBROUTINE turb_visfluxeddypatch

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: TURB_VisFluxEddyPatch.F90,v $
! Revision 1.14  2008/12/06 08:44:42  mtcampbe
! Updated license.
!
! Revision 1.13  2008/11/19 22:17:55  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.12  2006/08/19 15:40:40  mparmar
! Removed bf2bg and moved region%mixt%bGradFace to patch%mixt%gradFace
!
! Revision 1.11  2006/03/13 03:43:23  wasistho
! defined aero coeffs
!
! Revision 1.10  2004/10/22 23:20:39  wasistho
! collect esg1 and esg2+3 into statistics
!
! Revision 1.9  2004/08/04 22:07:19  wasistho
! bugfixed: ac0 and ac1 are common to flo and flu
!
! Revision 1.8  2004/08/02 23:06:36  wasistho
! mv libfloflu/viscousFluxEddy(Patch) to rocflo/RFLO_viscousFlux(Patch)
!
! Revision 1.7  2004/08/02 21:06:17  wasistho
! moved averaging coeff local variables to within ifdef RFLO
!
! Revision 1.6  2004/08/02 19:34:18  wasistho
! changed grid%avgCo to grid%c2fCo
!
! Revision 1.5  2004/07/31 00:53:38  wasistho
! replaced cell2face midpoint by linear averaging
!
! Revision 1.4  2004/03/27 02:16:42  wasistho
! compiled with Rocflu
!
! Revision 1.3  2004/03/25 04:40:41  wasistho
! prepared for RFLU
!
! Revision 1.2  2004/03/24 03:37:03  wasistho
! prepared for RFLU
!
! Revision 1.1  2004/03/05 04:37:00  wasistho
! changed nomenclature
!
! Revision 1.7  2004/02/26 21:22:23  wasistho
! changed turb%esg.. to global%esg..
!
! Revision 1.6  2004/01/21 03:42:18  wasistho
! get rid of grad(T)*n = 0 for adiabatic wall
!
! Revision 1.5  2003/08/29 01:43:17  wasistho
! Added TARGET attribute to region variable, since pointers are cached into it
!
! Revision 1.4  2003/06/05 19:21:16  wasistho
! adiabatic switch
!
! Revision 1.3  2003/05/31 01:48:23  wasistho
! installed turb. wall layer model
!
! Revision 1.2  2003/05/15 02:57:06  jblazek
! Inlined index function.
!
! Revision 1.1  2002/10/14 23:55:30  wasistho
! Install Rocturb
!
!******************************************************************************