RFLO_ViscousFlux.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
!
! Description: this routine compute S_ij, while mu is given
!              mu can be mu_l or mu_tot=mu_l+mu_t with l=laminar, t=turbulent
!
! Input: region  = data of current region
!        indxMu  = TV component index for dynamic viscosity Mu 
!        indxTCo = TV component index for thermal conductivity TCo 
!        tv      = transport variables, mu and kappa
!
! Output: region%levels%mixt%diss = viscous fluxes added to the dissipation.
!
! Notes: indxMu can have value TV_MIXT_MUEL (laminar) or TV_MIXT_MUET (total)
!       indxTCo can have value TV_MIXT_TCOL (laminar) or TV_MIXT_TCOT (total)
!
!******************************************************************************
!
! $Id: RFLO_ViscousFlux.F90,v 1.3 2008/12/06 08:44:28 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE rflo_viscousflux( region,indxMu,indxTCo,tv )

  USE moddatatypes
  USE moddatastruct, ONLY : t_region
  USE modinterfaces, ONLY : rflo_getdimensphys, &
                            rflo_getcelloffset, rflo_getnodeoffset
  USE modinterfaces, ONLY : rflo_viscousfluxpatch

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

! ... parameters
  TYPE(t_region) :: region
  INTEGER        :: indxmu, indxtco
  REAL(RFREAL), POINTER :: tv(:,:)

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

! ... local variables
  INTEGER            :: ipcbeg, ipcend, jpcbeg, jpcend, kpcbeg, kpcend
  INTEGER            :: ilev, icoff, ijcoff, inoff, ijnoff, ijkc0, ijkc1, ijkn
  INTEGER, PARAMETER :: idir=1, jdir=2, kdir=3

  REAL(RFREAL)          :: oo3, beta, muf, kpaf, velf(3), div
  REAL(RFREAL)          :: ux, uy, uz, vx, vy, vz, wx, wy, wz, tx, ty, tz, &
                           tgradf, fd(4)
  REAL(RFREAL)          :: sface(3), sij(3,3)
  REAL(RFREAL), POINTER :: avgco(:,:), sf(:,:), diss(:,:)
  REAL(RFREAL), POINTER :: dv(:,:), grad(:,:)

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

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

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

  ilev =  region%currLevel
  dv   => region%levels(ilev)%mixt%dv
  diss => region%levels(ilev)%mixt%diss
  
  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 )

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

  oo3  = 1.0_rfreal/3.0_rfreal
  beta = region%mixtInput%betrk(region%irkStep)

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

  CALL computeflux( idir )
  CALL computeflux( jdir )
  CALL computeflux( kdir )

! fluxes through boundaries ---------------------------------------------------

  DO ipatch=1,region%nPatches
    CALL rflo_viscousfluxpatch( region,region%levels(ilev)%patches(ipatch), &
                                indxmu,indxtco,tv )
  ENDDO

! 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 ---------------------------------------------------

    IF (ijk==idir) THEN
      ibeg = ipcbeg+1
      iend = ipcend
      jbeg = jpcbeg
      jend = jpcend
      kbeg = kpcbeg
      kend = kpcend
      iadd = -1
      jadd = 0
      kadd = 0
      grad  => region%levels(ilev)%mixt%gradi
      sf    => region%levels(ilev)%grid%si
      avgco => region%levels(ilev)%grid%c2fCoI
    ELSEIF (ijk==jdir) THEN
      ibeg = ipcbeg
      iend = ipcend
      jbeg = jpcbeg+1
      jend = jpcend
      kbeg = kpcbeg
      kend = kpcend
      iadd = 0
      jadd = -1
      kadd = 0
      grad  => region%levels(ilev)%mixt%gradj
      sf    => region%levels(ilev)%grid%sj
      avgco => region%levels(ilev)%grid%c2fCoJ
    ELSEIF (ijk==kdir) THEN
      ibeg = ipcbeg
      iend = ipcend
      jbeg = jpcbeg
      jend = jpcend
      kbeg = kpcbeg+1
      kend = kpcend
      iadd = 0
      jadd = 0
      kadd = -1
      grad  => region%levels(ilev)%mixt%gradk
      sf    => region%levels(ilev)%grid%sk
      avgco => region%levels(ilev)%grid%c2fCoK
    ENDIF

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

    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)

          muf    = ac0*tv(indxmu      ,ijkc0)+ac1*tv(indxmu      ,ijkc1)
          kpaf   = ac0*tv(indxtco     ,ijkc0)+ac1*tv(indxtco     ,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)
          ux  = grad(gr_mixt_ux,ijkn)
          uy  = grad(gr_mixt_uy,ijkn)
          uz  = grad(gr_mixt_uz,ijkn)
          vx  = grad(gr_mixt_vx,ijkn)
          vy  = grad(gr_mixt_vy,ijkn)
          vz  = grad(gr_mixt_vz,ijkn)
          wx  = grad(gr_mixt_wx,ijkn)
          wy  = grad(gr_mixt_wy,ijkn)
          wz  = grad(gr_mixt_wz,ijkn)
          tx  = grad(gr_mixt_tx,ijkn)
          ty  = grad(gr_mixt_ty,ijkn)
          tz  = grad(gr_mixt_tz,ijkn)

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

          div = oo3*(ux+vy+wz)

          sij(1,1) = 2.0_rfreal*(ux-div)
          sij(1,2) = uy+vx
          sij(1,3) = uz+wx

          sij(2,1) = sij(1,2)
          sij(2,2) = 2.0_rfreal*(vy-div)
          sij(2,3) = vz+wy

          sij(3,1) = sij(1,3)
          sij(3,2) = sij(2,3)
          sij(3,3) = 2.0_rfreal*(wz-div)

          fd(1) = muf*(sij(1,1)*sface(1)+sij(1,2)*sface(2)+sij(1,3)*sface(3))
          fd(2) = muf*(sij(2,1)*sface(1)+sij(2,2)*sface(2)+sij(2,3)*sface(3))
          fd(3) = muf*(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)) + kpaf*tgradf

          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

          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
          
        ENDDO   ! i
      ENDDO     ! j
    ENDDO       ! k

  END SUBROUTINE computeflux

END SUBROUTINE rflo_viscousflux

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: RFLO_ViscousFlux.F90,v $
! Revision 1.3  2008/12/06 08:44:28  mtcampbe
! Updated license.
!
! Revision 1.2  2008/11/19 22:17:39  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.1  2004/11/29 20:51:40  wasistho
! lower to upper case
!
! Revision 1.1  2004/08/02 23:13:59  wasistho
! mv libfloflu/viscousFluxEddy(Patch) to rocflo/RFLO_viscousFlux(Patch)
!
! Revision 1.14  2003/12/04 03:23:07  haselbac
! Removed RFLU code segments, moved to RFLU_ViscousFluxes
!
! Revision 1.13  2003/11/20 16:40:35  mdbrandy
! Backing out RocfluidMP changes from 11-17-03
!
! Revision 1.10  2003/05/15 02:57:02  jblazek
! Inlined index function.
!
! Revision 1.9  2002/10/27 18:49:21  haselbac
! Removed tabs
!
! Revision 1.8  2002/09/20 22:22:35  jblazek
! Finalized integration into GenX.
!
! Revision 1.7  2002/09/09 14:07:23  haselbac
! mixtInput now under regions, added patch routine, bug fix for sf
!
! Revision 1.6  2002/09/05 17:40:20  jblazek
! Variable global moved into regions().
!
! Revision 1.5  2002/08/01 01:39:24  wasistho
! Included RFLU capability and made generic
!
! Revision 1.4  2002/07/27 08:13:37  wasistho
! prepared for rocturb implementation
!
! Revision 1.3  2002/07/22 22:59:11  jblazek
! Some more clean up.
!
! Revision 1.2  2002/07/19 23:38:42  wasistho
! made compliant with CODING RULE
!
!******************************************************************************