TURB_RansSAVisFlux.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 SA viscous flux: nutot*d_j(tilde[nu])
!
! Description: this routine compute d_j(tilde[nu]), while nutot is given
!              nutot = nu_l+nu_t with l=laminar, t=turbulent
!
! Input: region  = data of current region
!
! Output: region%levels%turb%diss = viscous flux added to SA dissipation.
!
! Notes: none.
!
!******************************************************************************
!
! $Id: TURB_RansSAVisFlux.F90,v 1.12 2008/12/06 08:44:42 mtcampbe Exp $
!
! Copyright: (c) 2003 by the University of Illinois
!
!******************************************************************************

SUBROUTINE turb_ranssavisflux( region )

  USE moddatatypes
  USE modbndpatch, ONLY: t_patch
  USE moddatastruct, ONLY : t_region
#ifdef RFLO
  USE modinterfaces, ONLY : rflo_getdimensphys, rflo_getcelloffset, &
                            rflo_getnodeoffset, rflo_calcgradvector

#include "Indexing.h"
#endif
#ifdef RFLU
  USE rflu_moddifferentiationfaces, ONLY: rflu_computegradfaceswrapper
  USE rflu_moddifferentiationbfaces, ONLY: rflu_computegradbfaceswrapper
  USE modinterfaces, ONLY: rflu_decideneedbgradface
#endif
  USE turb_modinterfaces, ONLY : turb_ranssavisfluxpatch
  USE moderror
  USE modparameters
  USE turb_modparameters
  IMPLICIT NONE

! ... parameters
#ifdef RFLO
  TYPE(t_region) :: region
#endif
#ifdef RFLU
  TYPE(t_region), POINTER :: region
#endif

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

! ... local variables
  INTEGER :: ibegv, iendv, ibegg, iendg, ijkn, ijkc0, ijkc1
  REAL(RFREAL) :: cb2, opcb2, rsigma, beta, rhoa, mua, rnua
  REAL(RFREAL) :: nuf, nuc, nutilx, nutily, nutilz, fd
  REAL(RFREAL) :: sface(3)
  REAL(RFREAL), POINTER :: cv(:,:), tv(:,:), tcv(:,:), tdiss(:,:)
  REAL(RFREAL), POINTER :: vol(:)

#ifdef RFLO
  INTEGER :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER :: ilev, icoff, ijcoff, inoff, ijnoff
  REAL(RFREAL), POINTER :: avgco(:,:), sf(:,:), grad(:,:)
#endif
#ifdef RFLU
  INTEGER, POINTER      :: f2c(:,:)
  REAL(RFREAL), POINTER :: fn(:,:), grad(:,:,:)
  TYPE(t_patch), POINTER :: ppatch
#endif

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

  CALL registerfunction( region%global,'TURB_RansSAVisFlux',&
  'TURB_RansSAVisFlux.F90' )

! get dimensions and pointers ------------------------------------------------

#ifdef RFLO
  ilev   =  region%currLevel
  cv     => region%levels(ilev)%mixt%cv
  tv     => region%levels(ilev)%mixt%tv
  tcv    => region%levels(ilev)%turb%cv
  tdiss  => region%levels(ilev)%turb%diss
  vol    => region%levels(ilev)%grid%vol
  
  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 )
  ibegv = cv_sa_nutil
  iendv = cv_sa_nutil
  ibegg = gr_sa_nutilx
  iendg = gr_sa_nutilz
#endif
#ifdef RFLU
  cv     => region%mixt%cv
  tv     => region%mixt%tv
  tcv    => region%turb%cv
  tdiss  => region%turb%diss  
  vol    => region%grid%vol
  ibegv = cv_sa_nutil
  iendv = cv_sa_nutil
  ibegg = gr_sa_nutilx
  iendg = gr_sa_nutilx
#endif

! get needed quantities

  cb2    = region%turbInput%const(mc_sa_cb2)
  rsigma = region%turbInput%const(mc_sa_rsig)
  beta   = region%mixtInput%betrk(region%irkStep)*rsigma
  opcb2  = 1._rfreal+cb2

! get gradients of tilde[nu]

#ifdef RFLO
  CALL rflo_calcgradvector( region,ibegv,iendv,ibegg,iendg, &
                                   region%levels(ilev)%turb%cv, &
                                   region%levels(ilev)%turb%gradi, &
                                   region%levels(ilev)%turb%gradj, &
                                   region%levels(ilev)%turb%gradk )
#endif
#ifdef RFLU
  CALL rflu_computegradfaceswrapper( region,ibegv,iendv,ibegg,iendg, &
                                     region%turb%cv,region%turb%gradi )
  DO ipatch = 1,region%grid%nPatches
    ppatch => region%patches(ipatch)
                                  
    IF ( rflu_decideneedbgradface(region,ppatch) .EQV. .true. ) THEN
      CALL rflu_computegradbfaceswrapper( region,ppatch,ibegv,iendv,ibegg,iendg, &
                                          region%turb%cv,region%turb%bGradi )
    END IF ! RFLU_DecideNeedBGradFace
  END DO ! iPatch                                 
#endif

! interior fluxes -------------------------------------------------------------

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

! fluxes through boundaries ---------------------------------------------------
#ifdef RFLO  
  DO ipatch=1,region%nPatches
    CALL turb_ranssavisfluxpatch( region,region%levels(ilev)%patches(ipatch) )
  ENDDO
#endif
#ifdef RFLU
  DO ipatch = 1,region%grid%nPatches
    CALL turb_ranssavisfluxpatch( region,region%patches(ipatch) )
  END DO ! ipatch
#endif

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

  CALL deregisterfunction( region%global )

! =============================================================================
!   Flux computation subroutines
! =============================================================================

CONTAINS

  SUBROUTINE computeflux( ijk )

! ... parameters
    INTEGER   :: ijk

! ... local variables
    INTEGER   :: ibeg,iend,jbeg,jend,kbeg,kend, iadd,jadd,kadd
    REAL(RFREAL) :: ac0, ac1

! - 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
      grad  => region%levels(ilev)%turb%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
      grad  => region%levels(ilev)%turb%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
      grad  => region%levels(ilev)%turb%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
    grad => region%turb%gradi
    fn   => region%grid%fn
#endif

! -- flux in ijk-direction (except through boundary) --------------------------

#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)
          ac0   = avgco(2,ijkn)
          ac1   = avgco(1,ijkn)
          sface(1)= sf(xcoord,ijkn)
          sface(2)= sf(ycoord,ijkn)
          sface(3)= sf(zcoord,ijkn)
#endif
#ifdef RFLU
    ac0 = 0.5_rfreal
    ac1 = 0.5_rfreal
    DO ijkn = ibeg,iend
          ijkc0 = f2c(1,ijkn)
          ijkc1 = f2c(2,ijkn)
          sface(1)= fn(xcoord,ijkn)*fn(xyzmag,ijkn)
          sface(2)= fn(ycoord,ijkn)*fn(xyzmag,ijkn)
          sface(3)= fn(zcoord,ijkn)*fn(xyzmag,ijkn)           
#endif
          rhoa = ac0*cv(cv_mixt_dens ,ijkc0)+ac1*cv(cv_mixt_dens ,ijkc1)
          mua  = ac0*tv(tv_mixt_muel ,ijkc0)+ac1*tv(tv_mixt_muel ,ijkc1)
          rnua = ac0*tcv(cv_sa_nutil ,ijkc0)+ac1*tcv(cv_sa_nutil ,ijkc1)
          nuf  = (rnua + mua)/rhoa 
          nuc  = (tcv(cv_sa_nutil ,ijkc0) + tv(tv_mixt_muel ,ijkc0))/ &
                  cv(cv_mixt_dens ,ijkc0) 

#ifdef RFLO
          nutilx = grad(gr_sa_nutilx,ijkn)
          nutily = grad(gr_sa_nutily,ijkn)
          nutilz = grad(gr_sa_nutilz,ijkn)
#endif
#ifdef RFLU
          nutilx = grad(xcoord,gr_sa_nutilx,ijkn)
          nutily = grad(ycoord,gr_sa_nutilx,ijkn)
          nutilz = grad(zcoord,gr_sa_nutilx,ijkn)
#endif

          fd = beta*(nuf*opcb2-nuc*cb2)* &
               (nutilx*sface(1)+nutily*sface(2)+nutilz*sface(3))

          tdiss(cv_sa_nutil,ijkc0) = tdiss(cv_sa_nutil,ijkc0) + fd
          tdiss(cv_sa_nutil,ijkc1) = tdiss(cv_sa_nutil,ijkc1) - fd

#ifdef RFLO
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k
#else
    ENDDO       ! ijkN
#endif

  END SUBROUTINE computeflux

END SUBROUTINE turb_ranssavisflux

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: TURB_RansSAVisFlux.F90,v $
! Revision 1.12  2008/12/06 08:44:42  mtcampbe
! Updated license.
!
! Revision 1.11  2008/11/19 22:17:54  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.10  2006/08/19 15:40:38  mparmar
! Added use of RFLU_DecideNeedBGradFace
!
! Revision 1.9  2006/04/07 15:06:06  haselbac
! Bug fix: Incorrect ifs
!
! Revision 1.8  2006/04/07 14:56:02  haselbac
! Adapted to changes in f and bf grad routines
!
! Revision 1.7  2005/12/20 20:43:51  wasistho
! adapted to changing in Rocflu on face gradient routines
!
! Revision 1.6  2004/08/02 23:08:31  wasistho
! shift location of lines defining ac0 and ac1
!
! Revision 1.5  2004/08/02 21:55:39  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/24 03:37:02  wasistho
! prepared for RFLU
!
! Revision 1.2  2004/03/20 03:28:29  wasistho
! prepared for RFLU
!
! Revision 1.1  2004/03/05 04:37:00  wasistho
! changed nomenclature
!
! Revision 1.5  2003/10/25 22:07:26  wasistho
! modified non-conservative diffusion term
!
! Revision 1.4  2003/10/21 01:34:19  wasistho
! loop in computation of gradCell
!
! Revision 1.3  2003/10/20 20:28:21  wasistho
! made consistent with compressible SA formulation
!
! Revision 1.2  2003/10/10 20:35:06  wasistho
! multiplied SA viscous fluxes by 1/sigma through beta
!
! Revision 1.1  2003/10/07 02:17:03  wasistho
! initial installation of RaNS-SA and DES
!
!
!******************************************************************************