TURB_VisFluxEddy.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 with
!          mu defined at cell faces
!
! Description: The fluxes are computed in three consecutive sweeps (i, j and k
!              direction) at interior faces. The fluxes at the boundary patches
!              are then treated. The energy model contributions are also
!              treated on the fly.
!
! Input: region  = data of current region
!
! Output: mixt%diss = total viscous fluxes added to dissipation
!
! Notes: This flux routine is used for the fixed and dynamic Smagorinsky model.
!        This routine is similar to subroutine RFLO_ViscousFlux for the 
!        computation of viscous-fluxes based on mu defined at cell centers. 
!
!******************************************************************************
!
! $Id: TURB_VisFluxEddy.F90,v 1.11 2008/12/06 08:44:42 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_visfluxeddy( region )

!$startcopy TURB_VFluxHybrid
  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modglobal, ONLY     : t_global
  USE modturbulence, ONLY : t_turb
!$endcopy
  USE turb_modinterfaces, ONLY : turb_lesesgmodel4, turb_visfluxeddypatch, &
                                 turb_wlmtauwallmapping, turb_wlmupdate, &
                                 turb_wlmfluxpatch, turb_statccollector
!$startcopy TURB_VFluxHybrid
#ifdef RFLO
  USE modinterfaces, ONLY : rflo_getdimensphys, &
                            rflo_getcelloffset, rflo_getnodeoffset

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

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

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

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

  INTEGER :: ijkc0,ijkc1,ijkn, indcp, ibn,ien
  REAL(RFREAL)          :: one6th,beta,engmodel,rprt,cpprt,muel,tcol,tcot
  REAL(RFREAL)          :: muf,tx,ty,tz,tgradf
  REAL(RFREAL)          :: velf(3),fd(4),sface(3),sij(3,3)
  REAL(RFREAL)          :: velc0(3),velc1(3),tfd(4),fd4,ucitauij,taukk
  REAL(RFREAL), POINTER :: diss(:,:),tv(:,:),gv(:,:),muet(:,:),srate(:,:)
  REAL(RFREAL), POINTER :: trace(:),vol(:)
  LOGICAL               :: dowlm

#ifdef RFLO
  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ilev,icoff,ijcoff,inoff,ijnoff
  REAL(RFREAL), POINTER :: avgco(:,:), sf(:,:), dv(:,:), grad(:,:)
#endif
#ifdef RFLU
  INTEGER, POINTER      :: f2c(:,:)
  REAL(RFREAL)          :: rdens0, rdens1
  REAL(RFREAL), POINTER :: fn(:,:), cv(:,:), grad(:,:,:)
#endif
!$endcopy
!******************************************************************************

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

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

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

#ifdef RFLO
  ilev    =  region%currLevel
  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
  vol     => region%levels(ilev)%grid%vol
  turb    => region%levels(ilev)%turb

  CALL rflo_getdimensphys( region,ilev,ipcbeg,ipcend, &
                           jpcbeg,jpcend,kpcbeg,kpcend )
  CALL rflo_getcelloffset( region,ilev,icoff,ijcoff )
  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
#endif
#ifdef RFLU
  cv      => region%mixt%cv
  tv      => region%mixt%tv
  gv      => region%mixt%gv
  diss    => region%mixt%diss  
  muet    => region%turb%mueT
  trace   => region%turb%trace
  vol     => region%grid%vol
  turb    => region%turb
#endif

! get coefficients and parameters

  one6th = 1._rfreal/6._rfreal

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

#ifdef RFLO
  rprt  = 1._rfreal/region%levels(ilev)%mixt%prTurb
  indcp = region%levels(ilev)%mixt%indCp
#endif
#ifdef RFLU
  rprt  = 1._rfreal/region%mixtInput%prTurb
  indcp = region%mixtInput%indCp
#endif

! compute interior fluxes through I, J and K faces

#ifdef STATS
! if desired, collect quantities of interest; for that, allocate work space
  ibn = lbound( trace,1 )
  ien = ubound( trace,1 )
  ALLOCATE( turb%stwork(2,ibn:ien) )
  turb%stwork = 0._rfreal
#endif

  CALL computefluxtot( diri )
#ifdef RFLO
  CALL computefluxtot( dirj )
  CALL computefluxtot( dirk )
#endif

! get fluxes through boundaries
#ifdef RFLO  
  DO ipatch=1,region%nPatches
    patch => region%levels(ilev)%patches(ipatch)
#endif
#ifdef RFLU  
  DO ipatch=1,region%grid%nPatches
    patch => region%patches(ipatch)
#endif

    dowlm = .false.
#ifdef RFLO
    IF (patch%bcType>=bc_noslipwall .AND. &
        patch%bcType<=bc_noslipwall+bc_range) THEN   ! my boundary type
      IF (patch%valBola%switches(wlm_input_model) /= wlm_model_nomodel) THEN
        dowlm = .true.
      ENDIF
    ENDIF
#endif
!$endcopy

    IF (dowlm) THEN
      CALL turb_wlmupdate( region,patch )
      CALL turb_wlmtauwallmapping( region,patch )  
      CALL turb_wlmfluxpatch( region,patch )
    ELSE
      CALL turb_visfluxeddypatch( region,patch )
    ENDIF
  ENDDO

!$startcopy TURB_VFluxHybrid
#ifdef STATS
! collect statistics and deallocate work space used

  IF ((region%turbInput%nSt > 0) .AND. (region%turbInput%engModel/=off)) THEN
    turb%stwork(:,:) = one6th*turb%stwork(:,:)
    CALL turb_statccollector( region,1,2,turb%stwork )
  ENDIF
  DEALLOCATE( turb%stwork )
#endif

! if active, get contribution of energy subgrid model alpha_4
  IF (region%turbInput%engModel/=off) THEN
    CALL turb_lesesgmodel4( region )
  ENDIF

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

  CALL deregisterfunction( global )

! =============================================================================
!   Flux computation subroutine
! =============================================================================

CONTAINS

  SUBROUTINE computefluxtot( ijk )

! ... parameters
    INTEGER   :: ijk

! ... local variables
    INTEGER   :: ibeg,iend,jbeg,jend,kbeg,kend, iadd,jadd,kadd
    REAL(RFREAL) :: ac0, ac1
!$endcopy
!pert TURB 
    REAL(RFREAL) :: modstrain
!endpert

!$startcopy TURB_VFluxHybrid
! - Set limits and pointers ---------------------------------------------------

#ifdef RFLO
    IF (ijk==diri) THEN
      ibeg = ipcbeg+1
      iend = ipcend
      jbeg = jpcbeg
      jend = jpcend
      kbeg = kpcbeg
      kend = kpcend
      iadd = -1
      jadd = 0
      kadd = 0
      srate => region%levels(ilev)%turb%mISij
      grad  => region%levels(ilev)%mixt%gradi
      sf    => region%levels(ilev)%grid%si
      avgco => region%levels(ilev)%grid%c2fCoI
    ELSEIF (ijk==dirj) THEN
      ibeg = ipcbeg
      iend = ipcend
      jbeg = jpcbeg+1
      jend = jpcend
      kbeg = kpcbeg
      kend = kpcend
      iadd = 0
      jadd = -1
      kadd = 0
      srate => region%levels(ilev)%turb%mJSij
      grad  => region%levels(ilev)%mixt%gradj
      sf    => region%levels(ilev)%grid%sj
      avgco => region%levels(ilev)%grid%c2fCoJ
    ELSEIF (ijk==dirk) THEN
      ibeg = ipcbeg
      iend = ipcend
      jbeg = jpcbeg
      jend = jpcend
      kbeg = kpcbeg+1
      kend = kpcend
      iadd = 0
      jadd = 0
      kadd = -1
      srate => region%levels(ilev)%turb%mKSij
      grad  => region%levels(ilev)%mixt%gradk
      sf    => region%levels(ilev)%grid%sk
      avgco => region%levels(ilev)%grid%c2fCoK
    ENDIF
#endif
#ifdef RFLU
    ibeg =  1
    iend =  region%grid%nFaces
    f2c   => region%grid%f2c
    srate => region%turb%mISij 
    grad  => region%mixt%gradFace
    fn    => region%grid%fn
#endif

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

          ijkc0 = indijk(i     ,j     ,k     ,icoff,ijcoff)
          ijkc1 = indijk(i+iadd,j+jadd,k+kadd,icoff,ijcoff)
          ijkn  = indijk(i     ,j     ,k     ,inoff,ijnoff)
          sface(1)= sf(xcoord,ijkn)
          sface(2)= sf(ycoord,ijkn)
          sface(3)= sf(zcoord,ijkn)
          ac0 = avgco(2,ijkn)
          ac1 = avgco(1,ijkn)
#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__)
    ac0 = 0.5_rfreal
    ac1 = 0.5_rfreal
    DO ic=ibeg,iend
          ijkc0 = f2c(1,ic)
          ijkc1 = f2c(2,ic)
          ijkn  = ic
          sface(1)= fn(xcoord,ijkn)*fn(xyzmag,ijkn)
          sface(2)= fn(ycoord,ijkn)*fn(xyzmag,ijkn)
          sface(3)= fn(zcoord,ijkn)*fn(xyzmag,ijkn)           
          rdens0  = 1._rfreal/cv(cv_mixt_dens,ijkc0)
          rdens1  = 1._rfreal/cv(cv_mixt_dens,ijkc1)
#endif
          cpprt = (ac0*gv(gv_mixt_cp,ijkc0*indcp) + &
                   ac1*gv(gv_mixt_cp,ijkc1*indcp))*rprt
          muel  = ac0*tv(tv_mixt_muel,ijkc0)+ac1*tv(tv_mixt_muel,ijkc1)
          tcol  = ac0*tv(tv_mixt_tcol,ijkc0)+ac1*tv(tv_mixt_tcol,ijkc1)
          tcot  = cpprt*muet(ijk,ijkn)
#ifdef RFLO
          velc0(1)= dv(dv_mixt_uvel,ijkc0)
          velc0(2)= dv(dv_mixt_vvel,ijkc0)
          velc0(3)= dv(dv_mixt_wvel,ijkc0)
          velc1(1)= dv(dv_mixt_uvel,ijkc1)
          velc1(2)= dv(dv_mixt_vvel,ijkc1)
          velc1(3)= dv(dv_mixt_wvel,ijkc1)
          velf(1) = ac0*dv(dv_mixt_uvel,ijkc0)+ac1*dv(dv_mixt_uvel,ijkc1)
          velf(2) = ac0*dv(dv_mixt_vvel,ijkc0)+ac1*dv(dv_mixt_vvel,ijkc1)
          velf(3) = ac0*dv(dv_mixt_wvel,ijkc0)+ac1*dv(dv_mixt_wvel,ijkc1)
          tx  = grad(gr_mixt_tx,ijkn)
          ty  = grad(gr_mixt_ty,ijkn)
          tz  = grad(gr_mixt_tz,ijkn)
#endif
#ifdef RFLU
          velc0(1)= cv(cv_mixt_xmom,ijkc0)*rdens0
          velc0(2)= cv(cv_mixt_ymom,ijkc0)*rdens0
          velc0(3)= cv(cv_mixt_zmom,ijkc0)*rdens0
          velc1(1)= cv(cv_mixt_xmom,ijkc1)*rdens1
          velc1(2)= cv(cv_mixt_ymom,ijkc1)*rdens1
          velc1(3)= cv(cv_mixt_zmom,ijkc1)*rdens1
          velf(1) = ac0*velc0(1)+ac1*velc1(1)
          velf(2) = ac0*velc0(2)+ac1*velc1(2)
          velf(3) = ac0*velc0(3)+ac1*velc1(3)
          tx  = grad(xcoord,grf_mixt_temp,ijkn)
          ty  = grad(ycoord,grf_mixt_temp,ijkn)
          tz  = grad(zcoord,grf_mixt_temp,ijkn)
#endif

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

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

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

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

          fd(1) = muel*(sij(1,1)*sface(1)+sij(1,2)*sface(2)+sij(1,3)*sface(3))
          fd(2) = muel*(sij(2,1)*sface(1)+sij(2,2)*sface(2)+sij(2,3)*sface(3))
          fd(3) = muel*(sij(3,1)*sface(1)+sij(3,2)*sface(2)+sij(3,3)*sface(3))
          fd(4) = dot_product(fd(1:3),velf(1:3)) + tcol*tgradf
!$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,ijkn)*modstrain
          trace(ijkc0) = trace(ijkc0) + taukk
          trace(ijkc1) = trace(ijkc1) + taukk

          muf    = -muet(ijk,ijkn)/muel
          tfd(1) = muf*fd(1)
          tfd(2) = muf*fd(2)
          tfd(3) = muf*fd(3)

          fd(1) = fd(1)-tfd(1)
          fd(2) = fd(2)-tfd(2)
          fd(3) = fd(3)-tfd(3)

!$startcopy TURB_VFluxHybrid
          fd4      = fd(4)
          ucitauij = dot_product(tfd(1:3),velc0(1:3))
          tfd(4)   = ucitauij - tcot*tgradf
          fd(4)    = fd4 - engmodel*tfd(4)

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

          ucitauij = dot_product(tfd(1:3),velc1(1:3))
          tfd(4)   = ucitauij - tcot*tgradf
          fd(4)    = fd4 - engmodel*tfd(4)

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

#ifdef RFLO
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k
#endif
#ifdef RFLU
    ENDDO       ! iC
#endif

!$endcopy

  END SUBROUTINE computefluxtot

END SUBROUTINE turb_visfluxeddy

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: TURB_VisFluxEddy.F90,v $
! Revision 1.11  2008/12/06 08:44:42  mtcampbe
! Updated license.
!
! Revision 1.10  2008/11/19 22:17:54  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.9  2005/12/30 23:20:48  wasistho
! exclude rocflu from WLM treatment
!
! Revision 1.8  2004/10/22 23:20:25  wasistho
! collect esg1 and esg2+3 into statistics
!
! Revision 1.7  2004/08/02 23:06:32  wasistho
! mv libfloflu/viscousFluxEddy(Patch) to rocflo/RFLO_viscousFlux(Patch)
!
! Revision 1.6  2004/08/02 19:34:15  wasistho
! changed grid%avgCo to grid%c2fCo
!
! Revision 1.5  2004/07/31 00:53:33  wasistho
! replaced cell2face midpoint by linear averaging
!
! Revision 1.4  2004/06/03 02:21:37  wasistho
! shift location of endcopy to be the same as VFluxHybrid
!
! Revision 1.3  2004/03/27 02:16:42  wasistho
! compiled with Rocflu
!
! 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:16  wasistho
! changed turb%esg.. to global%esg..
!
! Revision 1.6  2003/08/29 01:43:13  wasistho
! Added TARGET attribute to region variable, since pointers are cached into it
!
! Revision 1.5  2003/06/05 19:21:07  wasistho
! implemented heat transfer model
!
! Revision 1.4  2003/05/31 01:48:23  wasistho
! installed turb. wall layer model
!
! Revision 1.3  2003/05/15 02:57:06  jblazek
! Inlined index function.
!
! Revision 1.2  2003/02/04 03:59:49  wasistho
! An Esgs-bug which caused fd(4) unsymmetry fixed
!
! Revision 1.1  2002/10/14 23:55:30  wasistho
! Install Rocturb
!
!******************************************************************************