RFLO_ViscousFluxPatch.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 a patch
!
! Description: this routine works in the same way as ViscousFluxEddy
!              but applied on region patches
!
! Input: region  = data of current region.
!        patch   = current patch.
!        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: none.
!
!******************************************************************************
!
! $Id: RFLO_ViscousFluxPatch.F90,v 1.4 2008/12/06 08:44:28 mtcampbe Exp $
!
! Copyright: (c) 2001 by the University of Illinois
!
!******************************************************************************

SUBROUTINE rflo_viscousfluxpatch( region,patch,indxMu,indxTCo,tv )

  USE moddatatypes
  USE modbndpatch, ONLY   : t_patch
  USE moddatastruct, ONLY : t_region
  USE modinterfaces, ONLY : rflo_getpatchindices, rflo_getpatchdirection, &
                            rflo_getcelloffset, rflo_getnodeoffset

#include "Indexing.h"

  USE moderror
  USE modparameters
  IMPLICIT NONE

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

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

! ... local variables
  INTEGER :: ilev, lbound, bctype
  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend
  INTEGER :: inode, jnode, knode, idir, jdir, kdir
  INTEGER :: icoff, ijcoff, inoff, ijnoff, ijkcb0, ijkcd, ijknb, acid0, acid1
  INTEGER :: n1, n2, noff, j2d, aerocoeff

  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)          :: sgn, sf(3), sij(3,3), ac0, ac1
  REAL(RFREAL)          :: rref, vref, rcfref, rchref
  REAL(RFREAL), POINTER :: avgco(:,:), diss(:,:)

  REAL(RFREAL), POINTER :: dv(:,:), grad(:,:), sface(:,:)

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

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

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

  bctype = patch%bcType

  ilev   = region%currLevel
  lbound = patch%lbound
  
  dv   => region%levels(ilev)%mixt%dv
  diss => region%levels(ilev)%mixt%diss  

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

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

  noff      = abs(patch%l1end-patch%l1beg) + 1
  aerocoeff = region%global%aeroCoeffs
  rref      = region%global%refDensity
  vref      = region%global%refVelocity
  
  rcfref = 2.0_rfreal/(rref*vref*vref)
  rchref = 2.0_rfreal/(rref*vref*vref*vref)

! 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

! get the appropriate face vector --------------------------------------------

  IF (lbound==1 .OR. lbound==2) THEN
    avgco => region%levels(ilev)%grid%c2fCoI
    sface => region%levels(ilev)%grid%si
    grad  => region%levels(ilev)%mixt%gradi
  ELSE IF (lbound==3 .OR. lbound==4) THEN
    avgco => region%levels(ilev)%grid%c2fCoJ
    sface => region%levels(ilev)%grid%sj
    grad  => region%levels(ilev)%mixt%gradj
  ELSE
    avgco => region%levels(ilev)%grid%c2fCoK
    sface => region%levels(ilev)%grid%sk
    grad  => region%levels(ilev)%mixt%gradk
  ENDIF

! 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
          ac0    = avgco(acid0,ijknb)
          ac1    = avgco(acid1,ijknb)

          muf  = ac0*tv(tv_mixt_muel,ijkcb0)+ac1*tv(tv_mixt_muel,ijkcd)
          kpaf = ac0*tv(tv_mixt_tcol,ijkcb0)+ac1*tv(tv_mixt_tcol,ijkcd)

          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)
                             
          ux  = grad(gr_mixt_ux,ijknb)
          uy  = grad(gr_mixt_uy,ijknb)
          uz  = grad(gr_mixt_uz,ijknb)
          vx  = grad(gr_mixt_vx,ijknb)
          vy  = grad(gr_mixt_vy,ijknb)
          vz  = grad(gr_mixt_vz,ijknb)
          wx  = grad(gr_mixt_wx,ijknb)
          wy  = grad(gr_mixt_wy,ijknb)
          wz  = grad(gr_mixt_wz,ijknb)
          tx  = grad(gr_mixt_tx,ijknb)
          ty  = grad(gr_mixt_ty,ijknb)
          tz  = grad(gr_mixt_tz,ijknb)

          sf(1)  = sgn*sface(xcoord,ijknb)
          sf(2)  = sgn*sface(ycoord,ijknb)
          sf(3)  = sgn*sface(zcoord,ijknb)

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

          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

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

          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 

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

  CALL deregisterfunction( region%global )

END SUBROUTINE rflo_viscousfluxpatch

!******************************************************************************
!
! RCS Revision history:
! 
! $Log: RFLO_ViscousFluxPatch.F90,v $
! Revision 1.4  2008/12/06 08:44:28  mtcampbe
! Updated license.
!
! Revision 1.3  2008/11/19 22:17:39  mtcampbe
! Added Illinois Open Source License/Copyright
!
! Revision 1.2  2006/03/13 03:42:28  wasistho
! defined aero coeffs
!
! 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.15  2003/12/04 03:23:09  haselbac
! Removed RFLU code segments, moved to RFLU_ViscousFluxesPatches
!
! Revision 1.14  2003/11/20 16:40:36  mdbrandy
! Backing out RocfluidMP changes from 11-17-03
!
! Revision 1.11  2003/05/15 02:57:02  jblazek
! Inlined index function.
!
! Revision 1.10  2002/09/27 00:57:09  jblazek
! Changed makefiles - no makelinks needed.
!
! Revision 1.9  2002/09/20 22:22:35  jblazek
! Finalized integration into GenX.
!
! Revision 1.8  2002/09/09 14:08:12  haselbac
! mixtInput now under regions, added some things, corrected a few bugs
!
! Revision 1.7  2002/09/05 17:40:20  jblazek
! Variable global moved into regions().
!
! Revision 1.5  2002/08/01 01:40:03  wasistho
! Made generic
!
! Revision 1.4  2002/07/27 08:13:52  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:48  wasistho
! made compliant with CODING RULE
!
!******************************************************************************